setelementstype

PURPOSE ^

SETELEMENTSTYPE - associate a solution type to each element

SYNOPSIS ^

function md=setelementstype(md,varargin)

DESCRIPTION ^

SETELEMENTSTYPE - associate a solution type to each element

   This routine works like plotmodel: it works with an even number of inputs
   'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options
   that must be followed by the corresponding exp file or flags list
   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
   If user wants every element outside the domain to be 
   setelementstyped, add '~' to the name of the domain file (ex: '~Pattyn.exp');
   an empty string '' will be considered as an empty domain
   a string 'all' will be considered as the entire domain

   Usage:
      md=setelementstype(md,varargin)

   Example:
      md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter');

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=setelementstype(md,varargin)
0002 %SETELEMENTSTYPE - associate a solution type to each element
0003 %
0004 %   This routine works like plotmodel: it works with an even number of inputs
0005 %   'hutter','macayeal','pattyn','stokes' and 'fill' are the possible options
0006 %   that must be followed by the corresponding exp file or flags list
0007 %   It can either be a domain file (argus type, .exp extension), or an array of element flags.
0008 %   If user wants every element outside the domain to be
0009 %   setelementstyped, add '~' to the name of the domain file (ex: '~Pattyn.exp');
0010 %   an empty string '' will be considered as an empty domain
0011 %   a string 'all' will be considered as the entire domain
0012 %
0013 %   Usage:
0014 %      md=setelementstype(md,varargin)
0015 %
0016 %   Example:
0017 %      md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter');
0018 
0019 %some checks on list of arguments
0020 if ((nargin<2) | (nargout~=1)),
0021     error('setelementstype error message');
0022 end
0023 
0024 if md.counter<2,
0025     error('only fully parameterized 2d models can be setelementstyped');
0026 end
0027 
0028 [hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
0029 
0030 %Flag the elements that has not been flagged as filltype
0031 if strcmpi(filltype,'hutter'),
0032     hutterflag(find(~macayealflag & ~pattynflag))=1;
0033 elseif strcmpi(filltype,'macayeal'),
0034     macayealflag(find(~hutterflag & ~pattynflag))=1;
0035 elseif strcmpi(filltype,'pattyn'),
0036     pattynflag(find(~hutterflag & ~macayealflag))=1;
0037 end
0038 
0039 %check that each element has at least one flag
0040 if any(hutterflag+ macayealflag+pattynflag==0),
0041     error('setelementstype error message: elements type not assigned, must be specified')
0042 end
0043 
0044 %check that each element has only one flag
0045 if any(hutterflag+ macayealflag+pattynflag>1),
0046     disp('setelementstype warning message: some elements have several types, higher order type is used for them')
0047     hutterflag(find(hutterflag & macayealflag))=0;
0048     macayealflag(find(macayealflag & pattynflag))=0;
0049     pattynflag(find(pattynflag & macayealflag))=0;
0050 end
0051 
0052 %Check that no pattyn or stokes for 2d mesh
0053 if strcmpi(md.type,'2d'),
0054     if any(stokesflag | pattynflag)
0055         error('setelementstype error message: stokes and pattyn elements no allowed in 2d mesh, extrude it first')
0056     end
0057 end
0058 
0059 %Modify stokesflag to get rid of elements contrained everywhere
0060 nonstokesflag=find(~stokesflag);                   %non stokes elements
0061 gridonstokes=ones(md.numberofgrids,1);
0062 gridonstokes(md.elements(nonstokesflag,:))=0;      %non stokes grids
0063 pos=find(md.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
0064 gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront
0065 
0066 %Rule out elements where all grids are on the border of the stokes grids (because those elements are completely
0067 %constrained
0068 pos=find(stokesflag);
0069 constrained_list=pos(find(sum(gridonstokes(md.elements(pos,:)),2)==0));
0070 stokesflag(constrained_list)=0;
0071 
0072 %add in model who is who
0073 md.elements_type=zeros(md.numberofelements,2);
0074 
0075 %1: Hutter elements
0076 gridonhutter=zeros(md.numberofgrids,1);
0077 gridonhutter(md.elements(find(hutterflag),:))=1;
0078 md.gridonhutter=gridonhutter;
0079 md.elements_type(find(hutterflag),1)=hutterenum();
0080 
0081 %2: MacAyeal elements
0082 gridonmacayeal=zeros(md.numberofgrids,1);
0083 gridonmacayeal(md.elements(find(macayealflag),:))=1;
0084 md.gridonmacayeal=gridonmacayeal;
0085 md.elements_type(find(macayealflag),1)=macayealenum();
0086 
0087 %3: Pattyn elements
0088 gridonpattyn=zeros(md.numberofgrids,1);
0089 gridonpattyn(md.elements(find(pattynflag),:))=1;
0090 md.gridonpattyn=gridonpattyn;
0091 md.elements_type(find(pattynflag),1)=pattynenum();
0092 
0093 %4: Stokes elements
0094 md.gridonstokes=gridonstokes;
0095 md.elements_type(find(stokesflag),2)=stokesenum();
0096 
0097 %5: None elements (non Stokes)
0098 md.elements_type(find(~stokesflag),2)=noneenum();
0099 
0100 %Create the border grids between Pattyn and MacAyeal and extrude them
0101 numgrids2d=md.numberofgrids2d;
0102 numlayers=md.numlayers;
0103 bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements
0104 
0105 %initialize and fill in penalties structure
0106 md.penalties=[];
0107 if ~isnan(bordergrids2d),
0108     penalties=[];
0109     for    i=1:numlayers,
0110         penalties=[penalties [bordergrids2d+md.numberofgrids2d*(i-1)]];
0111     end
0112     md.penalties=penalties;
0113 end
0114 
0115 %flag dead grids (strictly in MacAyeal and not on bed -> not used in diagnostic horiz)
0116 nonmacayeal_el=find(~macayealflag);                   %non macayeal elements
0117 deadgrids=ones(md.numberofgrids,1);
0118 deadgrids(md.elements(nonmacayeal_el,:))=0;          %non macayeal grids are not dead
0119 deadgrids(find(md.gridonbed))=0;                      %grid on bed are not dead
0120 md.deadgrids=deadgrids;
0121 end

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003