0001 function md=extrude(md,numlayers,extrusionexponent)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 if ((nargin~=3) | (nargout~=1)),
0020 extrudeusage();
0021 error('extrude error message');
0022 end
0023
0024 if md.counter<3,
0025 extrudeusage();
0026 error('only fully parameterized 2d models can be extruded');
0027 end
0028
0029 if md.counter>=4,
0030 error('This model has already been extruded!','s');
0031 end
0032
0033 if numlayers<2,
0034 extrudeusage();
0035 error('number of layers should be at least 2');
0036 end
0037
0038 if extrusionexponent<0,
0039 extrudeusage();
0040 error('extrusionexponent must be >0');
0041 end
0042
0043
0044
0045 x3d=md.x;
0046 y3d=md.y;
0047 z3d=md.bed;
0048 thickness3d=md.thickness;
0049 bed3d=md.bed;
0050
0051 for i=2:numlayers,
0052 x3d=[x3d; md.x];
0053 y3d=[y3d; md.y];
0054 z3d=[z3d; bed3d+thickness3d*(i/numlayers)^extrusionexponent];
0055 end
0056 number_grids3d=size(x3d,1);
0057
0058
0059 elements3d=[];
0060 for i=1:numlayers-1,
0061 elements3d=[elements3d;[md.elements+(i-1)*md.numberofgrids md.elements+i*md.numberofgrids]];
0062 end
0063 number_el3d=size(elements3d,1);
0064
0065
0066 lowergrids=NaN*ones(number_grids3d,1);
0067 uppergrids=NaN*ones(number_grids3d,1);
0068 lowergrids(md.numberofgrids+1:end)=1:(numlayers-1)*md.numberofgrids;
0069 uppergrids(1:(numlayers-1)*md.numberofgrids)=md.numberofgrids+1:number_grids3d;
0070 md.lowergrids=lowergrids;
0071 md.uppergrids=uppergrids;
0072
0073
0074 md.x2d=md.x;
0075 md.y2d=md.y;
0076 md.z2d=md.z;
0077 md.elements2d=md.elements;
0078 md.elements_type2d=md.elements_type;
0079 md.numberofelements2d=md.numberofelements;
0080 md.numberofgrids2d=md.numberofgrids;
0081
0082
0083 md.type='3d';
0084
0085
0086 md.elements=elements3d;
0087 md.x=x3d;
0088 md.y=y3d;
0089 md.z=z3d;
0090 md.numberofelements=number_el3d;
0091 md.numberofgrids=number_grids3d;
0092 md.numlayers=numlayers;
0093
0094
0095
0096
0097 md.drag=project3d(md,md.drag,'node',1);
0098
0099
0100 md.p=project3d(md,md.p,'element');
0101 md.q=project3d(md,md.q,'element');
0102
0103
0104 md.vx_obs=project3d(md,md.vx_obs,'node',md.numlayers);
0105 md.vy_obs=project3d(md,md.vy_obs,'node',md.numlayers);
0106 md.vel_obs=project3d(md,md.vel_obs,'node',md.numlayers);
0107 md.accumulation=project3d(md,md.accumulation,'node',md.numlayers);
0108 md.firn_layer=project3d(md,md.firn_layer,'node',md.numlayers);
0109
0110
0111 if ~isnan(md.vx),md.vx=project3d(md,md.vx,'node');end;
0112 if ~isnan(md.vy),md.vy=project3d(md,md.vy,'node');end;
0113 if ~isnan(md.vz),md.vz=project3d(md,md.vz,'node');end;
0114 if ~isnan(md.vel),md.vel=project3d(md,md.vel,'node');end;
0115 if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,md.surface_slopex,'node');end;
0116 if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,md.surface_slopey,'node');end;
0117 if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,md.bed_slopex,'node');end;
0118 if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,md.bed_slopey,'node');end;
0119
0120
0121 md.elementonbed=project3d(md,ones(md.numberofelements2d,1),'element',1);
0122 md.elementonsurface=project3d(md,ones(md.numberofelements2d,1),'element',md.numlayers-1);
0123 md.gridonbed=project3d(md,ones(md.numberofgrids2d,1),'node',1);
0124 md.gridonsurface=project3d(md,ones(md.numberofgrids2d,1),'node',md.numlayers);
0125
0126
0127 if ~isnan(md.elements_type)
0128 oldelements_type=md.elements_type2d;
0129 md.elements_type=zeros(number_el3d,2);
0130 md.elements_type(:,1)=project3d(md,oldelements_type(:,1),'element');
0131 md.elements_type(:,2)=project3d(md,oldelements_type(:,2),'element');
0132 md.gridonhutter=project3d(md,md.gridonhutter,'node');
0133 md.gridonmacayeal=project3d(md,md.gridonmacayeal,'node');
0134 md.gridonpattyn=project3d(md,md.gridonpattyn,'node');
0135 md.gridonstokes=project3d(md,md.gridonstokes,'node');
0136 end
0137
0138
0139 md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node');
0140 md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node');
0141
0142
0143
0144 oldsegmentonneumann_diag=md.segmentonneumann_diag;
0145 segmentonneumann_diag_layer1=[oldsegmentonneumann_diag(:,1:2) oldsegmentonneumann_diag(:,2)+md.numberofgrids2d oldsegmentonneumann_diag(:,1)+md.numberofgrids2d oldsegmentonneumann_diag(:,3)];
0146 segmentonneumann_diag=[];
0147 for i=1:numlayers-1,
0148 segmentonneumann_diag=[segmentonneumann_diag ;segmentonneumann_diag_layer1(:,1:4)+(i-1)*md.numberofgrids2d segmentonneumann_diag_layer1(:,5)+(i-1)*md.numberofelements2d ];
0149 end
0150
0151
0152 md.segmentonneumann_diag=segmentonneumann_diag;
0153 md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node');
0154 md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node',1);
0155
0156
0157
0158
0159 md.B=project3d(md,md.B,'node');
0160 md.n=project3d(md,md.n,'element');
0161
0162
0163 md.surface=project3d(md,md.surface,'node');
0164 md.thickness=project3d(md,md.thickness,'node');
0165 md.bed=project3d(md,md.bed,'node');
0166 md.gridonboundary=project3d(md,md.gridonboundary,'node');
0167 md.elementoniceshelf=project3d(md,md.elementoniceshelf,'element');
0168 md.gridoniceshelf=project3d(md,md.gridoniceshelf,'node');
0169 md.elementonicesheet=project3d(md,md.elementonicesheet,'element');
0170 md.gridonicesheet=project3d(md,md.gridonicesheet,'node');
0171
0172
0173 md.melting=project3d(md,md.melting,'node',1);
0174 md.observed_temperature=project3d(md,md.observed_temperature,'node');
0175 md.geothermalflux=project3d(md,md.geothermalflux,'node',1);
0176 md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers);
0177 md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers);
0178
0179
0180 pos=find(~md.gridondirichlet_thermal);
0181 md.dirichletvalues_thermal(pos)=NaN;
0182
0183
0184 md.connectivity=25;
0185
0186
0187 md.counter=4;
0188
0189 end
0190
0191 function extrudeusage(),
0192 disp('INPUT function md=extrude(md,numlayers,extrusionexponent) ');
0193 disp(' vertically extrude a 2d mesh and create corresponding 3d mesh, ');
0194 disp(' the vertical distribution follows a polynomial law:');
0195 disp(' extrusionexponent>1 refinement at the base');
0196 disp(' 0<extrusionexponent<1 refinement at the surface');
0197 disp(' extrusionexponent=1 linear extrusion');
0198 end