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

function md=extrude(md,numlayers,extrusionexponent)
%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

%some checks on list of arguments
if ((nargin~=3) | (nargout~=1)),
        extrudeusage();
        error('extrude error message');
end

if md.counter<3,
        extrudeusage();
        error('only fully parameterized 2d models can be extruded');
end

if md.counter>=4,
        error('This model has already been extruded!','s');
end

if numlayers<2,
        extrudeusage();
        error('number of layers should be at least 2');
end

if extrusionexponent<0,
        extrudeusage();
        error('extrusionexponent must be >0');
end

%Extrude the mesh
%Initialize with the 2d mesh
x3d=md.x; 
y3d=md.y;
z3d=md.bed;  %the lower grid is on the bed
thickness3d=md.thickness; %thickness and bed for these grids
bed3d=md.bed;
%Create the new layers
for i=2:numlayers,
        x3d=[x3d; md.x]; %build the grids of the other layers
        y3d=[y3d; md.y];
        z3d=[z3d; bed3d+thickness3d*(i/numlayers)^extrusionexponent]; %grids are distributed between bed and surface accordingly to the given exponent
end
number_grids3d=size(x3d,1); %number of 3d grids for the non extruded part of the mesh

%Extrude elements 
elements3d=[];
for i=1:numlayers-1,
        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
end
number_el3d=size(elements3d,1); %number of 3d grids for the non extruded part of the mesh

%Keep a trace of lower and upper grids
lowergrids=NaN*ones(number_grids3d,1);
uppergrids=NaN*ones(number_grids3d,1);
lowergrids(md.numberofgrids+1:end)=1:(numlayers-1)*md.numberofgrids;
uppergrids(1:(numlayers-1)*md.numberofgrids)=md.numberofgrids+1:number_grids3d;
md.lowergrids=lowergrids;
md.uppergrids=uppergrids;

%Save old mesh 
md.x2d=md.x;
md.y2d=md.y;
md.z2d=md.z;
md.elements2d=md.elements;
md.elements_type2d=md.elements_type;
md.numberofelements2d=md.numberofelements;
md.numberofgrids2d=md.numberofgrids;

%Update mesh type
md.type='3d';

%Build global 3d mesh 
md.elements=elements3d;
md.x=x3d;
md.y=y3d;
md.z=z3d;
md.numberofelements=number_el3d;
md.numberofgrids=number_grids3d;
md.numlayers=numlayers;

%Ok, now deal with the other fields from the 2d mesh:

%drag is limited to grids that are on the bedrock.
md.drag=project3d(md,md.drag,'node',1);

%p and q (same deal, except for element that are on the bedrock: )
md.p=project3d(md,md.p,'element');
md.q=project3d(md,md.q,'element');

%observations
md.vx_obs=project3d(md,md.vx_obs,'node');
md.vy_obs=project3d(md,md.vy_obs,'node');
md.vel_obs=project3d(md,md.vel_obs,'node');
md.accumulation=project3d(md,md.accumulation,'node');
md.firn_layer=project3d(md,md.firn_layer,'node',md.numlayers);

%results
if ~isnan(md.vx),md.vx=project3d(md,md.vx,'node');end;
if ~isnan(md.vy),md.vy=project3d(md,md.vy,'node');end;
if ~isnan(md.vz),md.vz=project3d(md,md.vz,'node');end;
if ~isnan(md.vel),md.vel=project3d(md,md.vel,'node');end;
if ~isnan(md.surface_slopex),md.surface_slopex=project3d(md,md.surface_slopex,'node');end;
if ~isnan(md.surface_slopey),md.surface_slopey=project3d(md,md.surface_slopey,'node');end;
if ~isnan(md.bed_slopex),md.bed_slopex=project3d(md,md.bed_slopex,'node');end;
if ~isnan(md.bed_slopey),md.bed_slopey=project3d(md,md.bed_slopey,'node');end;

%bedinfo and surface info
md.elementonbed=project3d(md,ones(md.numberofelements2d,1),'element',1);
md.elementonsurface=project3d(md,ones(md.numberofelements2d,1),'element',md.numlayers-1);
md.gridonbed=project3d(md,ones(md.numberofgrids2d,1),'node',1);
md.gridonsurface=project3d(md,ones(md.numberofgrids2d,1),'node',md.numlayers);

%elementstype
if ~isnan(md.elements_type)
        oldelements_type=md.elements_type2d;
        md.elements_type=zeros(number_el3d,2);
        md.elements_type(:,1)=project3d(md,oldelements_type(:,1),'element');
        md.elements_type(:,2)=project3d(md,oldelements_type(:,2),'element');
        md.gridonhutter=project3d(md,md.gridonhutter,'node');
        md.gridonmacayeal=project3d(md,md.gridonmacayeal,'node');
        md.gridonpattyn=project3d(md,md.gridonpattyn,'node');
        md.gridonstokes=project3d(md,md.gridonstokes,'node');

        %dead grids
        md.deadgrids=ones(md.numberofgrids,1);
        md.deadgrids(md.elements(md.elements_type(:,1)~=macayealenum,:))=0;%non macayeal grids are not dead
        md.deadgrids(find(md.gridonbed))=0;%grids from elements on bed are not dead
end

%boundary conditions
md.gridondirichlet_diag=project3d(md,md.gridondirichlet_diag,'node');
md.dirichletvalues_diag=project3d(md,md.dirichletvalues_diag,'node'); 

%Extrusion of Neumann BC
%in 3d, segmentonnumann is: [grid1 grid2 grid3 grid4 element]
oldsegmentonneumann_diag=md.segmentonneumann_diag;
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 
segmentonneumann_diag=[];
for i=1:numlayers-1,
        segmentonneumann_diag=[segmentonneumann_diag ;segmentonneumann_diag_layer1(:,1:4)+(i-1)*md.numberofgrids2d segmentonneumann_diag_layer1(:,5)+(i-1)*md.numberofelements2d ];
end

%plug into md
md.segmentonneumann_diag=segmentonneumann_diag;
md.gridondirichlet_prog=project3d(md,md.gridondirichlet_prog,'node');
md.dirichletvalues_prog=project3d(md,md.dirichletvalues_prog,'node');
%md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,3))];
%md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,3))];

%materials
md.B=project3d(md,md.B,'node');
md.n=project3d(md,md.n,'element');

%parameters
md.surface=project3d(md,md.surface,'node');
md.thickness=project3d(md,md.thickness,'node');
md.bed=project3d(md,md.bed,'node');
md.gridonboundary=project3d(md,md.gridonboundary,'node');
md.elementoniceshelf=project3d(md,md.elementoniceshelf,'element');
md.gridoniceshelf=project3d(md,md.gridoniceshelf,'node');
md.elementonicesheet=project3d(md,md.elementonicesheet,'element');
md.gridonicesheet=project3d(md,md.gridonicesheet,'node');

%special for thermal modeling:
md.melting=project3d(md,md.melting,'node',1); 
md.observed_temperature=project3d(md,md.observed_temperature,'node'); 
md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux
md.gridondirichlet_thermal=project3d(md,md.gridondirichlet_thermal,'node',md.numlayers); %surface temperature
md.dirichletvalues_thermal=project3d(md,md.dirichletvalues_thermal,'node',md.numlayers); %surface temperature

%NaN the values that are not on an spc'd temperature.
pos=find(~md.gridondirichlet_thermal);
md.dirichletvalues_thermal(pos)=NaN;

%update of connectivity matrix
md.connectivity=25;

%augment counter  keeping track of what has been done to this model
md.counter=4;

end

function extrudeusage(),
disp('INPUT function md=extrude(md,numlayers,extrusionexponent) ');
disp('     vertically extrude a 2d mesh and create corresponding 3d mesh,  ');
disp('     the vertical distribution follows a polynomial law:');
disp('          extrusionexponent>1  refinement at the base');
disp('        0<extrusionexponent<1  refinement at the surface');
disp('          extrusionexponent=1  linear extrusion');
end