extrude

PURPOSE ^

EXTRUDE - vertically extrude a 2d mesh

SYNOPSIS ^

function md=extrude(md,numlayers,extrusionexponent)

DESCRIPTION ^

EXTRUDE - vertically extrude a 2d mesh

   vertically extrude a 2d mesh and create corresponding 3d mesh following the extrusionexponent, 
   the vertical distribution follows a polynomial law:');
     extrusionexponent>1    refinement at the base');
   0<extrusionexponent<1    refinement at the surface');
     extrusionexponent=1    linear extrusion');

   Usage:
      md=extrude(md,numlayers,extrusionexponent);

   Example:
      md=extrude(md,8,3);

   See also: MODELEXTRACT, COLLAPSE

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function md=extrude(md,numlayers,extrusionexponent)
0002 %EXTRUDE - vertically extrude a 2d mesh
0003 %
0004 %   vertically extrude a 2d mesh and create corresponding 3d mesh following the extrusionexponent,
0005 %   the vertical distribution follows a polynomial law:');
0006 %     extrusionexponent>1    refinement at the base');
0007 %   0<extrusionexponent<1    refinement at the surface');
0008 %     extrusionexponent=1    linear extrusion');
0009 %
0010 %   Usage:
0011 %      md=extrude(md,numlayers,extrusionexponent);
0012 %
0013 %   Example:
0014 %      md=extrude(md,8,3);
0015 %
0016 %   See also: MODELEXTRACT, COLLAPSE
0017 
0018 %some checks on list of arguments
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 %Extrude the mesh
0044 %Initialize with the 2d mesh
0045 x3d=md.x; 
0046 y3d=md.y;
0047 z3d=md.bed;  %the lower grid is on the bed
0048 thickness3d=md.thickness; %thickness and bed for these grids
0049 bed3d=md.bed;
0050 %Create the new layers
0051 for i=2:numlayers,
0052     x3d=[x3d; md.x]; %build the grids of the other layers
0053     y3d=[y3d; md.y];
0054     z3d=[z3d; bed3d+thickness3d*(i/numlayers)^extrusionexponent]; %grids are distributed between bed and surface accordingly to the given exponent
0055 end
0056 number_grids3d=size(x3d,1); %number of 3d grids for the non extruded part of the mesh
0057 
0058 %Extrude elements
0059 elements3d=[];
0060 for i=1:numlayers-1,
0061     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
0062 end
0063 number_el3d=size(elements3d,1); %number of 3d grids for the non extruded part of the mesh
0064 
0065 %Keep a trace of lower and upper grids
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 %Save old mesh
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 %Update mesh type
0083 md.type='3d';
0084 
0085 %Build global 3d mesh
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 %Ok, now deal with the other fields from the 2d mesh:
0095 
0096 %drag is limited to grids that are on the bedrock.
0097 md.drag=project3d(md,md.drag,'node',1);
0098 
0099 %p and q (same deal, except for element that are on the bedrock: )
0100 md.p=project3d(md,md.p,'element');
0101 md.q=project3d(md,md.q,'element');
0102 
0103 %observations
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 %results
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 %bedinfo and surface info
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 %elementstype
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 %boundary conditions
0139 md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node');
0140 md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node'); 
0141 
0142 %Extrusion of Neumann BC
0143 %in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element]
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)]; %Add two columns on the first layer
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 %plug into md
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 %md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))];
0156 %md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))];
0157 
0158 %materials
0159 md.B=project3d(md,md.B,'node');
0160 md.n=project3d(md,md.n,'element');
0161 
0162 %parameters
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 %special for thermal modeling:
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); %bedrock only gets geothermal flux
0176 md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers); %surface temperature
0177 md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers); %surface temperature
0178 
0179 %NaN the values that are not on an spc'd temperature.
0180 pos=find(~md.gridondirichlet_thermal);
0181 md.dirichletvalues_thermal(pos)=NaN;
0182 
0183 %update of connectivity matrix
0184 md.connectivity=25;
0185 
0186 %augment counter  keeping track of what has been done to this model
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

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