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

function md=extrude(md,varargin)
%EXTRUDE - vertically extrude a 2d mesh
%
%   vertically extrude a 2d mesh and create corresponding 3d mesh.
%   The vertical distribution can:
%    - follow a polynomial law
%    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
%    - be discribed by a list of coefficients (between 0 and 1)
%   
%
%   Usage:
%      md=extrude(md,numlayers,extrusionexponent);
%      md=extrude(md,numlayers,lowerexponent,upperexponent);
%      md=extrude(md,listofcoefficients);
%
%   Example:
%      md=extrude(md,8,3);
%      md=extrude(md,8,3,2);
%      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
%
%   See also: MODELEXTRACT, COLLAPSE

%some checks on list of arguments
if ((nargin>4) | (nargin<2) | (nargout~=1)),
        help extrude;
        error('extrude error message');
end

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

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

%Extrude the mesh
if nargin==2, %list of coefficients
        list=varargin{1};
        if any(list<0) | any(list>1),
                error('extrusioncoefficients must be between 0 and 1');
        end
        extrusionlist=sort(unique([list(:);0;1]));
        numlayers=length(extrusionlist);
elseif nargin==3, %one polynomial law
        if varargin{2}<=0,
                help extrude;
                error('extrusionexponent must be >=0');
        end
        numlayers=varargin{1};
        extrusionlist=((0:1:numlayers-1)/(numlayers-1)).^varargin{2};
elseif nargin==4, %two polynomial laws
        numlayers=varargin{1};
        lowerexp=varargin{2};
        upperexp=varargin{3};

        if varargin{2}<=0 | varargin{3}<=0,
                help extrude;
                error('lower and upper extrusionexponents must be >=0');
        end

        lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
        upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;
        extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));

end

if numlayers<2,
        disp('number of layers should be at least 2. returning initial model...');
        return
end

%Initialize with the 2d mesh
x3d=[]; 
y3d=[];
z3d=[];  %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=1:numlayers,
        x3d=[x3d; md.x]; 
        y3d=[y3d; md.y];
        %grids are distributed between bed and surface accordingly to the given exponent
        z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; 
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;

%same for lower and upper elements
lowerelements=NaN*ones(number_el3d,1);
upperelements=NaN*ones(number_el3d,1);
lowerelements(md.numberofelements+1:end)=1:(numlayers-2)*md.numberofelements;
upperelements(1:(numlayers-2)*md.numberofelements)=md.numberofelements+1:(numlayers-1)*md.numberofelements;
md.lowerelements=lowerelements;
md.upperelements=upperelements;



%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_coefficient is limited to grids that are on the bedrock.
md.drag_coefficient=project3d(md,md.drag_coefficient,'node',1);

%p and q (same deal, except for element that are on the bedrock: )
md.drag_p=project3d(md,md.drag_p,'element');
md.drag_q=project3d(md,md.drag_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.vx_bal=project3d(md,md.vx_bal,'node');
md.vy_bal=project3d(md,md.vy_bal,'node');
md.vel_bal=project3d(md,md.vel_bal,'node');
md.vel_obs_raw=project3d(md,md.vel_obs_raw,'node');
md.accumulation_rate=project3d(md,md.accumulation_rate,'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.temperature),md.temperature=project3d(md,md.temperature,'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)~=MacAyealFormulationEnum,:))=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.spcvelocity=project3d(md,md.spcvelocity,'node');
md.spctemperature=project3d(md,md.spctemperature,'node',md.numlayers);
md.spcthickness=project3d(md,md.spcthickness,'node');

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

%plug into md
md.pressureload=pressureload;

%materials
md.rheology_B=project3d(md,md.rheology_B,'node');
md.rheology_n=project3d(md,md.rheology_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');
md.elementonwater=project3d(md,md.elementonwater,'element');
md.gridonwater=project3d(md,md.gridonwater,'node');
if ~isnan(md.weights),md.weights=project3d(md,md.weights,'node');end;

%Put lithostatic pressure is there is an existing pressure
if ~isnan(md.pressure),
        md.pressure=md.g*md.rho_ice*(md.surface-md.z);
end

%special for thermal modeling:
md.melting_rate=project3d(md,md.melting_rate,'node',1); 
md.observed_temperature=project3d(md,md.observed_temperature,'node'); 
md.geothermalflux=project3d(md,md.geothermalflux,'node',1); %bedrock only gets geothermal flux

%increase connectivity if less than 25:
if md.connectivity<=25,
        md.connectivity=100;
end

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