source: issm/trunk/src/m/classes/public/extrude.m@ 3811

Last change on this file since 3811 was 3811, checked in by Eric.Larour, 15 years ago

Added upper and lower elements neighboors in Penta element

File size: 8.2 KB
RevLine 
[2269]1function md=extrude(md,varargin)
[1]2%EXTRUDE - vertically extrude a 2d mesh
3%
[2269]4% vertically extrude a 2d mesh and create corresponding 3d mesh.
5% The vertical distribution can:
6% - follow a polynomial law
7% - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
8% - be discribed by a list of coefficients (between 0 and 1)
9%
[1]10%
11% Usage:
12% md=extrude(md,numlayers,extrusionexponent);
[2269]13% md=extrude(md,numlayers,lowerexponent,upperexponent);
14% md=extrude(md,listofcoefficients);
[1]15%
16% Example:
17% md=extrude(md,8,3);
[2269]18% md=extrude(md,8,3,2);
19% md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
[1]20%
21% See also: MODELEXTRACT, COLLAPSE
22
23%some checks on list of arguments
[2269]24if ((nargin>4) | (nargin<2) | (nargout~=1)),
[1241]25 help extrude;
[1]26 error('extrude error message');
27end
28
29if md.counter<3,
[1241]30 help extrude;
[1]31 error('only fully parameterized 2d models can be extruded');
32end
33
34if md.counter>=4,
35 error('This model has already been extruded!','s');
36end
37
[2269]38%Extrude the mesh
39if nargin==2, %list of coefficients
40 list=varargin{1};
41 if any(list<0) | any(list>1),
42 error('extrusioncoefficients must be between 0 and 1');
43 end
44 extrusionlist=sort(unique([list(:);0;1]));
45 numlayers=length(extrusionlist);
46elseif nargin==3, %one polynomial law
47 if varargin{2}<=0,
48 help extrude;
49 error('extrusionexponent must be >=0');
50 end
51 numlayers=varargin{1};
52 extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2};
53elseif nargin==4, %two polynomial laws
54 numlayers=varargin{1};
55 lowerexp=varargin{2};
56 upperexp=varargin{3};
57
58 if varargin{2}<=0 | varargin{3}<=0,
59 help extrude;
60 error('lower and upper extrusionexponents must be >=0');
61 end
62
63 lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
64 upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;
65 extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));
66
67end
68
[1]69if numlayers<2,
[1996]70 disp('number of layers should be at least 2. returning initial model...');
71 return
[1]72end
73
74%Initialize with the 2d mesh
[2269]75x3d=[];
76y3d=[];
77z3d=[]; %the lower grid is on the bed
[1]78thickness3d=md.thickness; %thickness and bed for these grids
79bed3d=md.bed;
[512]80
[1]81%Create the new layers
[2269]82for i=1:numlayers,
[512]83 x3d=[x3d; md.x];
[1]84 y3d=[y3d; md.y];
[512]85 %grids are distributed between bed and surface accordingly to the given exponent
[2269]86 z3d=[z3d; bed3d+thickness3d*extrusionlist(i)];
[1]87end
88number_grids3d=size(x3d,1); %number of 3d grids for the non extruded part of the mesh
89
90%Extrude elements
91elements3d=[];
92for i=1:numlayers-1,
93 elements3d=[elements3d;[md.elements+(i-1)*md.numberofgrids md.elements+i*md.numberofgrids]]; %Create the elements of the 3d mesh for the non extruded part
94end
95number_el3d=size(elements3d,1); %number of 3d grids for the non extruded part of the mesh
96
97%Keep a trace of lower and upper grids
98lowergrids=NaN*ones(number_grids3d,1);
99uppergrids=NaN*ones(number_grids3d,1);
100lowergrids(md.numberofgrids+1:end)=1:(numlayers-1)*md.numberofgrids;
101uppergrids(1:(numlayers-1)*md.numberofgrids)=md.numberofgrids+1:number_grids3d;
102md.lowergrids=lowergrids;
103md.uppergrids=uppergrids;
104
[3811]105%same for lower and upper elements
106lowerelements=NaN*ones(number_el3d,1);
107upperelements=NaN*ones(number_el3d,1);
108lowerelements(md.numberofelements+1:end)=1:(numlayers-2)*md.numberofelements;
109upperelements(1:(numlayers-2)*md.numberofelements)=md.numberofelements+1:(numlayers-1)*md.numberofelements;
110md.lowerelements=lowerelements;
111md.upperelements=upperelements;
112
113
114
[1]115%Save old mesh
116md.x2d=md.x;
117md.y2d=md.y;
118md.z2d=md.z;
119md.elements2d=md.elements;
120md.elements_type2d=md.elements_type;
121md.numberofelements2d=md.numberofelements;
122md.numberofgrids2d=md.numberofgrids;
123
124%Update mesh type
125md.type='3d';
126
127%Build global 3d mesh
128md.elements=elements3d;
129md.x=x3d;
130md.y=y3d;
131md.z=z3d;
132md.numberofelements=number_el3d;
133md.numberofgrids=number_grids3d;
134md.numlayers=numlayers;
135
136%Ok, now deal with the other fields from the 2d mesh:
137
[3760]138%drag_coefficient is limited to grids that are on the bedrock.
139md.drag_coefficient=project3d(md,md.drag_coefficient,'node',1);
[1]140
141%p and q (same deal, except for element that are on the bedrock: )
[3760]142md.drag_p=project3d(md,md.drag_p,'element');
143md.drag_q=project3d(md,md.drag_q,'element');
[1]144
145%observations
[510]146md.vx_obs=project3d(md,md.vx_obs,'node');
147md.vy_obs=project3d(md,md.vy_obs,'node');
148md.vel_obs=project3d(md,md.vel_obs,'node');
[2307]149md.vx_bal=project3d(md,md.vx_bal,'node');
150md.vy_bal=project3d(md,md.vy_bal,'node');
151md.vel_bal=project3d(md,md.vel_bal,'node');
[2252]152md.vel_obs_raw=project3d(md,md.vel_obs_raw,'node');
[3760]153md.accumulation_rate=project3d(md,md.accumulation_rate,'node');
[1]154md.firn_layer=project3d(md,md.firn_layer,'node',md.numlayers);
155
156%results
157if ~isnan(md.vx),md.vx=project3d(md,md.vx,'node');end;
158if ~isnan(md.vy),md.vy=project3d(md,md.vy,'node');end;
159if ~isnan(md.vz),md.vz=project3d(md,md.vz,'node');end;
160if ~isnan(md.vel),md.vel=project3d(md,md.vel,'node');end;
[1240]161if ~isnan(md.temperature),md.temperature=project3d(md,md.temperature,'node');end;
[1]162if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,md.surface_slopex,'node');end;
163if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,md.surface_slopey,'node');end;
164if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,md.bed_slopex,'node');end;
165if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,md.bed_slopey,'node');end;
166
167%bedinfo and surface info
168md.elementonbed=project3d(md,ones(md.numberofelements2d,1),'element',1);
169md.elementonsurface=project3d(md,ones(md.numberofelements2d,1),'element',md.numlayers-1);
170md.gridonbed=project3d(md,ones(md.numberofgrids2d,1),'node',1);
171md.gridonsurface=project3d(md,ones(md.numberofgrids2d,1),'node',md.numlayers);
172
173%elementstype
174if ~isnan(md.elements_type)
175 oldelements_type=md.elements_type2d;
176 md.elements_type=zeros(number_el3d,2);
177 md.elements_type(:,1)=project3d(md,oldelements_type(:,1),'element');
178 md.elements_type(:,2)=project3d(md,oldelements_type(:,2),'element');
179 md.gridonhutter=project3d(md,md.gridonhutter,'node');
180 md.gridonmacayeal=project3d(md,md.gridonmacayeal,'node');
181 md.gridonpattyn=project3d(md,md.gridonpattyn,'node');
182 md.gridonstokes=project3d(md,md.gridonstokes,'node');
[440]183
184 %dead grids
185 md.deadgrids=ones(md.numberofgrids,1);
[1651]186 md.deadgrids(md.elements(md.elements_type(:,1)~=MacAyealFormulationEnum,:))=0;%non macayeal grids are not dead
[440]187 md.deadgrids(find(md.gridonbed))=0;%grids from elements on bed are not dead
[1]188end
189
190%boundary conditions
[1760]191md.spcvelocity=project3d(md,md.spcvelocity,'node');
[1772]192md.spctemperature=project3d(md,md.spctemperature,'node',md.numlayers);
[1761]193md.spcthickness=project3d(md,md.spcthickness,'node');
[1]194
195%Extrusion of Neumann BC
196%in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element]
[3099]197pressureload_layer1=[md.pressureload(:,1:2) md.pressureload(:,2)+md.numberofgrids2d md.pressureload(:,1)+md.numberofgrids2d md.pressureload(:,3:4)]; %Add two columns on the first layer
[1755]198pressureload=[];
[1]199for i=1:numlayers-1,
[3099]200 pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.numberofgrids2d pressureload_layer1(:,5)+(i-1)*md.numberofelements2d pressureload_layer1(:,6)];
[1]201end
202
203%plug into md
[1755]204md.pressureload=pressureload;
[1]205
206%materials
[3760]207md.rheology_B=project3d(md,md.rheology_B,'node');
208md.rheology_n=project3d(md,md.rheology_n,'element');
[1]209
210%parameters
211md.surface=project3d(md,md.surface,'node');
212md.thickness=project3d(md,md.thickness,'node');
213md.bed=project3d(md,md.bed,'node');
214md.gridonboundary=project3d(md,md.gridonboundary,'node');
215md.elementoniceshelf=project3d(md,md.elementoniceshelf,'element');
216md.gridoniceshelf=project3d(md,md.gridoniceshelf,'node');
217md.elementonicesheet=project3d(md,md.elementonicesheet,'element');
218md.gridonicesheet=project3d(md,md.gridonicesheet,'node');
[1241]219md.elementonwater=project3d(md,md.elementonwater,'element');
220md.gridonwater=project3d(md,md.gridonwater,'node');
[3210]221if ~isnan(md.weights),md.weights=project3d(md,md.weights,'node');end;
[1]222
[2521]223%Put lithostatic pressure is there is an existing pressure
224if ~isnan(md.pressure),
225 md.pressure=md.g*md.rho_ice*(md.surface-md.z);
226end
227
[1]228%special for thermal modeling:
[3760]229md.melting_rate=project3d(md,md.melting_rate,'node',1);
[1]230md.observed_temperature=project3d(md,md.observed_temperature,'node');
231md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux
232
[867]233%increase connectivity if less than 25:
[2188]234if md.connectivity<=25,
235 md.connectivity=100;
[867]236end
237
[1]238%augment counter keeping track of what has been done to this model
239md.counter=4;
Note: See TracBrowser for help on using the repository browser.