Index: /issm/trunk-jpl/src/m/classes/mesh2d.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh2d.m	(revision 16287)
+++ /issm/trunk-jpl/src/m/classes/mesh2d.m	(revision 16287)
@@ -0,0 +1,106 @@
+%MESH class definition
+%
+%   Usage:
+%      mesh=mesh();
+
+classdef mesh
+	properties (SetAccess=public) 
+		x                           = NaN;
+		y                           = NaN;
+		elements                    = NaN
+		numberofelements            = 0;
+		numberofvertices            = 0;
+		numberofedges               = 0;
+
+		lat                         = NaN
+		long                        = NaN
+		hemisphere                  = NaN
+
+		vertexonboundary            = NaN
+
+		edges                       = NaN
+		segments                    = NaN
+		segmentmarkers              = NaN
+		vertexconnectivity          = NaN
+		elementconnectivity         = NaN
+		average_vertex_connectivity = 0;
+
+		extractedvertices           = NaN
+		extractedelements           = NaN
+	end
+	methods
+		function obj = mesh(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%the connectivity is the averaged number of nodes linked to a
+			%given node through an edge. This connectivity is used to initially
+			%allocate memory to the stiffness matrix. A value of 16 seems to
+			%give a good memory/time ration. This value can be checked in
+			%trunk/test/Miscellaneous/runme.m
+			obj.average_vertex_connectivity=25;
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+
+			md = checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'mesh.elements','NaN',1,'>',0,'values',1:md.mesh.numberofvertices);
+			md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements 3]);
+			if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.mesh.elements(:)))));
+				md = checkmessage(md,'orphan nodes have been found. Check the mesh outline');
+			end
+			md = checkfield(md,'mesh.numberofelements','>',0);
+			md = checkfield(md,'mesh.numberofvertices','>',0);
+			md = checkfield(md,'mesh.vertexonbed','size',[md.mesh.numberofvertices 1],'values',[0 1]);
+			md = checkfield(md,'mesh.average_vertex_connectivity','>=',9,'message','''mesh.average_vertex_connectivity'' should be at least 9 in 2d');
+
+			switch(solution),
+				case ThermalSolutionEnum(),
+					md = checkmessage(md,'thermal not supported for 2d mesh');
+			end
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   2d Mesh:')); 
+
+			disp(sprintf('\n      Elements and vertices:'));
+			fielddisplay(obj,'numberofelements','number of elements');
+			fielddisplay(obj,'numberofvertices','number of vertices');
+			fielddisplay(obj,'elements','vertex indices of the mesh elements');
+			fielddisplay(obj,'x','vertices x coordinate [m]');
+			fielddisplay(obj,'y','vertices y coordinate [m]');
+			fielddisplay(obj,'edges','edges of the 2d mesh (vertex1 vertex2 element1 element2)');
+			fielddisplay(obj,'numberofedges','number of edges of the 2d mesh');
+
+			disp(sprintf('\n      Properties:'));
+			fielddisplay(obj,'vertexonboundary','vertices on the boundary of the domain flag list');
+			fielddisplay(obj,'segments','edges on domain boundary (vertex1 vertex2 element)');
+			fielddisplay(obj,'segmentmarkers','number associated to each segment');
+			fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i');
+			fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i');
+			fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
+
+			disp(sprintf('\n      Extracted model:'));
+			fielddisplay(obj,'extractedvertices','vertices extracted from the model');
+			fielddisplay(obj,'extractedelements','elements extracted from the model');
+
+			disp(sprintf('\n      Projection:'));
+			fielddisplay(obj,'lat','vertices latitude [degrees]');
+			fielddisplay(obj,'long','vertices longitude [degrees]');
+			fielddisplay(obj,'hemisphere','Indicate hemisphere ''n'' or ''s'' ');
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'object',obj,'fieldname','x','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','y','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
+		end % }}}
+	end
+end
Index: /issm/trunk-jpl/src/m/classes/mesh2d.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh2d.py	(revision 16287)
+++ /issm/trunk-jpl/src/m/classes/mesh2d.py	(revision 16287)
@@ -0,0 +1,107 @@
+import numpy
+from fielddisplay import fielddisplay
+from EnumDefinitions import *
+from checkfield import *
+from MatlabFuncs import *
+
+class mesh(object):
+	"""
+	MESH class definition
+
+	   Usage:
+	      mesh=mesh();
+	"""
+
+	def __init__(self): # {{{
+		self.x                           = float('NaN');
+		self.y                           = float('NaN');
+		self.elements                    = float('NaN');
+		self.numberofelements            = 0;
+		self.numberofvertices            = 0;
+		self.numberofedges               = 0;
+		
+		self.lat                         = float('NaN');
+		self.long                        = float('NaN');
+		self.hemisphere                  = float('NaN');
+
+		self.vertexonboundary            = float('NaN');
+		self.edges                       = float('NaN');
+		self.segments                    = float('NaN');
+		self.segmentmarkers              = float('NaN');
+		self.vertexconnectivity          = float('NaN');
+		self.elementconnectivity         = float('NaN');
+		self.average_vertex_connectivity = 0;
+
+		self.extractedvertices           = float('NaN');
+		self.extractedelements           = float('NaN');
+
+		#set defaults
+		self.setdefaultparameters()
+
+		#}}}
+	def __repr__(self): # {{{
+		string="   2d Mesh:" 
+
+		string="%s\n%s"%(string,"\n      Elements and vertices:")
+		string="%s\n%s"%(string,fielddisplay(self,"numberofelements","number of elements"))
+		string="%s\n%s"%(string,fielddisplay(self,"numberofvertices","number of vertices"))
+		string="%s\n%s"%(string,fielddisplay(self,"elements","vertex indices of the mesh elements"))
+		string="%s\n%s"%(string,fielddisplay(self,"x","vertices x coordinate [m]"))
+		string="%s\n%s"%(string,fielddisplay(self,"y","vertices y coordinate [m]"))
+		string="%s\n%s"%(string,fielddisplay(self,"edges","edges of the 2d mesh (vertex1 vertex2 element1 element2)"))
+		string="%s\n%s"%(string,fielddisplay(self,"numberofedges","number of edges of the 2d mesh"))
+
+		string="%s%s"%(string,"\n\n      Properties:")
+		string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
+		string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)"))
+		string="%s\n%s"%(string,fielddisplay(self,"segmentmarkers","number associated to each segment"))
+		string="%s\n%s"%(string,fielddisplay(self,"vertexconnectivity","list of vertices connected to vertex_i"))
+		string="%s\n%s"%(string,fielddisplay(self,"elementconnectivity","list of vertices connected to element_i"))
+		string="%s\n%s"%(string,fielddisplay(self,"average_vertex_connectivity","average number of vertices connected to one vertex"))
+
+		string="%s%s"%(string,"\n\n      Extracted model:")
+		string="%s\n%s"%(string,fielddisplay(self,"extractedvertices","vertices extracted from the model"))
+		string="%s\n%s"%(string,fielddisplay(self,"extractedelements","elements extracted from the model"))
+
+		string="%s%s"%(string,"\n\n      Projection:")
+		string="%s\n%s"%(string,fielddisplay(self,"lat","vertices latitude [degrees]"))
+		string="%s\n%s"%(string,fielddisplay(self,"long","vertices longitude [degrees]"))
+		string="%s\n%s"%(string,fielddisplay(self,"hemisphere","Indicate hemisphere 'n' or 's'"))
+		return string
+		#}}}
+	def setdefaultparameters(self): # {{{
+		
+		#the connectivity is the averaged number of nodes linked to a
+		#given node through an edge. This connectivity is used to initially
+		#allocate memory to the stiffness matrix. A value of 16 seems to
+		#give a good memory/time ration. This value can be checked in
+		#trunk/test/Miscellaneous/runme.m
+		self.average_vertex_connectivity=25
+
+		return self
+	#}}}
+	def checkconsistency(self,md,solution,analyses):    # {{{
+
+		md = checkfield(md,'mesh.x','NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'mesh.y','NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'mesh.z','NaN',1,'size',[md.mesh.numberofvertices])
+		md = checkfield(md,'mesh.elements','NaN',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1))
+		md = checkfield(md,'mesh.elements','size',[md.mesh.numberofelements,3])
+		if numpy.any(numpy.logical_not(ismember(numpy.arange(1,md.mesh.numberofvertices+1),md.mesh.elements))):
+			md.checkmessage("orphan nodes have been found. Check the mesh outline")
+		md = checkfield(md,'mesh.numberofelements','>',0)
+		md = checkfield(md,'mesh.numberofvertices','>',0)
+		md = checkfield(md,'mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d")
+		if solution==ThermalSolutionEnum():
+			md.checkmessage("thermal not supported for 2d mesh")
+
+		return md
+	# }}}
+	def marshall(self,md,fid):    # {{{
+		WriteData(fid,'object',self,'fieldname','x','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','y','format','DoubleMat','mattype',1)
+		WriteData(fid,'object',self,'fieldname','elements','format','DoubleMat','mattype',2)
+		WriteData(fid,'object',self,'fieldname','numberofelements','format','Integer')
+		WriteData(fid,'object',self,'fieldname','numberofvertices','format','Integer')
+		WriteData(fid,'object',self,'fieldname','average_vertex_connectivity','format','Integer')
+	# }}}
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 16287)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 16287)
@@ -0,0 +1,1204 @@
+%MODEL class definition
+%
+%   Usage:
+%      md = model(varargin)
+
+classdef model
+	properties (SetAccess=public) %Model fields
+		% {{{
+		%Careful here: no other class should be used as default value this is a bug of matlab
+		mesh             = 0;
+		mask             = 0;
+
+		geometry         = 0;
+		constants        = 0;
+		surfaceforcings  = 0;
+		basalforcings    = 0;
+		materials        = 0;
+		damage           = 0;
+		friction         = 0;
+		flowequation     = 0;
+		timestepping     = 0;
+		initialization   = 0;
+		rifts            = 0;
+
+		debug            = 0;
+		verbose          = 0;
+		settings         = 0;
+		toolkits         = 0;
+		cluster          = 0;
+
+		balancethickness = 0;
+		stressbalance       = 0;
+		groundingline    = 0;
+		hydrology        = 0;
+		masstransport       = 0;
+		thermal          = 0;
+		steadystate      = 0;
+		transient        = 0;
+		gia              = 0;
+
+		autodiff         = 0;
+		flaim            = 0;
+		inversion        = 0;
+		qmu              = 0;
+
+		results          = 0;
+		radaroverlay     = 0;
+		miscellaneous    = 0;
+		private          = 0;
+
+		%}}}
+	end
+	methods (Static)
+		function md = loadobj(md) % {{{
+			% This function is directly called by matlab when a model object is
+			% loaded. If the input is a struct it is an old version of model and
+			% old fields must be recovered (make sure they are in the deprecated
+			% model properties)
+
+			if verLessThan('matlab','7.9'),
+				disp('Warning: your matlab version is old and there is a risk that load does not work correctly');
+				disp('         if the model is not loaded correctly, rename temporarily loadobj so that matlab does not use it');
+
+				% This is a Matlab bug: all the fields of md have their default value
+				% Example of error message:
+				% Warning: Error loading an object of class 'model':
+				% Undefined function or method 'exist' for input arguments of type 'cell'
+				%
+				% This has been fixed in MATLAB 7.9 (R2009b) and later versions
+			end
+
+			if isstruct(md)
+				disp('Recovering model object from a previous version');
+				md = structtomodel(model,md);
+			end
+
+			%2012 August 4th
+			if isa(md.materials,'materials'),
+				disp('Recovering old materials');
+				if numel(md.materials.rheology_Z)==1 & isnan(md.materials.rheology_Z),
+					md.materials=matice(md.materials);
+				else
+					md.materials=matdamageice(md.materials);
+				end
+			end
+			%2013 April 12
+			if numel(md.stressbalance.loadingforce==1)
+				md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+			end
+			%2013 April 17
+			if isa(md.hydrology,'hydrology'),
+				disp('Recovering old hydrology class');
+				md.hydrology=hydrologyshreve(md.materials);
+			end
+		end% }}}
+	end
+	methods
+		function md = model(varargin) % {{{
+
+			switch nargin
+				case 0
+					md=setdefaultparameters(md);
+				otherwise
+					error('model constructor error message: 0 of 1 argument only in input.');
+				end
+		end
+		%}}}
+		function md = checkmessage(md,string) % {{{
+			if(nargout~=1) error('wrong usage, model must be an output'); end
+			disp(['model not consistent: ' string]);
+			md.private.isconsistent=false;
+		end
+		%}}}
+		function md = collapse(md)% {{{
+			%COLLAPSE - collapses a 3d mesh into a 2d mesh
+			%
+			%   This routine collapses a 3d model into a 2d model
+			%   and collapses all the fileds of the 3d model by
+			%   taking their depth-averaged values
+			%
+			%   Usage:
+			%      md=collapse(md)
+			%
+			%   See also: EXTRUDE, MODELEXTRACT
+
+			%Check that the model is really a 3d model
+			if ~md.mesh.dimension==3,
+				error('collapse error message: only 3d mesh can be collapsed')
+			end
+
+			%Start with changing alle the fields from the 3d mesh 
+
+			%drag is limited to nodes that are on the bedrock.
+			md.friction.coefficient=project2d(md,md.friction.coefficient,1);
+
+			%p and q (same deal, except for element that are on the bedrock: )
+			md.friction.p=project2d(md,md.friction.p,1);
+			md.friction.q=project2d(md,md.friction.q,1);
+
+			%observations
+			if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
+			if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
+			if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
+			if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end;
+			if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end;
+			if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end;
+			if ~isnan(md.surfaceforcings.mass_balance),
+				md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); 
+			end;
+			if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end;
+
+			%results
+			if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
+			if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
+			if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
+			if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;
+			if ~isnan(md.initialization.temperature),md.initialization.temperature=DepthAverage(md,md.initialization.temperature);end;
+
+			%gia
+			if ~isnan(md.gia.mantle_viscosity), md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); end
+			if ~isnan(md.gia.lithosphere_thickness), md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); end
+
+			%bedinfo and surface info
+			md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1);
+			md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1);
+			md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1);
+			md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1);
+
+			%elementstype
+			if ~isnan(md.flowequation.element_equation)
+				md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
+				md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
+				md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
+				md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
+				md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
+			end	
+
+			%boundary conditions
+			md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
+			md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
+			md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
+			md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
+			md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
+			md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
+			md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
+
+			%materials
+			md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
+			md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
+			
+			%damage: 
+			md.damage.D=DepthAverage(md,md.damage.D);
+
+			%special for thermal modeling:
+			md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); 
+			md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux
+
+			%update of connectivity matrix
+			md.mesh.average_vertex_connectivity=25;
+
+			%Collapse the mesh
+			nodes2d=md.mesh.numberofvertices2d;
+			elements2d=md.mesh.numberofelements2d;
+
+			%parameters
+			md.geometry.surface=project2d(md,md.geometry.surface,1);
+			md.geometry.thickness=project2d(md,md.geometry.thickness,1);
+			md.geometry.bed=project2d(md,md.geometry.bed,1);
+			md.geometry.bathymetry=project2d(md,md.geometry.bathymetry,1);
+			md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
+			md.mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
+			md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
+			md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
+
+			%lat long
+			if numel(md.mesh.lat) ==md.mesh.numberofvertices,  md.mesh.lat=project2d(md,md.mesh.lat,1); end
+			if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
+
+			%Initialize with the 2d mesh
+			md.mesh.x=md.mesh.x2d;
+			md.mesh.y=md.mesh.y2d;
+			md.mesh.z=zeros(size(md.mesh.x2d));
+			md.mesh.numberofvertices=md.mesh.numberofvertices2d;
+			md.mesh.numberofelements=md.mesh.numberofelements2d;
+			md.mesh.elements=md.mesh.elements2d;
+
+			%Keep a trace of lower and upper nodes
+			md.mesh.lowervertex=NaN;
+			md.mesh.uppervertex=NaN;
+			md.mesh.lowerelements=NaN;
+			md.mesh.upperelements=NaN;
+
+			%Remove old mesh 
+			md.mesh.x2d=NaN;
+			md.mesh.y2d=NaN;
+			md.mesh.elements2d=NaN;
+			md.mesh.numberofelements2d=md.mesh.numberofelements;
+			md.mesh.numberofvertices2d=md.mesh.numberofvertices;
+			md.mesh.numberoflayers=0;
+
+			%Update mesh type
+			md.mesh.dimension=2;
+		end % }}}
+		function md2 = extract(md,area) % {{{
+			%extract - extract a model according to an Argus contour or flag list
+			%
+			%   This routine extracts a submodel from a bigger model with respect to a given contour
+			%   md must be followed by the corresponding exp file or flags list
+			%   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
+			%   If user wants every element outside the domain to be 
+			%   extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
+			%   an empty string '' will be considered as an empty domain
+			%   a string 'all' will be considered as the entire domain
+			%
+			%   Usage:
+			%      md2=extract(md,area);
+			%
+			%   Examples:
+			%      md2=extract(md,'Domain.exp');
+			%
+			%   See also: EXTRUDE, COLLAPSE
+
+			%copy model
+			md1=md;
+
+			%some checks
+			if ((nargin~=2) | (nargout~=1)),
+				help extract
+				error('extract error message: bad usage');
+			end
+
+			%get elements that are inside area
+			flag_elem=FlagElements(md1,area);
+			if ~any(flag_elem),
+				error('extracted model is empty');
+			end
+
+			%kick out all elements with 3 dirichlets
+			spc_elem=find(~flag_elem);
+			spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
+			flag=ones(md1.mesh.numberofvertices,1);
+			flag(spc_node)=0;
+			pos=find(sum(flag(md1.mesh.elements),2)==0);
+			flag_elem(pos)=0;
+
+			%extracted elements and nodes lists
+			pos_elem=find(flag_elem);
+			pos_node=sort(unique(md1.mesh.elements(pos_elem,:)));
+
+			%keep track of some fields
+			numberofvertices1=md1.mesh.numberofvertices;
+			numberofelements1=md1.mesh.numberofelements;
+			numberofvertices2=length(pos_node);
+			numberofelements2=length(pos_elem);
+			flag_node=zeros(numberofvertices1,1);
+			flag_node(pos_node)=1;
+
+			%Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
+			Pelem=zeros(numberofelements1,1);
+			Pelem(pos_elem)=[1:numberofelements2]';
+			Pnode=zeros(numberofvertices1,1);
+			Pnode(pos_node)=[1:numberofvertices2]';
+
+			%renumber the elements (some nodes won't exist anymore)
+			elements_1=md1.mesh.elements;
+			elements_2=elements_1(pos_elem,:);
+			elements_2(:,1)=Pnode(elements_2(:,1));
+			elements_2(:,2)=Pnode(elements_2(:,2));
+			elements_2(:,3)=Pnode(elements_2(:,3));
+			if md1.mesh.dimension==3,
+				elements_2(:,4)=Pnode(elements_2(:,4));
+				elements_2(:,5)=Pnode(elements_2(:,5));
+				elements_2(:,6)=Pnode(elements_2(:,6));
+			end
+
+			%OK, now create the new model!
+
+			%take every field from model
+			md2=md1;
+
+			%automatically modify fields
+
+			%loop over model fields
+			model_fields=fields(md1);
+			for i=1:length(model_fields),
+				%get field
+				field=md1.(model_fields{i});
+				fieldsize=size(field);
+				if isobject(field), %recursive call
+					object_fields=fields(md1.(model_fields{i}));
+					for j=1:length(object_fields),
+						%get field
+						field=md1.(model_fields{i}).(object_fields{j});
+						fieldsize=size(field);
+						%size = number of nodes * n
+						if fieldsize(1)==numberofvertices1
+							md2.(model_fields{i}).(object_fields{j})=field(pos_node,:);
+						elseif (fieldsize(1)==numberofvertices1+1)
+							md2.(model_fields{i}).(object_fields{j})=[field(pos_node,:); field(end,:)];
+						%size = number of elements * n
+						elseif fieldsize(1)==numberofelements1
+							md2.(model_fields{i}).(object_fields{j})=field(pos_elem,:);
+						end
+					end
+				else
+					%size = number of nodes * n
+					if fieldsize(1)==numberofvertices1
+						md2.(model_fields{i})=field(pos_node,:);
+					elseif (fieldsize(1)==numberofvertices1+1)
+						md2.(model_fields{i})=[field(pos_node,:); field(end,:)];
+					%size = number of elements * n
+					elseif fieldsize(1)==numberofelements1
+						md2.(model_fields{i})=field(pos_elem,:);
+					end
+				end
+			end
+
+			%modify some specific fields
+
+			%Mesh
+			md2.mesh.numberofelements=numberofelements2;
+			md2.mesh.numberofvertices=numberofvertices2;
+			md2.mesh.elements=elements_2;
+
+			%mesh.uppervertex mesh.lowervertex
+			if md1.mesh.dimension==3
+				md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
+				pos=find(~isnan(md2.mesh.uppervertex));
+				md2.mesh.uppervertex(pos)=Pnode(md2.mesh.uppervertex(pos));
+
+				md2.mesh.lowervertex=md1.mesh.lowervertex(pos_node);
+				pos=find(~isnan(md2.mesh.lowervertex));
+				md2.mesh.lowervertex(pos)=Pnode(md2.mesh.lowervertex(pos));
+
+				md2.mesh.upperelements=md1.mesh.upperelements(pos_elem);
+				pos=find(~isnan(md2.mesh.upperelements));
+				md2.mesh.upperelements(pos)=Pelem(md2.mesh.upperelements(pos));
+
+				md2.mesh.lowerelements=md1.mesh.lowerelements(pos_elem);
+				pos=find(~isnan(md2.mesh.lowerelements));
+				md2.mesh.lowerelements(pos)=Pelem(md2.mesh.lowerelements(pos));
+			end
+
+			%Initial 2d mesh 
+			if md1.mesh.dimension==3
+				flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
+				pos_elem_2d=find(flag_elem_2d);
+				flag_node_2d=flag_node(1:md1.mesh.numberofvertices2d);
+				pos_node_2d=find(flag_node_2d);
+
+				md2.mesh.numberofelements2d=length(pos_elem_2d);
+				md2.mesh.numberofvertices2d=length(pos_node_2d);
+				md2.mesh.elements2d=md1.mesh.elements2d(pos_elem_2d,:);
+				md2.mesh.elements2d(:,1)=Pnode(md2.mesh.elements2d(:,1));
+				md2.mesh.elements2d(:,2)=Pnode(md2.mesh.elements2d(:,2));
+				md2.mesh.elements2d(:,3)=Pnode(md2.mesh.elements2d(:,3));
+
+				md2.mesh.x2d=md1.mesh.x(pos_node_2d);
+				md2.mesh.y2d=md1.mesh.y(pos_node_2d);
+			end
+
+			%Edges
+			if size(md2.mesh.edges,2)>1, %do not use ~isnan because there are some NaNs...
+				%renumber first two columns
+				pos=find(md2.mesh.edges(:,4)~=-1);
+				md2.mesh.edges(:  ,1)=Pnode(md2.mesh.edges(:,1));
+				md2.mesh.edges(:  ,2)=Pnode(md2.mesh.edges(:,2));
+				md2.mesh.edges(:  ,3)=Pelem(md2.mesh.edges(:,3));
+				md2.mesh.edges(pos,4)=Pelem(md2.mesh.edges(pos,4));
+				%remove edges when the 2 vertices are not in the domain.
+				md2.mesh.edges=md2.mesh.edges(find(md2.mesh.edges(:,1) & md2.mesh.edges(:,2)),:);
+				%Replace all zeros by -1 in the last two columns
+				pos=find(md2.mesh.edges(:,3)==0);
+				md2.mesh.edges(pos,3)=-1;
+				pos=find(md2.mesh.edges(:,4)==0);
+				md2.mesh.edges(pos,4)=-1;
+				%Invert -1 on the third column with last column (Also invert first two columns!!)
+				pos=find(md2.mesh.edges(:,3)==-1);
+				md2.mesh.edges(pos,3)=md2.mesh.edges(pos,4);
+				md2.mesh.edges(pos,4)=-1;
+				values=md2.mesh.edges(pos,2);
+				md2.mesh.edges(pos,2)=md2.mesh.edges(pos,1);
+				md2.mesh.edges(pos,1)=values;
+				%Finally remove edges that do not belong to any element
+				pos=find(md2.mesh.edges(:,3)==-1 & md2.mesh.edges(:,4)==-1);
+				md2.mesh.edges(pos,:)=[];
+			end
+
+			%Penalties
+			if ~isnan(md2.stressbalance.vertex_pairing),
+				for i=1:size(md1.stressbalance.vertex_pairing,1);
+					md2.stressbalance.vertex_pairing(i,:)=Pnode(md1.stressbalance.vertex_pairing(i,:));
+				end
+				md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing(find(md2.stressbalance.vertex_pairing(:,1)),:);
+			end
+			if ~isnan(md2.masstransport.vertex_pairing),
+				for i=1:size(md1.masstransport.vertex_pairing,1);
+					md2.masstransport.vertex_pairing(i,:)=Pnode(md1.masstransport.vertex_pairing(i,:));
+				end
+				md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing(find(md2.masstransport.vertex_pairing(:,1)),:);
+			end
+
+			%recreate segments
+			if md1.mesh.dimension==2
+				md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
+				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
+				md2.mesh.segments=contourenvelope(md2);
+				md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+			else
+				%First do the connectivity for the contourenvelope in 2d
+				md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
+				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
+				md2.mesh.segments=contourenvelope(md2);
+				md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+				md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
+				%Then do it for 3d as usual
+				md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
+				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
+			end
+
+			%Boundary conditions: Dirichlets on new boundary
+			%Catch the elements that have not been extracted
+			orphans_elem=find(~flag_elem);
+			orphans_node=unique(md1.mesh.elements(orphans_elem,:))';
+			%Figure out which node are on the boundary between md2 and md1
+			nodestoflag1=intersect(orphans_node,pos_node);
+			nodestoflag2=Pnode(nodestoflag1);
+			if numel(md1.stressbalance.spcvx)>1 & numel(md1.stressbalance.spcvy)>2 & numel(md1.stressbalance.spcvz)>2,
+				if numel(md1.inversion.vx_obs)>1 & numel(md1.inversion.vy_obs)>1
+					md2.stressbalance.spcvx(nodestoflag2)=md2.inversion.vx_obs(nodestoflag2); 
+					md2.stressbalance.spcvy(nodestoflag2)=md2.inversion.vy_obs(nodestoflag2);
+				else
+					md2.stressbalance.spcvx(nodestoflag2)=NaN;
+					md2.stressbalance.spcvy(nodestoflag2)=NaN;
+					disp(' ')
+					disp('!! extract warning: spc values should be checked !!')
+					disp(' ')
+				end
+				%put 0 for vz
+				md2.stressbalance.spcvz(nodestoflag2)=0;
+			end
+			if ~isnan(md1.thermal.spctemperature),
+				md2.thermal.spctemperature(nodestoflag2,1)=1;
+			end
+
+			%Results fields
+			if isstruct(md1.results),
+				md2.results=struct();
+				solutionfields=fields(md1.results);
+				for i=1:length(solutionfields),
+					if isstruct(md1.results.(solutionfields{i}))
+						%get subfields
+						solutionsubfields=fields(md1.results.(solutionfields{i}));
+						for j=1:length(solutionsubfields),
+							field=md1.results.(solutionfields{i}).(solutionsubfields{j});
+							if length(field)==numberofvertices1,
+								md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_node);
+							elseif length(field)==numberofelements1,
+								md2.results.(solutionfields{i}).(solutionsubfields{j})=field(pos_elem);
+							else
+								md2.results.(solutionfields{i}).(solutionsubfields{j})=field;
+							end
+						end
+					else
+						field=md1.results.(solutionfields{i});
+						if length(field)==numberofvertices1,
+							md2.results.(solutionfields{i})=field(pos_node);
+						elseif length(field)==numberofelements1,
+							md2.results.(solutionfields{i})=field(pos_elem);
+						else
+							md2.results.(solutionfields{i})=field;
+						end
+					end
+				end
+			end
+
+			%Keep track of pos_node and pos_elem
+			md2.mesh.extractedvertices=pos_node;
+			md2.mesh.extractedelements=pos_elem;
+		end % }}}
+		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
+				clist=varargin{1};
+				if any(clist<0) | any(clist>1),
+					error('extrusioncoefficients must be between 0 and 1');
+				end
+				extrusionlist=sort(unique([clist(:);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;
+            if ~isnan(md.initialization.watercolumn),md.initialization.watercolumn=project3d(md,'vector',md.initialization.watercolumn,'type','node','layer',1);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.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node');
+			md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node');
+			md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node');
+
+			%boundary conditions
+			md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node');
+			md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node');
+			md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node');
+			md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
+			md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node');
+			md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
+			md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node');
+			md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
+			md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
+
+			%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');
+				
+			%damage
+			md.damage.D=project3d(md,'vector',md.damage.D,'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.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node');
+			md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_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.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node');
+			md.mask.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'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
+			end % }}}
+		function md = structtomodel(md,structmd) % {{{
+
+			if ~isstruct(structmd) error('input model is not a structure'); end
+
+			%loaded model is a struct, initialize output and recover all fields
+			md = structtoobj(model,structmd);
+
+			%Old field now classes
+			if (isfield(structmd,'timestepping') & isnumeric(md.timestepping)), md.timestepping=timestepping(); end
+			if (isfield(structmd,'mask') & isnumeric(md.mask)),md.mask=mask(); end
+
+			%Field name change
+			if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
+			if isfield(structmd,'p'), md.friction.p=structmd.p; end
+			if isfield(structmd,'q'), md.friction.q=structmd.p; end
+			if isfield(structmd,'melting'), md.basalforcings.melting_rate=structmd.melting; end
+			if isfield(structmd,'melting_rate'), md.basalforcings.melting_rate=structmd.melting_rate; end
+			if isfield(structmd,'accumulation'), md.surfaceforcings.mass_balance=structmd.accumulation; end
+			if isfield(structmd,'numberofgrids'), md.mesh.numberofvertices=structmd.numberofgrids; end
+			if isfield(structmd,'numberofgrids2d'), md.mesh.numberofvertices2d=structmd.numberofgrids2d; end
+			if isfield(structmd,'uppergrids'), md.mesh.uppervertex=structmd.uppergrids; end
+			if isfield(structmd,'lowergrids'), md.mesh.lowervertex=structmd.lowergrids; end
+			if isfield(structmd,'gridonbed'), md.mesh.vertexonbed=structmd.gridonbed; end
+			if isfield(structmd,'gridonsurface'), md.mesh.vertexonsurface=structmd.gridonsurface; end
+			if isfield(structmd,'extractedgrids'), md.mesh.extractedvertices=structmd.extractedgrids; end
+			if isfield(structmd,'gridonboundary'), md.mesh.vertexonboundary=structmd.gridonboundary; end
+			if isfield(structmd,'petscoptions') & ~isempty(structmd.petscoptions), md.toolkits=structmd.petscoptions; end
+			if isfield(structmd,'g'), md.constants.g=structmd.g; end
+			if isfield(structmd,'yts'), md.constants.yts=structmd.yts; end
+			if isfield(structmd,'surface_mass_balance'), md.surfaceforcings.mass_balance=structmd.surface_mass_balance; end
+			if isfield(structmd,'basal_melting_rate'), md.basalforcings.melting_rate=structmd.basal_melting_rate; end
+			if isfield(structmd,'basal_melting_rate_correction'), md.basalforcings.melting_rate_correction=structmd.basal_melting_rate_correction; end
+			if isfield(structmd,'geothermalflux'), md.basalforcings.geothermalflux=structmd.geothermalflux; end
+			if isfield(structmd,'drag'), md.friction.coefficient=structmd.drag; end
+			if isfield(structmd,'drag_coefficient'), md.friction.coefficient=structmd.drag_coefficient; end
+			if isfield(structmd,'drag_p'), md.friction.p=structmd.drag_p; end
+			if isfield(structmd,'drag_q'), md.friction.q=structmd.drag_q; end
+			if isfield(structmd,'riftproperties'), %old implementation
+				md.rifts=rifts();
+				md.rifts.riftproperties=structmd.riftproperties; 
+				md.rifts.riftstruct=structmd.rifts;
+				md.rifts.riftproperties=structmd.riftinfo;
+			end
+			if isfield(structmd,'bamg'), md.private.bamg=structmd.bamg; end
+			if isfield(structmd,'lowmem'), md.settings.lowmem=structmd.lowmem; end
+			if isfield(structmd,'io_gather'), md.settings.io_gather=structmd.io_gather; end
+			if isfield(structmd,'spcwatercolumn'), md.hydrology.spcwatercolumn=structmd.spcwatercolumn; end
+			if isfield(structmd,'hydro_n'), md.hydrology.n=structmd.hydro_n; end
+			if isfield(structmd,'hydro_p'), md.hydrology.p=structmd.hydro_p; end
+			if isfield(structmd,'hydro_q'), md.hydrology.q=structmd.hydro_q; end
+			if isfield(structmd,'hydro_CR'), md.hydrology.CR=structmd.hydro_CR; end
+			if isfield(structmd,'hydro_kn'), md.hydrology.kn=structmd.hydro_kn; end
+			if isfield(structmd,'spctemperature'), md.thermal.spctemperature=structmd.spctemperature; end
+			if isfield(structmd,'min_thermal_constraints'), md.thermal.penalty_threshold=structmd.min_thermal_constraints; end
+			if isfield(structmd,'artificial_diffusivity'), md.thermal.stabilization=structmd.artificial_diffusivity; end
+			if isfield(structmd,'max_nonlinear_iterations'), md.thermal.maxiter=structmd.max_nonlinear_iterations; end
+			if isfield(structmd,'stabilize_constraints'), md.thermal.penalty_lock=structmd.stabilize_constraints; end
+			if isfield(structmd,'penalty_offset'), md.thermal.penalty_factor=structmd.penalty_offset; end
+			if isfield(structmd,'name'), md.miscellaneous.name=structmd.name; end
+			if isfield(structmd,'notes'), md.miscellaneous.notes=structmd.notes; end
+			if isfield(structmd,'dummy'), md.miscellaneous.dummy=structmd.dummy; end
+			if isfield(structmd,'dt'), md.timestepping.time_step=structmd.dt; end
+			if isfield(structmd,'ndt'), md.timestepping.final_time=structmd.ndt; end
+			if isfield(structmd,'time_adapt'), md.timestepping.time_adapt=structmd.time_adapt; end
+			if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end
+			if isfield(structmd,'spcthickness'), md.masstransport.spcthickness=structmd.spcthickness; end
+			if isfield(structmd,'artificial_diffusivity'), md.masstransport.stabilization=structmd.artificial_diffusivity; end
+			if isfield(structmd,'hydrostatic_adjustment'), md.masstransport.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end
+			if isfield(structmd,'penalties'), md.masstransport.vertex_pairing=structmd.penalties; end
+			if isfield(structmd,'penalty_offset'), md.masstransport.penalty_factor=structmd.penalty_offset; end
+			if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
+			if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
+			if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
+			if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
+			if isfield(structmd,'rheology_Z'), md.damage.D=(1-structmd.rheology_Z); end
+			if isfield(structmd,'spcthickness'), md.balancethickness.spcthickness=structmd.spcthickness; end
+			if isfield(structmd,'artificial_diffusivity'), md.balancethickness.stabilization=structmd.artificial_diffusivity; end
+			if isfield(structmd,'dhdt'), md.balancethickness.thickening_rate=structmd.dhdt; end
+			if isfield(structmd,'isSIA'), md.flowequation.isSIA=structmd.isSIA; end
+			if isfield(structmd,'isFS'), md.flowequation.isFS=structmd.isFS; end
+			if isfield(structmd,'elements_type'), md.flowequation.element_equation=structmd.elements_type; end
+			if isfield(structmd,'vertices_type'), md.flowequation.vertex_equation=structmd.vertices_type; end
+			if isfield(structmd,'eps_rel'), md.steadystate.reltol=structmd.eps_rel; end
+			if isfield(structmd,'max_steadystate_iterations'), md.steadystate.maxiter=structmd.max_steadystate_iterations; end
+			if isfield(structmd,'isdiagnostic'), md.transient.isstressbalance=structmd.isdiagnostic; end
+			if isfield(structmd,'isprognostic'), md.transient.ismasstransport=structmd.isprognostic; end
+			if isfield(structmd,'isthermal'), md.transient.isthermal=structmd.isthermal; end
+			if isfield(structmd,'control_analysis'), md.inversion.iscontrol=structmd.control_analysis; end
+			if isfield(structmd,'weights'), md.inversion.cost_functions_coefficients=structmd.weights; end
+			if isfield(structmd,'nsteps'), md.inversion.nsteps=structmd.nsteps; end
+			if isfield(structmd,'maxiter_per_step'), md.inversion.maxiter_per_step=structmd.maxiter_per_step; end
+			if isfield(structmd,'cm_min'), md.inversion.min_parameters=structmd.cm_min; end
+			if isfield(structmd,'cm_max'), md.inversion.max_parameters=structmd.cm_max; end
+			if isfield(structmd,'vx_obs'), md.inversion.vx_obs=structmd.vx_obs; end
+			if isfield(structmd,'vy_obs'), md.inversion.vy_obs=structmd.vy_obs; end
+			if isfield(structmd,'vel_obs'), md.inversion.vel_obs=structmd.vel_obs; end
+			if isfield(structmd,'thickness_obs'), md.inversion.thickness_obs=structmd.thickness_obs; end
+			if isfield(structmd,'vx'), md.initialization.vx=structmd.vx; end
+			if isfield(structmd,'vy'), md.initialization.vy=structmd.vy; end
+			if isfield(structmd,'vz'), md.initialization.vz=structmd.vz; end
+			if isfield(structmd,'vel'), md.initialization.vel=structmd.vel; end
+			if isfield(structmd,'pressure'), md.initialization.pressure=structmd.pressure; end
+			if isfield(structmd,'temperature'), md.initialization.temperature=structmd.temperature; end
+			if isfield(structmd,'waterfraction'), md.initialization.waterfraction=structmd.waterfraction; end
+			if isfield(structmd,'watercolumn'), md.initialization.watercolumn=structmd.watercolumn; end
+			if isfield(structmd,'surface'), md.geometry.surface=structmd.surface; end
+			if isfield(structmd,'bed'), md.geometry.bed=structmd.bed; end
+			if isfield(structmd,'thickness'), md.geometry.thickness=structmd.thickness; end
+			if isfield(structmd,'bathymetry'), md.geometry.bathymetry=structmd.bathymetry; end
+			if isfield(structmd,'thickness_coeff'), md.geometry.hydrostatic_ratio=structmd.thickness_coeff; end
+			if isfield(structmd,'connectivity'), md.mesh.average_vertex_connectivity=structmd.connectivity; end
+			if isfield(structmd,'extractednodes'), md.mesh.extractedvertices=structmd.extractednodes; end
+			if isfield(structmd,'extractedelements'), md.mesh.extractedelements=structmd.extractedelements; end
+			if isfield(structmd,'nodeonboundary'), md.mesh.vertexonboundary=structmd.nodeonboundary; end
+			if isfield(structmd,'hemisphere'), md.mesh.hemisphere=structmd.hemisphere; end
+			if isfield(structmd,'lat'), md.mesh.lat=structmd.lat; end
+			if isfield(structmd,'long'), md.mesh.long=structmd.long; end
+			if isfield(structmd,'segments'), md.mesh.segments=structmd.segments; end
+			if isfield(structmd,'segmentmarkers'), md.mesh.segmentmarkers=structmd.segmentmarkers; end
+			if isfield(structmd,'dim'), md.mesh.dimension=structmd.dim; end
+			if isfield(structmd,'numlayers'), md.mesh.numberoflayers=structmd.numlayers; end
+			if isfield(structmd,'numberofelements'), md.mesh.numberofelements=structmd.numberofelements; end
+			if isfield(structmd,'numberofvertices'), md.mesh.numberofvertices=structmd.numberofvertices; end
+			if isfield(structmd,'numberofnodes'), md.mesh.numberofvertices=structmd.numberofnodes; end
+			if isfield(structmd,'numberofedges'), md.mesh.numberofedges=structmd.numberofedges; end
+			if isfield(structmd,'numberofelements2d'), md.mesh.numberofelements2d=structmd.numberofelements2d; end
+			if isfield(structmd,'numberofnodes2d'), md.mesh.numberofvertices2d=structmd.numberofnodes2d; end
+			if isfield(structmd,'nodeconnectivity'), md.mesh.vertexconnectivity=structmd.nodeconnectivity; end
+			if isfield(structmd,'elementconnectivity'), md.mesh.elementconnectivity=structmd.elementconnectivity; end
+			if isfield(structmd,'uppernodes'), md.mesh.uppervertex=structmd.uppernodes; end
+			if isfield(structmd,'lowernodes'), md.mesh.lowervertex=structmd.lowernodes; end
+			if isfield(structmd,'upperelements'), md.mesh.upperelements=structmd.upperelements; end
+			if isfield(structmd,'lowerelements'), md.mesh.lowerelements=structmd.lowerelements; end
+			if isfield(structmd,'elementonbed'), md.mesh.elementonbed=structmd.elementonbed; end
+			if isfield(structmd,'elementonsurface'), md.mesh.elementonsurface=structmd.elementonsurface; end
+			if isfield(structmd,'nodeonsurface'), md.mesh.vertexonsurface=structmd.nodeonsurface; end
+			if isfield(structmd,'nodeonbed'), md.mesh.vertexonbed=structmd.nodeonbed; end
+			if isfield(structmd,'elements2d'), md.mesh.elements2d=structmd.elements2d; end
+			if isfield(structmd,'y2d'), md.mesh.y2d=structmd.y2d; end
+			if isfield(structmd,'x2d'), md.mesh.x2d=structmd.x2d; end
+			if isfield(structmd,'elements'), md.mesh.elements=structmd.elements; end
+			if isfield(structmd,'edges'), 
+				md.mesh.edges=structmd.edges; 
+				md.mesh.edges(isnan(md.mesh.edges))=-1;
+			end
+			if isfield(structmd,'y'), md.mesh.y=structmd.y; end
+			if isfield(structmd,'x'), md.mesh.x=structmd.x; end
+			if isfield(structmd,'z'), md.mesh.z=structmd.z; end
+			if isfield(structmd,'mask'), md.flaim.criterion=structmd.mask; end
+			if isfield(structmd,'diagnostic_ref'), md.stressbalance.referential=structmd.diagnostic_ref; end
+			if isfield(structmd,'npart'); md.qmu.numberofpartitions=structmd.npart; end
+			if isfield(structmd,'part'); md.qmu.partition=structmd.part; end
+
+			%Field changes
+			if (isfield(structmd,'type') & ischar(structmd.type)), 
+				if strcmpi(structmd.type,'2d'), md.mesh.dimension=2; end
+				if strcmpi(structmd.type,'3d'), md.mesh.dimension=3; end
+			end
+			if isnumeric(md.verbose),
+				md.verbose=verbose;
+			end
+
+			if isfield(structmd,'spcvelocity'), 
+				md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+				md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+				md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+				pos=find(structmd.spcvelocity(:,1)); md.stressbalance.spcvx(pos)=structmd.spcvelocity(pos,4); 
+				pos=find(structmd.spcvelocity(:,2)); md.stressbalance.spcvy(pos)=structmd.spcvelocity(pos,5); 
+				pos=find(structmd.spcvelocity(:,3)); md.stressbalance.spcvz(pos)=structmd.spcvelocity(pos,6); 
+			end
+			if isfield(structmd,'spcvx'), 
+				md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1);
+				pos=find(~isnan(structmd.spcvx)); md.stressbalance.spcvx(pos)=structmd.spcvx(pos); 
+			end
+			if isfield(structmd,'spcvy'),
+				md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1);
+				pos=find(~isnan(structmd.spcvy)); md.stressbalance.spcvy(pos)=structmd.spcvy(pos);     
+			end
+			if isfield(structmd,'spcvz'),
+				md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1);
+				pos=find(~isnan(structmd.spcvz)); md.stressbalance.spcvz(pos)=structmd.spcvz(pos);     
+			end
+			if isfield(structmd,'pressureload'),
+				if ~isempty(structmd.pressureload) & ismember(structmd.pressureload(end,end),[118 119 120]),
+					pos=find(structmd.pressureload(:,end)==120); md.stressbalance.icefront(pos,end)=0;
+					pos=find(structmd.pressureload(:,end)==118); md.stressbalance.icefront(pos,end)=1;
+					pos=find(structmd.pressureload(:,end)==119); md.stressbalance.icefront(pos,end)=2;
+				end
+			end
+			if isfield(structmd,'elements_type') & structmd.elements_type(end,end)>50,
+				pos=find(structmd.elements_type==59); md.flowequation.element_equation(pos,end)=0;
+				pos=find(structmd.elements_type==55); md.flowequation.element_equation(pos,end)=1;
+				pos=find(structmd.elements_type==56); md.flowequation.element_equation(pos,end)=2;
+				pos=find(structmd.elements_type==60); md.flowequation.element_equation(pos,end)=3;
+				pos=find(structmd.elements_type==62); md.flowequation.element_equation(pos,end)=4;
+				pos=find(structmd.elements_type==57); md.flowequation.element_equation(pos,end)=5;
+				pos=find(structmd.elements_type==58); md.flowequation.element_equation(pos,end)=6;
+				pos=find(structmd.elements_type==61); md.flowequation.element_equation(pos,end)=7;
+			end
+			if isfield(structmd,'vertices_type') & structmd.vertices_type(end,end)>50,
+				pos=find(structmd.vertices_type==59); md.flowequation.vertex_equation(pos,end)=0;
+				pos=find(structmd.vertices_type==55); md.flowequation.vertex_equation(pos,end)=1;
+				pos=find(structmd.vertices_type==56); md.flowequation.vertex_equation(pos,end)=2;
+				pos=find(structmd.vertices_type==60); md.flowequation.vertex_equation(pos,end)=3;
+				pos=find(structmd.vertices_type==62); md.flowequation.vertex_equation(pos,end)=4;
+				pos=find(structmd.vertices_type==57); md.flowequation.vertex_equation(pos,end)=5;
+				pos=find(structmd.vertices_type==58); md.flowequation.vertex_equation(pos,end)=6;
+				pos=find(structmd.vertices_type==61); md.flowequation.vertex_equation(pos,end)=7;
+			end
+			if isfield(structmd,'rheology_law') & isnumeric(structmd.rheology_law),
+				if (structmd.rheology_law==272), md.materials.rheology_law='None';      end
+				if (structmd.rheology_law==368), md.materials.rheology_law='Paterson';  end
+				if (structmd.rheology_law==369), md.materials.rheology_law='Arrhenius'; end
+			end
+			if isfield(structmd,'groundingline_migration') & isnumeric(structmd.groundingline_migration),
+				if (structmd.groundingline_migration==272), md.groundingline.migration='None';      end
+				if (structmd.groundingline_migration==273), md.groundingline.migration='AgressiveMigration';  end
+				if (structmd.groundingline_migration==274), md.groundingline.migration='SoftMigration'; end
+			end
+			if isfield(structmd,'control_type') & isnumeric(structmd.control_type),
+				if (structmd.control_type==143), md.inversion.control_parameters={'FrictionCoefficient'}; end
+				if (structmd.control_type==190), md.inversion.control_parameters={'RheologyBbar'}; end
+				if (structmd.control_type==147), md.inversion.control_parameters={'Thickeningrate'}; end
+			end
+			if isfield(structmd,'cm_responses') & ismember(structmd.cm_responses(end,end),[165:170 383 388 389]),
+				pos=find(structmd.cm_responses==166); md.inversion.cost_functions(pos)=101;
+				pos=find(structmd.cm_responses==167); md.inversion.cost_functions(pos)=102;
+				pos=find(structmd.cm_responses==168); md.inversion.cost_functions(pos)=103;
+				pos=find(structmd.cm_responses==169); md.inversion.cost_functions(pos)=104;
+				pos=find(structmd.cm_responses==170); md.inversion.cost_functions(pos)=105;
+				pos=find(structmd.cm_responses==165); md.inversion.cost_functions(pos)=201;
+				pos=find(structmd.cm_responses==389); md.inversion.cost_functions(pos)=501;
+				pos=find(structmd.cm_responses==388); md.inversion.cost_functions(pos)=502;
+				pos=find(structmd.cm_responses==382); md.inversion.cost_functions(pos)=503;
+			end
+
+			if isfield(structmd,'artificial_diffusivity') & structmd.artificial_diffusivity==2,
+					md.thermal.stabilization=2;
+					md.masstransport.stabilization=1;
+					md.balancethickness.stabilization=1;
+			end
+			if isnumeric(md.masstransport.hydrostatic_adjustment)
+				if md.masstransport.hydrostatic_adjustment==269,
+					md.masstransport.hydrostatic_adjustment='Incremental';
+				else
+					md.masstransport.hydrostatic_adjustment='Absolute';
+				end
+			end
+
+			%New fields
+			if ~isfield(structmd,'upperelements');
+				md.mesh.upperelements=transpose(1:md.mesh.numberofelements)+md.mesh.numberofelements2d;
+				md.mesh.upperelements(end-md.mesh.numberofelements2d+1:end)=NaN;
+			end
+			if ~isfield(structmd,'lowerelements');
+				md.mesh.lowerelements=transpose(1:md.mesh.numberofelements)-md.mesh.numberofelements2d;
+				md.mesh.lowerelements(1:md.mesh.numberofelements2d)=NaN;
+			end
+			if ~isfield(structmd,'diagnostic_ref');
+				md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
+			end
+			if ~isfield(structmd,'loadingforce');
+				md.stressbalance.loadingforce=0*ones(md.mesh.numberofvertices,3);
+			end
+
+			%2013 August 9
+			if isfield(structmd,'prognostic') & isa(structmd.prognostic,'prognostic'),
+				disp('Recovering old prognostic class');
+				md.masstransport=masstransport(structmd.prognostic);
+			end
+			%2013 August 9
+			if isfield(structmd,'diagnostic') & (isa(structmd.diagnostic,'diagnostic') || isa(structmd.diagnostic,'stressbalance')),
+				disp('Recovering old diagnostic class');
+				md.stressbalance=stressbalance(structmd.diagnostic);
+			end
+		end% }}}
+		function md = setdefaultparameters(md) % {{{
+
+			%initialize subclasses
+			md.mesh             = mesh();
+			md.mask             = mask();
+			md.constants        = constants();
+			md.geometry         = geometry();
+			md.initialization   = initialization();
+			md.surfaceforcings  = surfaceforcings();
+			md.basalforcings    = basalforcings();
+			md.friction         = friction();
+			md.rifts            = rifts();
+			md.timestepping     = timestepping();
+			md.groundingline    = groundingline();
+			md.materials        = matice();
+			md.damage           = damage();
+			md.flowequation     = flowequation();
+			md.debug            = debug();
+			md.verbose          = verbose();
+			md.settings         = settings();
+			md.toolkits         = toolkits();
+			md.cluster          = generic();
+			md.balancethickness = balancethickness();
+			md.stressbalance       = stressbalance();
+			md.hydrology        = hydrologyshreve();
+			md.masstransport       = masstransport();
+			md.thermal          = thermal();
+			md.steadystate      = steadystate();
+			md.transient        = transient();
+			md.gia              = gia();
+			md.autodiff         = autodiff();
+			md.flaim            = flaim();
+			md.inversion        = inversion();
+			md.qmu              = qmu();
+			md.radaroverlay     = radaroverlay();
+			md.results          = struct();
+			md.miscellaneous    = miscellaneous();
+			md.private          = private();
+		end
+		%}}}
+		function disp(obj) % {{{
+			disp(sprintf('%19s: %-22s -- %s','mesh'            ,['[1x1 ' class(obj.mesh) ']'],'mesh properties'));
+			disp(sprintf('%19s: %-22s -- %s','mask'            ,['[1x1 ' class(obj.mask) ']'],'defines grounded and floating elements'));
+			disp(sprintf('%19s: %-22s -- %s','geometry'        ,['[1x1 ' class(obj.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
+			disp(sprintf('%19s: %-22s -- %s','constants'       ,['[1x1 ' class(obj.constants) ']'],'physical constants'));
+			disp(sprintf('%19s: %-22s -- %s','surfaceforcings' ,['[1x1 ' class(obj.surfaceforcings) ']'],'surface forcings'));
+			disp(sprintf('%19s: %-22s -- %s','basalforcings'   ,['[1x1 ' class(obj.basalforcings) ']'],'bed forcings'));
+			disp(sprintf('%19s: %-22s -- %s','materials'       ,['[1x1 ' class(obj.materials) ']'],'material properties'));
+			disp(sprintf('%19s: %-22s -- %s','damage'          ,['[1x1 ' class(obj.damage) ']'],'damage propagation laws'));
+			disp(sprintf('%19s: %-22s -- %s','friction'        ,['[1x1 ' class(obj.friction) ']'],'basal friction/drag properties'));
+			disp(sprintf('%19s: %-22s -- %s','flowequation'    ,['[1x1 ' class(obj.flowequation) ']'],'flow equations'));
+			disp(sprintf('%19s: %-22s -- %s','timestepping'    ,['[1x1 ' class(obj.timestepping) ']'],'time stepping for transient models'));
+			disp(sprintf('%19s: %-22s -- %s','initialization'  ,['[1x1 ' class(obj.initialization) ']'],'initial guess/state'));
+			disp(sprintf('%19s: %-22s -- %s','rifts'           ,['[1x1 ' class(obj.rifts) ']'],'rifts properties'));
+			disp(sprintf('%19s: %-22s -- %s','debug'           ,['[1x1 ' class(obj.debug) ']'],'debugging tools (valgrind, gprof)'));
+			disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve'));
+			disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(obj.settings) ']'],'settings properties'));
+			disp(sprintf('%19s: %-22s -- %s','toolkits'          ,['[1x1 ' class(obj.toolkits) ']'],'PETSc options for each solution'));
+			disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)'));
+			disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution'));
+			disp(sprintf('%19s: %-22s -- %s','stressbalance'      ,['[1x1 ' class(obj.stressbalance) ']'],'parameters for stressbalance solution'));
+			disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution'));
+			disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution'));
+			disp(sprintf('%19s: %-22s -- %s','masstransport'      ,['[1x1 ' class(obj.masstransport) ']'],'parameters for masstransport solution'));
+			disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution'));
+			disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution'));
+			disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution'));
+			disp(sprintf('%19s: %-22s -- %s','gia'       ,['[1x1 ' class(obj.gia) ']'],'parameters for gia solution'));
+			disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters'));
+			disp(sprintf('%19s: %-22s -- %s','flaim'           ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters'));
+			disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(obj.inversion) ']'],'parameters for inverse methods'));
+			disp(sprintf('%19s: %-22s -- %s','qmu'             ,['[1x1 ' class(obj.qmu) ']'],'dakota properties'));
+			disp(sprintf('%19s: %-22s -- %s','results'         ,['[1x1 ' class(obj.results) ']'],'model results'));
+			disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay'));
+			disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields'));
+		end % }}}
+		function memory(obj) % {{{
+
+		disp(sprintf('\nMemory imprint:\n'));
+
+		fields=properties('model');
+		mem=0;
+
+		for i=1:length(fields),
+			field=obj.(fields{i});
+			s=whos('field'); 
+			mem=mem+s.bytes/1e6;
+			disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
+		end
+		disp(sprintf('%19s--%10s','--------------','--------------'));
+		disp(sprintf('%19s: %g Mb','Total',mem));
+		end % }}}
+		function netcdf(obj,filename) % {{{
+		%NETCDF - save model as netcdf
+		%
+		%   Usage:
+		%      netcdf(md,filename)
+		%
+		%   Example:
+		%      netcdf(md,'model.nc');
+
+		disp('Saving model as NetCDF');
+		%1. Create NetCDF file
+		ncid=netcdf.create(filename,'CLOBBER');
+		netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
+		netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' obj.miscellaneous.name ')']);
+		netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
+		netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
+
+		%Preallocate variable id, needed to write variables in netcdf file
+		var_id=zeros(1000,1);%preallocate
+
+		for step=1:2,
+			counter=0;
+			[var_id,counter]=structtonc(ncid,'md',obj,0,var_id,counter,step);
+			if step==1, netcdf.endDef(ncid); end
+		end
+
+		if counter>1000,
+			warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
+		end
+
+		netcdf.close(ncid)
+		end % }}}
+		function xylim(obj) % {{{
+
+			xlim([min(obj.mesh.x) max(obj.mesh.x)]);
+			ylim([min(obj.mesh.y) max(obj.mesh.y)])
+		end % }}}
+		function md=upload(md) % {{{
+		%the goal of this routine is to upload the model onto a server, and to empty it.
+		%So first, save the model with a unique name and upload the file to the server: 
+		random_part=fix(rand(1)*10000);
+		id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part)  '-' getenv('USER') '-' oshostname() '.upload']; 
+		eval(['save ' id ' md']);
+
+		%Now, upload the file: 
+		issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
+
+		%Now, empty this model of everything except settings, and record name of file we just uploaded!
+		settings_back=md.settings;
+		md=model();
+		md.settings=settings_back;
+		md.settings.upload_filename=id;
+
+		%get locally rid of file that was uploaded
+		eval(['delete ' id]);
+
+		end % }}}
+		function md=download(md) % {{{
+
+		%the goal of this routine is to download the internals of the current model from a server, because 
+		%this model is empty, except for the settings which tell us where to go and find this model!
+
+		%Download the file: 
+		issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
+
+		name=md.settings.upload_filename;
+
+		%Now, load this model: 
+		md=loadmodel(md.settings.upload_filename);
+
+		%get locally rid of file that was downloaded
+		eval(['delete ' name]);
+
+		end % }}}
+	end
+ end
Index: /issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.py	(revision 16287)
+++ /issm/trunk-jpl/src/m/classes/model.py	(revision 16287)
@@ -0,0 +1,677 @@
+#module imports {{{
+import numpy
+import copy
+import sys
+from mesh import mesh
+from mask import mask
+from geometry import geometry
+from constants import constants
+from surfaceforcings import surfaceforcings
+from basalforcings import basalforcings
+from matice import matice
+from damage import damage
+from friction import friction
+from flowequation import flowequation
+from timestepping import timestepping
+from initialization import initialization
+from rifts import rifts
+from debug import debug
+from verbose import verbose
+from settings import settings
+from toolkits import toolkits
+from generic import generic
+from balancethickness import balancethickness
+from stressbalance import stressbalance
+from groundingline import groundingline
+from hydrologyshreve import hydrologyshreve
+from masstransport import masstransport
+from thermal import thermal
+from steadystate import steadystate
+from transient import transient
+from gia import gia
+from autodiff import autodiff
+from flaim import flaim
+from inversion import inversion
+from qmu import qmu
+from results import results
+from radaroverlay import radaroverlay
+from miscellaneous import miscellaneous
+from private import private
+from EnumDefinitions import *
+from mumpsoptions import *
+from iluasmoptions import *
+from project3d import *
+from FlagElements import *
+from NodeConnectivity import *
+from ElementConnectivity import *
+from contourenvelope import *
+from PythonFuncs import *
+#}}}
+
+class model(object):
+	#properties
+	def __init__(self):#{{{
+		self.mesh             = mesh()
+		self.mask             = mask()
+		self.geometry         = geometry()
+		self.constants        = constants()
+		self.surfaceforcings  = surfaceforcings()
+		self.basalforcings    = basalforcings()
+		self.materials        = matice()
+		self.damage           = damage()
+		self.friction         = friction()
+		self.flowequation     = flowequation()
+		self.timestepping     = timestepping()
+		self.initialization   = initialization()
+		self.rifts            = rifts()
+
+		self.debug            = debug()
+		self.verbose          = verbose()
+		self.settings         = settings()
+		self.toolkits         = toolkits()
+		self.cluster          = generic()
+
+		self.balancethickness = balancethickness()
+		self.stressbalance       = stressbalance()
+		self.groundingline    = groundingline()
+		self.hydrology        = hydrologyshreve()
+		self.masstransport       = masstransport()
+		self.thermal          = thermal()
+		self.steadystate      = steadystate()
+		self.transient        = transient()
+		self.gia              = gia()
+
+		self.autodiff         = autodiff()
+		self.flaim            = flaim()
+		self.inversion        = inversion()
+		self.qmu              = qmu()
+
+		self.results          = results()
+		self.radaroverlay     = radaroverlay()
+		self.miscellaneous    = miscellaneous()
+		self.private          = private()
+		#}}}
+	def properties(self):    # {{{
+		# ordered list of properties since vars(self) is random
+		return ['mesh',\
+		        'mask',\
+		        'geometry',\
+		        'constants',\
+		        'surfaceforcings',\
+		        'basalforcings',\
+		        'materials',\
+		        'damage',\
+		        'friction',\
+		        'flowequation',\
+		        'timestepping',\
+		        'initialization',\
+		        'rifts',\
+		        'debug',\
+		        'verbose',\
+		        'settings',\
+		        'toolkits',\
+		        'cluster',\
+		        'balancethickness',\
+		        'stressbalance',\
+		        'groundingline',\
+		        'hydrology',\
+		        'masstransport',\
+		        'thermal',\
+		        'steadystate',\
+		        'transient',\
+				  'gia',\
+		        'autodiff',\
+		        'flaim',\
+		        'inversion',\
+		        'qmu',\
+		        'results',\
+		        'radaroverlay',\
+		        'miscellaneous',\
+		        'private']
+	# }}}
+	def __repr__(obj): #{{{
+		#print "Here %s the number: %d" % ("is", 37)
+		string="%19s: %-22s -- %s" % ("mesh","[%s,%s]" % ("1x1",obj.mesh.__class__.__name__),"mesh properties")
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("mask","[%s,%s]" % ("1x1",obj.mask.__class__.__name__),"defines grounded and floating elements"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("geometry","[%s,%s]" % ("1x1",obj.geometry.__class__.__name__),"surface elevation, bedrock topography, ice thickness,..."))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("constants","[%s,%s]" % ("1x1",obj.constants.__class__.__name__),"physical constants"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("surfaceforcings","[%s,%s]" % ("1x1",obj.surfaceforcings.__class__.__name__),"surface forcings"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("basalforcings","[%s,%s]" % ("1x1",obj.basalforcings.__class__.__name__),"bed forcings"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("materials","[%s,%s]" % ("1x1",obj.materials.__class__.__name__),"material properties"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("damage","[%s,%s]" % ("1x1",obj.damage.__class__.__name__),"damage propagation laws"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("friction","[%s,%s]" % ("1x1",obj.friction.__class__.__name__),"basal friction/drag properties"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("flowequation","[%s,%s]" % ("1x1",obj.flowequation.__class__.__name__),"flow equations"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("timestepping","[%s,%s]" % ("1x1",obj.timestepping.__class__.__name__),"time stepping for transient models"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("initialization","[%s,%s]" % ("1x1",obj.initialization.__class__.__name__),"initial guess/state"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("rifts","[%s,%s]" % ("1x1",obj.rifts.__class__.__name__),"rifts properties"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("debug","[%s,%s]" % ("1x1",obj.debug.__class__.__name__),"debugging tools (valgrind, gprof)"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("verbose","[%s,%s]" % ("1x1",obj.verbose.__class__.__name__),"verbosity level in solve"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("settings","[%s,%s]" % ("1x1",obj.settings.__class__.__name__),"settings properties"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("toolkits","[%s,%s]" % ("1x1",obj.toolkits.__class__.__name__),"PETSc options for each solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("cluster","[%s,%s]" % ("1x1",obj.cluster.__class__.__name__),"cluster parameters (number of cpus...)"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("balancethickness","[%s,%s]" % ("1x1",obj.balancethickness.__class__.__name__),"parameters for balancethickness solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("stressbalance","[%s,%s]" % ("1x1",obj.stressbalance.__class__.__name__),"parameters for stressbalance solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("groundingline","[%s,%s]" % ("1x1",obj.groundingline.__class__.__name__),"parameters for groundingline solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("hydrology","[%s,%s]" % ("1x1",obj.hydrology.__class__.__name__),"parameters for hydrology solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("masstransport","[%s,%s]" % ("1x1",obj.masstransport.__class__.__name__),"parameters for masstransport solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("thermal","[%s,%s]" % ("1x1",obj.thermal.__class__.__name__),"parameters for thermal solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("steadystate","[%s,%s]" % ("1x1",obj.steadystate.__class__.__name__),"parameters for steadystate solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("transient","[%s,%s]" % ("1x1",obj.transient.__class__.__name__),"parameters for transient solution"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("autodiff","[%s,%s]" % ("1x1",obj.autodiff.__class__.__name__),"automatic differentiation parameters"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("flaim","[%s,%s]" % ("1x1",obj.flaim.__class__.__name__),"flaim parameters"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("inversion","[%s,%s]" % ("1x1",obj.inversion.__class__.__name__),"parameters for inverse methods"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("qmu","[%s,%s]" % ("1x1",obj.qmu.__class__.__name__),"dakota properties"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("results","[%s,%s]" % ("1x1",obj.results.__class__.__name__),"model results"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("radaroverlay","[%s,%s]" % ("1x1",obj.radaroverlay.__class__.__name__),"radar image for plot overlay"))
+		string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("miscellaneous","[%s,%s]" % ("1x1",obj.miscellaneous.__class__.__name__),"miscellaneous fields"))
+		return string
+	# }}}
+	def checkmessage(self,string):    # {{{
+		print ("model not consistent: %s" % string)
+		self.private.isconsistent=False
+		return self
+	# }}}
+	def extract(md,area):    # {{{
+		"""
+		extract - extract a model according to an Argus contour or flag list
+
+		   This routine extracts a submodel from a bigger model with respect to a given contour
+		   md must be followed by the corresponding exp file or flags list
+		   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
+		   If user wants every element outside the domain to be 
+		   extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
+		   an empty string '' will be considered as an empty domain
+		   a string 'all' will be considered as the entire domain
+
+		   Usage:
+		      md2=extract(md,area);
+
+		   Examples:
+		      md2=extract(md,'Domain.exp');
+
+		   See also: EXTRUDE, COLLAPSE
+		"""
+
+		#copy model
+		md1=copy.deepcopy(md)
+
+		#get elements that are inside area
+		flag_elem=FlagElements(md1,area)
+		if not numpy.any(flag_elem):
+			raise RuntimeError("extracted model is empty")
+
+		#kick out all elements with 3 dirichlets
+		spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
+		spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
+		flag=numpy.ones(md1.mesh.numberofvertices)
+		flag[spc_node]=0
+		pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]
+		flag_elem[pos]=0
+
+		#extracted elements and nodes lists
+		pos_elem=numpy.nonzero(flag_elem)[0]
+		pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
+
+		#keep track of some fields
+		numberofvertices1=md1.mesh.numberofvertices
+		numberofelements1=md1.mesh.numberofelements
+		numberofvertices2=numpy.size(pos_node)
+		numberofelements2=numpy.size(pos_elem)
+		flag_node=numpy.zeros(numberofvertices1)
+		flag_node[pos_node]=1
+
+		#Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
+		Pelem=numpy.zeros(numberofelements1,int)
+		Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
+		Pnode=numpy.zeros(numberofvertices1,int)
+		Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
+
+		#renumber the elements (some node won't exist anymore)
+		elements_1=copy.deepcopy(md1.mesh.elements)
+		elements_2=elements_1[pos_elem,:]
+		elements_2[:,0]=Pnode[elements_2[:,0]-1]
+		elements_2[:,1]=Pnode[elements_2[:,1]-1]
+		elements_2[:,2]=Pnode[elements_2[:,2]-1]
+		if md1.mesh.dimension==3:
+			elements_2[:,3]=Pnode[elements_2[:,3]-1]
+			elements_2[:,4]=Pnode[elements_2[:,4]-1]
+			elements_2[:,5]=Pnode[elements_2[:,5]-1]
+
+		#OK, now create the new model!
+
+		#take every field from model
+		md2=copy.deepcopy(md1)
+
+		#automatically modify fields
+
+		#loop over model fields
+		model_fields=vars(md1)
+		for fieldi in model_fields:
+			#get field
+			field=getattr(md1,fieldi)
+			fieldsize=numpy.shape(field)
+			if hasattr(field,'__dict__') and not ismember(fieldi,['results'])[0]:    #recursive call
+				object_fields=vars(field)
+				for fieldj in object_fields:
+					#get field
+					field=getattr(getattr(md1,fieldi),fieldj)
+					fieldsize=numpy.shape(field)
+					if len(fieldsize):
+						#size = number of nodes * n
+						if   fieldsize[0]==numberofvertices1:
+							setattr(getattr(md2,fieldi),fieldj,field[pos_node,:])
+						elif fieldsize[0]==numberofvertices1+1:
+							setattr(getattr(md2,fieldi),fieldj,numpy.vstack((field[pos_node,:],field[-1,:])))
+						#size = number of elements * n
+						elif fieldsize[0]==numberofelements1:
+							setattr(getattr(md2,fieldi),fieldj,field[pos_elem,:])
+			else:
+				if len(fieldsize):
+					#size = number of nodes * n
+					if   fieldsize[0]==numberofvertices1:
+						setattr(md2,fieldi,field[pos_node,:])
+					elif fieldsize[0]==numberofvertices1+1:
+						setattr(md2,fieldi,numpy.hstack((field[pos_node,:],field[-1,:])))
+					#size = number of elements * n
+					elif fieldsize[0]==numberofelements1:
+						setattr(md2,fieldi,field[pos_elem,:])
+
+		#modify some specific fields
+
+		#Mesh
+		md2.mesh.numberofelements=numberofelements2
+		md2.mesh.numberofvertices=numberofvertices2
+		md2.mesh.elements=elements_2
+
+		#mesh.uppervertex mesh.lowervertex
+		if md1.mesh.dimension==3:
+			md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]
+			md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1]
+
+			md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0]
+			md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1]
+
+			md2.mesh.upperelements=md1.mesh.upperelements[pos_elem]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0]
+			md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1]
+
+			md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0]
+			md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1]
+
+		#Initial 2d mesh 
+		if md1.mesh.dimension==3:
+			flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)]
+			pos_elem_2d=numpy.nonzero(flag_elem_2d)[0]
+			flag_node_2d=flag_node[numpy.arange(0,md1.mesh.numberofvertices2d)]
+			pos_node_2d=numpy.nonzero(flag_node_2d)[0]
+
+			md2.mesh.numberofelements2d=numpy.size(pos_elem_2d)
+			md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
+			md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
+			md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1]
+			md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1]
+			md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1]
+
+			md2.mesh.x2d=md1.mesh.x[pos_node_2d]
+			md2.mesh.y2d=md1.mesh.y[pos_node_2d]
+
+		#Edges
+		if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
+			#renumber first two columns
+			pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
+			md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
+			md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
+			md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2]-1]
+			md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
+			#remove edges when the 2 vertices are not in the domain.
+			md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
+			#Replace all zeros by -1 in the last two columns
+			pos=numpy.nonzero(md2.mesh.edges[:,2]==0)[0]
+			md2.mesh.edges[pos,2]=-1
+			pos=numpy.nonzero(md2.mesh.edges[:,3]==0)[0]
+			md2.mesh.edges[pos,3]=-1
+			#Invert -1 on the third column with last column (Also invert first two columns!!)
+			pos=numpy.nonzero(md2.mesh.edges[:,2]==-1)[0]
+			md2.mesh.edges[pos,2]=md2.mesh.edges[pos,3]
+			md2.mesh.edges[pos,3]=-1
+			values=md2.mesh.edges[pos,1]
+			md2.mesh.edges[pos,1]=md2.mesh.edges[pos,0]
+			md2.mesh.edges[pos,0]=values
+			#Finally remove edges that do not belong to any element
+			pos=numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,1]==-1,md2.mesh.edges[:,2]==-1))[0]
+			md2.mesh.edges=numpy.delete(md2.mesh.edges,pos,axis=0)
+
+		#Penalties
+		if numpy.any(numpy.logical_not(numpy.isnan(md2.stressbalance.vertex_pairing))):
+			for i in xrange(numpy.size(md1.stressbalance.vertex_pairing,axis=0)):
+				md2.stressbalance.vertex_pairing[i,:]=Pnode[md1.stressbalance.vertex_pairing[i,:]]
+			md2.stressbalance.vertex_pairing=md2.stressbalance.vertex_pairing[numpy.nonzero(md2.stressbalance.vertex_pairing[:,0])[0],:]
+		if numpy.any(numpy.logical_not(numpy.isnan(md2.masstransport.vertex_pairing))):
+			for i in xrange(numpy.size(md1.masstransport.vertex_pairing,axis=0)):
+				md2.masstransport.vertex_pairing[i,:]=Pnode[md1.masstransport.vertex_pairing[i,:]]
+			md2.masstransport.vertex_pairing=md2.masstransport.vertex_pairing[numpy.nonzero(md2.masstransport.vertex_pairing[:,0])[0],:]
+
+		#recreate segments
+		if md1.mesh.dimension==2:
+			[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
+			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
+			md2.mesh.segments=contourenvelope(md2)
+			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
+			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
+		else:
+			#First do the connectivity for the contourenvelope in 2d
+			[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d)
+			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
+			md2.mesh.segments=contourenvelope(md2)
+			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
+			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
+			md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
+			#Then do it for 3d as usual
+			[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
+			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
+
+		#Boundary conditions: Dirichlets on new boundary
+		#Catch the elements that have not been extracted
+		orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
+		orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
+		#Figure out which node are on the boundary between md2 and md1
+		nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
+		nodestoflag2=Pnode[nodestoflag1].astype(int)-1
+		if numpy.size(md1.stressbalance.spcvx)>1 and numpy.size(md1.stressbalance.spcvy)>2 and numpy.size(md1.stressbalance.spcvz)>2:
+			if numpy.size(md1.inversion.vx_obs)>1 and numpy.size(md1.inversion.vy_obs)>1:
+				md2.stressbalance.spcvx[nodestoflag2]=md2.inversion.vx_obs[nodestoflag2] 
+				md2.stressbalance.spcvy[nodestoflag2]=md2.inversion.vy_obs[nodestoflag2]
+			else:
+				md2.stressbalance.spcvx[nodestoflag2]=float('NaN')
+				md2.stressbalance.spcvy[nodestoflag2]=float('NaN')
+				print "\n!! extract warning: spc values should be checked !!\n\n"
+			#put 0 for vz
+			md2.stressbalance.spcvz[nodestoflag2]=0
+		if numpy.any(numpy.logical_not(numpy.isnan(md1.thermal.spctemperature))):
+			md2.thermal.spctemperature[nodestoflag2,0]=1
+
+		#Results fields
+		if md1.results:
+			md2.results=results()
+			for solutionfield,field in md1.results.__dict__.iteritems():
+				if   isinstance(field,list):
+					setattr(md2.results,solutionfield,[])
+					#get time step
+					for i,fieldi in enumerate(field):
+						if isinstance(fieldi,results) and fieldi:
+							getattr(md2.results,solutionfield).append(results())
+							fieldr=getattr(md2.results,solutionfield)[i]
+							#get subfields
+							for solutionsubfield,subfield in fieldi.__dict__.iteritems():
+								if   numpy.size(subfield)==numberofvertices1:
+									setattr(fieldr,solutionsubfield,subfield[pos_node])
+								elif numpy.size(subfield)==numberofelements1:
+									setattr(fieldr,solutionsubfield,subfield[pos_elem])
+								else:
+									setattr(fieldr,solutionsubfield,subfield)
+						else:
+							getattr(md2.results,solutionfield).append(None)
+				elif isinstance(field,results):
+					setattr(md2.results,solutionfield,results())
+					if isinstance(field,results) and field:
+						fieldr=getattr(md2.results,solutionfield)
+						#get subfields
+						for solutionsubfield,subfield in field.__dict__.iteritems():
+							if   numpy.size(subfield)==numberofvertices1:
+								setattr(fieldr,solutionsubfield,subfield[pos_node])
+							elif numpy.size(subfield)==numberofelements1:
+								setattr(fieldr,solutionsubfield,subfield[pos_elem])
+							else:
+								setattr(fieldr,solutionsubfield,subfield)
+
+		#Keep track of pos_node and pos_elem
+		md2.mesh.extractedvertices=pos_node+1
+		md2.mesh.extractedelements=pos_elem+1
+
+		return md2
+	# }}}
+	def extrude(md,*args):    # {{{
+		"""
+		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 len(args)>3 or len(args)<1:
+			raise RuntimeError("extrude error message")
+
+		#Extrude the mesh
+		if   len(args)==1:    #list of coefficients
+			clist=args[0]
+			if any(clist<0) or any(clist>1):
+				raise TypeError("extrusioncoefficients must be between 0 and 1")
+			clist.extend([0.,1.])
+			clist.sort()
+			extrusionlist=list(set(clist))
+			numlayers=len(extrusionlist)
+
+		elif len(args)==2:    #one polynomial law
+			if args[1]<=0:
+				raise TypeError("extrusionexponent must be >=0")
+			numlayers=args[0]
+			extrusionlist=(numpy.arange(0.,float(numlayers-1)+1.,1.)/float(numlayers-1))**args[1]
+
+		elif len(args)==3:    #two polynomial laws
+			numlayers=args[0]
+			lowerexp=args[1]
+			upperexp=args[2]
+
+			if args[1]<=0 or args[2]<=0:
+				raise TypeError("lower and upper extrusionexponents must be >=0")
+
+			lowerextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**lowerexp/2.
+			upperextrusionlist=(numpy.arange(0.,1.+2./float(numlayers-1),2./float(numlayers-1)))**upperexp/2.
+			extrusionlist=numpy.unique(numpy.concatenate((lowerextrusionlist,1.-upperextrusionlist)))
+
+		if numlayers<2:
+			raise TypeError("number of layers should be at least 2")
+		if md.mesh.dimension==3:
+			raise TypeError("Cannot extrude a 3d mesh (extrude cannot be called more than once)")
+
+		#Initialize with the 2d mesh
+		x3d=numpy.empty((0))
+		y3d=numpy.empty((0))
+		z3d=numpy.empty((0))    #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 in xrange(numlayers):
+			x3d=numpy.concatenate((x3d,md.mesh.x))
+			y3d=numpy.concatenate((y3d,md.mesh.y))
+			#nodes are distributed between bed and surface accordingly to the given exponent
+			z3d=numpy.concatenate((z3d,(bed3d+thickness3d*extrusionlist[i]).reshape(-1)))
+		number_nodes3d=numpy.size(x3d)    #number of 3d nodes for the non extruded part of the mesh
+
+		#Extrude elements 
+		elements3d=numpy.empty((0,6),int)
+		for i in xrange(numlayers-1):
+			elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices))))    #Create the elements of the 3d mesh for the non extruded part
+		number_el3d=numpy.size(elements3d,axis=0)    #number of 3d nodes for the non extruded part of the mesh
+
+		#Keep a trace of lower and upper nodes
+		mesh.lowervertex=-1*numpy.ones(number_nodes3d,int)
+		mesh.uppervertex=-1*numpy.ones(number_nodes3d,int)
+		mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
+		mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
+		md.mesh.lowervertex=mesh.lowervertex
+		md.mesh.uppervertex=mesh.uppervertex
+
+		#same for lower and upper elements
+		mesh.lowerelements=-1*numpy.ones(number_el3d,int)
+		mesh.upperelements=-1*numpy.ones(number_el3d,int)
+		mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
+		mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
+		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 not numpy.any(numpy.isnan(md.initialization.vx)):
+			md.initialization.vx=project3d(md,'vector',md.initialization.vx,'type','node')
+		if not numpy.any(numpy.isnan(md.initialization.vy)):
+			md.initialization.vy=project3d(md,'vector',md.initialization.vy,'type','node')
+		if not numpy.any(numpy.isnan(md.initialization.vz)):
+			md.initialization.vz=project3d(md,'vector',md.initialization.vz,'type','node')
+		if not numpy.any(numpy.isnan(md.initialization.vel)):
+			md.initialization.vel=project3d(md,'vector',md.initialization.vel,'type','node')
+		if not numpy.any(numpy.isnan(md.initialization.temperature)):
+			md.initialization.temperature=project3d(md,'vector',md.initialization.temperature,'type','node')
+		if not numpy.any(numpy.isnan(md.initialization.waterfraction)):
+			md.initialization.waterfraction=project3d(md,'vector',md.initialization.waterfraction,'type','node')
+
+		#bedinfo and surface info
+		md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1)
+		md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1)
+		md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
+		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
+
+		#elementstype
+		if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
+			oldelements_type=md.flowequation.element_equation
+			md.flowequation.element_equation=numpy.zeros(number_el3d,int)
+			md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
+
+		#verticestype
+		if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
+			oldvertices_type=md.flowequation.vertex_equation
+			md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int)
+			md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
+
+		md.flowequation.borderSSA=project3d(md,'vector',md.flowequation.borderSSA,'type','node')
+		md.flowequation.borderHO=project3d(md,'vector',md.flowequation.borderHO,'type','node')
+		md.flowequation.borderFS=project3d(md,'vector',md.flowequation.borderFS,'type','node')
+
+		#boundary conditions
+		md.stressbalance.spcvx=project3d(md,'vector',md.stressbalance.spcvx,'type','node')
+		md.stressbalance.spcvy=project3d(md,'vector',md.stressbalance.spcvy,'type','node')
+		md.stressbalance.spcvz=project3d(md,'vector',md.stressbalance.spcvz,'type','node')
+		md.thermal.spctemperature=project3d(md,'vector',md.thermal.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',float('NaN'))
+		md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node')
+		md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node')
+		md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node')
+		md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node')
+		md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
+
+		#connectivity
+		md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
+		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
+		for i in xrange(1,numlayers-1):
+			md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \
+				=md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
+		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=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')
+
+		#damage
+		md.damage.D=project3d(md,'vector',md.damage.D,'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.gia.mantle_viscosity=project3d(md,'vector',md.gia.mantle_viscosity,'type','node')
+		md.gia.lithosphere_thickness=project3d(md,'vector',md.gia.lithosphere_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.ice_levelset=project3d(md,'vector',md.mask.ice_levelset,'type','node')
+		md.mask.groundedice_levelset=project3d(md,'vector',md.mask.groundedice_levelset,'type','node')
+		if not numpy.any(numpy.isnan(md.inversion.cost_functions_coefficients)):
+			md.inversion.cost_functions_coefficients=project3d(md,'vector',md.inversion.cost_functions_coefficients,'type','node');end;
+		if not numpy.any(numpy.isnan(md.inversion.min_parameters)):
+			md.inversion.min_parameters=project3d(md,'vector',md.inversion.min_parameters,'type','node')
+		if not numpy.any(numpy.isnan(md.inversion.max_parameters)):
+			md.inversion.max_parameters=project3d(md,'vector',md.inversion.max_parameters,'type','node')
+		if not numpy.any(numpy.isnan(md.qmu.partition)):
+			md.qmu.partition=project3d(md,'vector',numpy.transpose(md.qmu.partition),'type','node')
+		if(md.surfaceforcings.isdelta18o):
+			md.surfaceforcings.temperatures_lgm=project3d(md,'vector',md.surfaceforcings.temperatures_lgm,'type','node')
+		if(md.surfaceforcings.isdelta18o):
+			md.surfaceforcings.temperatures_presentday=project3d(md,'vector',md.surfaceforcings.temperatures_presentday,'type','node')
+		if(md.surfaceforcings.isdelta18o):
+			md.surfaceforcings.precipitations_presentday=project3d(md,'vector',md.surfaceforcings.precipitations_presentday,'type','node')
+
+		#Put lithostatic pressure if there is an existing pressure
+		if not numpy.any(numpy.isnan(md.initialization.pressure)):
+			md.initialization.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z.reshape(-1,1))
+
+		#special for thermal modeling:
+		md.basalforcings.melting_rate=project3d(md,'vector',md.basalforcings.melting_rate,'type','node','layer',1)
+		if not numpy.any(numpy.isnan(md.basalforcings.geothermalflux)):
+			md.basalforcings.geothermalflux=project3d(md,'vector',md.basalforcings.geothermalflux,'type','node','layer',1)    #bedrock only gets geothermal flux
+
+		#increase connectivity if less than 25:
+		if md.mesh.average_vertex_connectivity<=25:
+			md.mesh.average_vertex_connectivity=100
+
+		return md
+		# }}}
Index: sm/trunk-jpl/src/m/classes/planetmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/planetmesh.m	(revision 16286)
+++ 	(revision )
@@ -1,123 +1,0 @@
-%PLANETMESH class definition
-%
-%   Usage:
-%      planetmesh=planetmesh();
-
-classdef planetmesh
-	properties (SetAccess=public) 
-		x                           = NaN;
-		y                           = NaN;
-		z                           = NaN;
-		r                           = NaN;
-		theta                       = NaN;
-		phi                         = NaN
-		elements                    = NaN
-		dimension                   = 0;
-		numberoflayers              = 0;
-		numberofelements            = 0;
-		numberofvertices            = 0;
-
-		vertexconnectivity          = NaN
-		elementconnectivity         = NaN
-		average_vertex_connectivity = 0;
-	end
-	methods
-		function obj = planetmesh(varargin) % {{{
-			switch nargin
-				case 0
-					obj=setdefaultparameters(obj);
-				otherwise
-					error('constructor not supported');
-			end
-		end % }}}
-		function obj = setdefaultparameters(obj) % {{{
-
-			%the connectivity is the avergaded number of nodes linked to a
-			%given node through an edge. This connectivity is used to initially
-			%allocate memory to the stiffness matrix. A value of 16 seems to
-			%give a good memory/time ration. This value can be checked in
-			%trunk/test/Miscellaneous/runme.m
-			obj.average_vertex_connectivity=25;
-		end % }}}
-		function md = checkconsistency(obj,md,solution,analyses) % {{{
-
-			md = checkfield(md,'planetmesh.x','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.y','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.z','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.r','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.theta','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.phi','NaN',1,'size',[md.planetmesh.numberofvertices 1]);
-			md = checkfield(md,'planetmesh.elements','NaN',1,'>',0,'values',1:md.planetmesh.numberofvertices);
-			if(md.planetmesh.dimension==2),
-				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 3]);
-			else
-				md = checkfield(md,'planetmesh.elements','size',[md.planetmesh.numberofelements 6]);
-			end
-			if any(~ismember(1:md.planetmesh.numberofvertices,sort(unique(md.planetmesh.elements(:)))));
-				md = checkmessage(md,'orphan nodes have been found. Check the planetmesh outline');
-			end
-			md = checkfield(md,'planetmesh.dimension','values',[2 3]);
-			md = checkfield(md,'planetmesh.numberoflayers','>=',0);
-			md = checkfield(md,'planetmesh.numberofelements','>',0);
-			md = checkfield(md,'planetmesh.numberofvertices','>',0);
-			if (md.planetmesh.dimension==2),
-				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',9,'message','''planetmesh.average_vertex_connectivity'' should be at least 9 in 2d');
-			else
-				md = checkfield(md,'planetmesh.average_vertex_connectivity','>=',24,'message','''planetmesh.average_vertex_connectivity'' should be at least 24 in 3d');
-			end
-			md = checkfield(md,'planetmesh.elementconnectivity','size',[md.planetmesh.numberofelements 3],'NaN',1);
-
-			%Solution specific checks
-			switch(solution),
-				case MasstransportSolutionEnum(),
-					if md.masstransport.stabilization==3,
-						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
-					end
-				case TransientSolutionEnum(),
-					if md.transient.ismasstransport & md.masstransport.stabilization==3,
-						md = checkfield(md,'planetmesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d planetmeshes');
-					end
-				case ThermalSolutionEnum(),
-					md = checkfield(md,'planetmesh.dimension','values',3,'message','thermal solution only supported on 3d planetmeshes');
-			end
-		end % }}}
-		function disp(obj) % {{{
-			disp(sprintf('   Mesh:')); 
-
-			disp(sprintf('\n      Elements and vertices:'));
-			fielddisplay(obj,'numberofelements','number of elements');
-			fielddisplay(obj,'numberofvertices','number of vertices');
-			fielddisplay(obj,'elements','vertex indices of the mesh elements');
-			fielddisplay(obj,'x','vertices x coordinate [m]');
-			fielddisplay(obj,'y','vertices y coordinate [m]');
-			fielddisplay(obj,'z','vertices z coordinate [m]');
-			fielddisplay(obj,'r','vertices r coordinate [m]');
-			fielddisplay(obj,'theta','vertices theta coordinate [degrees]');
-			fielddisplay(obj,'phi','vertices phi coordinate [degrees]');
-
-			disp(sprintf('\n      Properties:'));
-			fielddisplay(obj,'dimension','planetmesh dimension (2d or 3d)');
-			fielddisplay(obj,'numberoflayers','number of extrusion layers');
-
-			fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i');
-			fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i');
-			fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
-
-		end % }}}
-		function marshall(obj,md,fid) % {{{
-			WriteData(fid,'object',obj,'fieldname','x','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','y','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','z','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','r','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','theta','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','phi','format','DoubleMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
-			WriteData(fid,'object',obj,'fieldname','dimension','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','numberoflayers','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','elementconnectivity','format','DoubleMat','mattype',3);
-			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
-		end % }}}
-	end
-end
Index: /issm/trunk-jpl/src/m/classes/spheremesh.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/spheremesh.m	(revision 16287)
+++ /issm/trunk-jpl/src/m/classes/spheremesh.m	(revision 16287)
@@ -0,0 +1,123 @@
+%SPHEREMESH class definition
+%
+%   Usage:
+%      spheremesh=spheremesh();
+
+classdef spheremesh
+	properties (SetAccess=public) 
+		x                           = NaN;
+		y                           = NaN;
+		z                           = NaN;
+		r                           = NaN;
+		theta                       = NaN;
+		phi                         = NaN
+		elements                    = NaN
+		dimension                   = 0;
+		numberoflayers              = 0;
+		numberofelements            = 0;
+		numberofvertices            = 0;
+
+		vertexconnectivity          = NaN
+		elementconnectivity         = NaN
+		average_vertex_connectivity = 0;
+	end
+	methods
+		function obj = spheremesh(varargin) % {{{
+			switch nargin
+				case 0
+					obj=setdefaultparameters(obj);
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function obj = setdefaultparameters(obj) % {{{
+
+			%the connectivity is the avergaded number of nodes linked to a
+			%given node through an edge. This connectivity is used to initially
+			%allocate memory to the stiffness matrix. A value of 16 seems to
+			%give a good memory/time ration. This value can be checked in
+			%trunk/test/Miscellaneous/runme.m
+			obj.average_vertex_connectivity=25;
+		end % }}}
+		function md = checkconsistency(obj,md,solution,analyses) % {{{
+
+			md = checkfield(md,'spheremesh.x','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.y','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.z','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.r','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.theta','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.phi','NaN',1,'size',[md.spheremesh.numberofvertices 1]);
+			md = checkfield(md,'spheremesh.elements','NaN',1,'>',0,'values',1:md.spheremesh.numberofvertices);
+			if(md.spheremesh.dimension==2),
+				md = checkfield(md,'spheremesh.elements','size',[md.spheremesh.numberofelements 3]);
+			else
+				md = checkfield(md,'spheremesh.elements','size',[md.spheremesh.numberofelements 6]);
+			end
+			if any(~ismember(1:md.spheremesh.numberofvertices,sort(unique(md.spheremesh.elements(:)))));
+				md = checkmessage(md,'orphan nodes have been found. Check the spheremesh outline');
+			end
+			md = checkfield(md,'spheremesh.dimension','values',[2 3]);
+			md = checkfield(md,'spheremesh.numberoflayers','>=',0);
+			md = checkfield(md,'spheremesh.numberofelements','>',0);
+			md = checkfield(md,'spheremesh.numberofvertices','>',0);
+			if (md.spheremesh.dimension==2),
+				md = checkfield(md,'spheremesh.average_vertex_connectivity','>=',9,'message','''spheremesh.average_vertex_connectivity'' should be at least 9 in 2d');
+			else
+				md = checkfield(md,'spheremesh.average_vertex_connectivity','>=',24,'message','''spheremesh.average_vertex_connectivity'' should be at least 24 in 3d');
+			end
+			md = checkfield(md,'spheremesh.elementconnectivity','size',[md.spheremesh.numberofelements 3],'NaN',1);
+
+			%Solution specific checks
+			switch(solution),
+				case MasstransportSolutionEnum(),
+					if md.masstransport.stabilization==3,
+						md = checkfield(md,'spheremesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d spheremeshes');
+					end
+				case TransientSolutionEnum(),
+					if md.transient.ismasstransport & md.masstransport.stabilization==3,
+						md = checkfield(md,'spheremesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d spheremeshes');
+					end
+				case ThermalSolutionEnum(),
+					md = checkfield(md,'spheremesh.dimension','values',3,'message','thermal solution only supported on 3d spheremeshes');
+			end
+		end % }}}
+		function disp(obj) % {{{
+			disp(sprintf('   Mesh:')); 
+
+			disp(sprintf('\n      Elements and vertices:'));
+			fielddisplay(obj,'numberofelements','number of elements');
+			fielddisplay(obj,'numberofvertices','number of vertices');
+			fielddisplay(obj,'elements','vertex indices of the mesh elements');
+			fielddisplay(obj,'x','vertices x coordinate [m]');
+			fielddisplay(obj,'y','vertices y coordinate [m]');
+			fielddisplay(obj,'z','vertices z coordinate [m]');
+			fielddisplay(obj,'r','vertices r coordinate [m]');
+			fielddisplay(obj,'theta','vertices theta coordinate [degrees]');
+			fielddisplay(obj,'phi','vertices phi coordinate [degrees]');
+
+			disp(sprintf('\n      Properties:'));
+			fielddisplay(obj,'dimension','spheremesh dimension (2d or 3d)');
+			fielddisplay(obj,'numberoflayers','number of extrusion layers');
+
+			fielddisplay(obj,'vertexconnectivity','list of vertices connected to vertex_i');
+			fielddisplay(obj,'elementconnectivity','list of vertices connected to element_i');
+			fielddisplay(obj,'average_vertex_connectivity','average number of vertices connected to one vertex');
+
+		end % }}}
+		function marshall(obj,md,fid) % {{{
+			WriteData(fid,'object',obj,'fieldname','x','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','y','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','z','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','r','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','theta','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','phi','format','DoubleMat','mattype',1);
+			WriteData(fid,'object',obj,'fieldname','elements','format','DoubleMat','mattype',2);
+			WriteData(fid,'object',obj,'fieldname','dimension','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberoflayers','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberofelements','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','numberofvertices','format','Integer');
+			WriteData(fid,'object',obj,'fieldname','elementconnectivity','format','DoubleMat','mattype',3);
+			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
+		end % }}}
+	end
+end