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)),
	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

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

if extrusionexponent<=0,
	help extrude;
	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=1:numlayers-1,
	x3d=[x3d; md.x]; 
	y3d=[y3d; md.y];
	%grids are distributed between bed and surface accordingly to the given exponent
	z3d=[z3d; bed3d+thickness3d*(i/(numlayers-1))^extrusionexponent]; 
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.vel_obs_raw=project3d(md,md.vel_obs_raw,'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.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]
oldpressureload=md.pressureload;
pressureload_layer1=[oldpressureload(:,1:2)  oldpressureload(:,2)+md.numberofgrids2d  oldpressureload(:,1)+md.numberofgrids2d  oldpressureload(:,3)]; %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 ];
end

%plug into md
md.pressureload=pressureload;

%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');
md.elementonwater=project3d(md,md.elementonwater,'element');
md.gridonwater=project3d(md,md.gridonwater,'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

%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;