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

%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,
	error('number of layers should be at least 2');
end
if md.mesh.dimension==3,
	error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
end

%Initialize with the 2d mesh
x3d=[]; 
y3d=[];
z3d=[];  %the lower node is on the bed
thickness3d=md.geometry.thickness; %thickness and bed for these nodes
bed3d=md.geometry.bed;

%Create the new layers
for i=1:numlayers,
	x3d=[x3d; md.mesh.x]; 
	y3d=[y3d; md.mesh.y];
	%nodes are distributed between bed and surface accordingly to the given exponent
	z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; 
end
number_nodes3d=size(x3d,1); %number of 3d nodes for the non extruded part of the mesh

%Extrude elements 
elements3d=[];
for i=1:numlayers-1,
	elements3d=[elements3d;[md.mesh.elements+(i-1)*md.mesh.numberofvertices md.mesh.elements+i*md.mesh.numberofvertices]]; %Create the elements of the 3d mesh for the non extruded part
end
number_el3d=size(elements3d,1); %number of 3d nodes for the non extruded part of the mesh

%Keep a trace of lower and upper nodes
mesh.lowervertex=NaN*ones(number_nodes3d,1);
mesh.uppervertex=NaN*ones(number_nodes3d,1);
mesh.lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
mesh.uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d;
md.mesh.lowervertex=mesh.lowervertex;
md.mesh.uppervertex=mesh.uppervertex;

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

%Save old mesh 
md.mesh.x2d=md.mesh.x;
md.mesh.y2d=md.mesh.y;
md.mesh.elements2d=md.mesh.elements;
md.mesh.numberofelements2d=md.mesh.numberofelements;
md.mesh.numberofvertices2d=md.mesh.numberofvertices;

%Update mesh type
md.mesh.dimension=3;

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

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

%lat long
md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');

%drag coefficient is limited to nodes that are on the bedrock.
md.friction.coefficient=project3d(md,'vector',md.friction.coefficient,'type','node','layer',1);

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

%observations
md.inversion.vx_obs=project3d(md,'vector',md.inversion.vx_obs,'type','node');
md.inversion.vy_obs=project3d(md,'vector',md.inversion.vy_obs,'type','node');
md.inversion.vel_obs=project3d(md,'vector',md.inversion.vel_obs,'type','node');
md.surfaceforcings.mass_balance=project3d(md,'vector',md.surfaceforcings.mass_balance,'type','node');
md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node');
md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node');

%results
if ~isnan(md.initialization.vx),md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node');end;
if ~isnan(md.initialization.vy),md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node');end;
if ~isnan(md.initialization.vz),md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node');end;
if ~isnan(md.initialization.vel),md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node');end;
if ~isnan(md.initialization.temperature),md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node');end;
if ~isnan(md.initialization.waterfraction),md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node');end;

%bedinfo and surface info
md.mesh.elementonbed=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',1);
md.mesh.elementonsurface=project3d(md,'vector',ones(md.mesh.numberofelements2d,1),'type','element','layer',md.mesh.numberoflayers-1);
md.mesh.vertexonbed=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);

%elementstype
if ~isnan(md.flowequation.element_equation)
	oldelements_type=md.flowequation.element_equation;
	md.flowequation.element_equation=zeros(number_el3d,1);
	md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element');
end

%verticestype
if ~isnan(md.flowequation.vertex_equation)
	oldvertices_type=md.flowequation.vertex_equation;
	md.flowequation.vertex_equation=zeros(number_nodes3d,1);
	md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node');
end
md.flowequation.bordermacayeal=project3d(md,'vector',md.flowequation.bordermacayeal,'type','node');
md.flowequation.borderpattyn=project3d(md,'vector',md.flowequation.borderpattyn,'type','node');
md.flowequation.borderstokes=project3d(md,'vector',md.flowequation.borderstokes,'type','node');

%boundary conditions
md.diagnostic.spcvx=project3d(md,'vector',md.diagnostic.spcvx,'type','node');
md.diagnostic.spcvy=project3d(md,'vector',md.diagnostic.spcvy,'type','node');
md.diagnostic.spcvz=project3d(md,'vector',md.diagnostic.spcvz,'type','node');
md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
md.prognostic.spcthickness=project3d(md,'vector',md.prognostic.spcthickness,'type','node');
md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node');

%in 3d, pressureload: [node1 node2 node3 node4 element]
pressureload_layer1=[md.diagnostic.icefront(:,1:2)  md.diagnostic.icefront(:,2)+md.mesh.numberofvertices2d  md.diagnostic.icefront(:,1)+md.mesh.numberofvertices2d  md.diagnostic.icefront(:,3:4)]; %Add two columns on the first layer 
pressureload=[];
for i=1:numlayers-1,
	pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)];
end
md.diagnostic.icefront=pressureload;

%connectivity
md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
for i=2:numlayers-1,
	md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
		=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
end
md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;

%materials
md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node');

%parameters
md.geometry.surface=project3d(md,'vector',md.geometry.surface,'type','node');
md.geometry.thickness=project3d(md,'vector',md.geometry.thickness,'type','node');
md.geometry.hydrostatic_ratio=project3d(md,'vector',md.geometry.hydrostatic_ratio,'type','node');
md.geometry.bed=project3d(md,'vector',md.geometry.bed,'type','node');
md.geometry.bathymetry=project3d(md,'vector',md.geometry.bathymetry,'type','node');
md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
md.mask.elementonfloatingice=project3d(md,'vector',md.mask.elementonfloatingice,'type','element');
md.mask.vertexonfloatingice=project3d(md,'vector',md.mask.vertexonfloatingice,'type','node');
md.mask.elementongroundedice=project3d(md,'vector',md.mask.elementongroundedice,'type','element');
md.mask.vertexongroundedice=project3d(md,'vector',md.mask.vertexongroundedice,'type','node');
md.mask.elementonwater=project3d(md,'vector',md.mask.elementonwater,'type','element');
md.mask.vertexonwater=project3d(md,'vector',md.mask.vertexonwater,'type','node');
if ~isnan(md.inversion.cost_functions_coefficients),md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end;
if ~isnan(md.inversion.min_parameters),md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node');end;
if ~isnan(md.inversion.max_parameters),md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node');end;
if ~isnan(md.qmu.partition),md.qmu.partition=project3d(md,'vector',md.qmu.partition','type','node');end
if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node');end
if(md.surfaceforcings.isdelta18o),md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node');end
if(md.surfaceforcings.isdelta18o),md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node');end

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

%special for thermal modeling:
md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1); 
if ~isnan(md.basalforcings.geothermalflux)
	md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1); %bedrock only gets geothermal flux
end

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