Index: /issm/trunk/src/m/classes/public/extrude.m
===================================================================
--- /issm/trunk/src/m/classes/public/extrude.m	(revision 2268)
+++ /issm/trunk/src/m/classes/public/extrude.m	(revision 2269)
@@ -1,21 +1,26 @@
-function md=extrude(md,numlayers,extrusionexponent)
+function md=extrude(md,varargin)
 %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');
+%   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~=3) | (nargout~=1)),
+if ((nargin>4) | (nargin<2) | (nargout~=1)),
 	help extrude;
 	error('extrude error message');
@@ -31,4 +36,35 @@
 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...');
@@ -36,23 +72,17 @@
 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
+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-1,
+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*(i/(numlayers-1))^extrusionexponent]; 
+	z3d=[z3d; bed3d+thickness3d*extrusionlist(i)]; 
 end
 number_grids3d=size(x3d,1); %number of 3d grids for the non extruded part of the mesh