Index: sm/trunk-jpl/src/m/classes/mesh.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh.m	(revision 17557)
+++ 	(revision )
@@ -1,210 +1,0 @@
-%MESH class definition
-%
-%   Usage:
-%      mesh=mesh();
-
-classdef mesh
-	properties (SetAccess=public) 
-		x                           = NaN;
-		y                           = NaN;
-		z                           = NaN
-		elements                    = NaN
-		dimension                   = 0;
-		numberoflayers              = 0;
-		numberofelements            = 0;
-		numberofvertices            = 0;
-		numberofedges               = 0;
-
-		lat                         = NaN
-		long                        = NaN
-		hemisphere                  = NaN
-
-		elementonbed                = NaN
-		elementonsurface            = NaN
-		vertexonbed                 = NaN
-		vertexonsurface             = NaN
-		lowerelements               = NaN
-		lowervertex                 = NaN
-		upperelements               = NaN
-		uppervertex                 = NaN
-		vertexonboundary            = NaN
-
-		edges                       = NaN
-		segments                    = NaN
-		segmentmarkers              = NaN
-		vertexconnectivity          = NaN
-		elementconnectivity         = NaN
-		average_vertex_connectivity = 0;
-
-		x2d                         = NaN
-		y2d                         = NaN
-		elements2d                  = NaN
-		numberofvertices2d          = 0;
-		numberofelements2d          = 0;
-
-		extractedvertices           = NaN
-		extractedelements           = NaN
-	end
-	methods (Static)
-		function obj = loadobj(obj) % {{{
-			% This function is directly called by matlab when a model object is
-			% loaded. Update old properties here
-
-			%2012 June 28th
-			if numel(obj.edges)>1 & any(isnan(obj.edges(:)))
-				disp('Update model edges from previous version');
-				obj.edges(isnan(obj.edges))=-1;
-			end
-
-		end% }}}
-	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,'fieldname','mesh.x','NaN',1,'size',[md.mesh.numberofvertices 1]);
-			md = checkfield(md,'fieldname','mesh.y','NaN',1,'size',[md.mesh.numberofvertices 1]);
-			md = checkfield(md,'fieldname','mesh.z','NaN',1,'size',[md.mesh.numberofvertices 1]);
-			md = checkfield(md,'fieldname','mesh.elements','NaN',1,'>',0,'values',1:md.mesh.numberofvertices);
-			if(md.mesh.dimension==2),
-				md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements 3]);
-			else
-				md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements 6]);
-			end
-			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,'fieldname','mesh.dimension','values',[2 3]);
-			md = checkfield(md,'fieldname','mesh.numberoflayers','>=',0);
-			md = checkfield(md,'fieldname','mesh.numberofelements','>',0);
-			md = checkfield(md,'fieldname','mesh.numberofvertices','>',0);
-			%no checks for numberofedges lat long and hemisphere
-			md = checkfield(md,'fieldname','mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]);
-			md = checkfield(md,'fieldname','mesh.elementonsurface','size',[md.mesh.numberofelements 1],'values',[0 1]);
-			md = checkfield(md,'fieldname','mesh.vertexonbed','size',[md.mesh.numberofvertices 1],'values',[0 1]);
-			md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices 1],'values',[0 1]);
-			if (md.mesh.dimension==3),
-				md = checkfield(md,'fieldname','mesh.z','>=',md.geometry.bed-10^-10,'message','''mesh.z'' lower than bedrock');
-				md = checkfield(md,'fieldname','mesh.z','<=',md.geometry.surface+10^-10,'message','''mesh.z'' higher than surface elevation');
-			end
-			if (md.mesh.dimension==2),
-				md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message','''mesh.average_vertex_connectivity'' should be at least 9 in 2d');
-			else
-				md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message','''mesh.average_vertex_connectivity'' should be at least 24 in 3d');
-			end
-
-			%Solution specific checks
-			switch(solution),
-				case MasstransportSolutionEnum(),
-					if md.masstransport.stabilization==3,
-						md = checkfield(md,'fieldname','mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
-					end
-				case BalancethicknessSolutionEnum(),
-					if md.balancethickness.stabilization==3,
-						md = checkfield(md,'fieldname','mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
-					end
-				case TransientSolutionEnum(),
-					if md.transient.ismasstransport & md.masstransport.stabilization==3,
-						md = checkfield(md,'fieldname','mesh.dimension','values',2,'message','Discontinuous Galerkin only supported for 2d meshes');
-					end
-				case ThermalSolutionEnum(),
-					md = checkfield(md,'fieldname','mesh.dimension','values',3,'message','thermal solution only supported on 3d meshes');
-			end
-		end % }}}
-		function disp(obj) % {{{
-			disp(sprintf('   Mesh:')); 
-
-			if obj.dimension==3,
-				disp(sprintf('\n      Elements and vertices of the original 2d mesh:'));
-				fielddisplay(obj,'numberofelements2d','number of elements');
-				fielddisplay(obj,'numberofvertices2d','number of vertices');
-				fielddisplay(obj,'elements2d','vertex indices of the mesh elements');
-				fielddisplay(obj,'x2d','vertices x coordinate [m]');
-				fielddisplay(obj,'y2d','vertices y coordinate [m]');
-
-				disp(sprintf('\n      Elements and vertices of the extruded 3d mesh:'));
-			else
-				disp(sprintf('\n      Elements and vertices:'));
-			end
-			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,'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,'dimension','mesh dimension');
-			fielddisplay(obj,'numberoflayers','number of extrusion layers');
-			fielddisplay(obj,'vertexonbed','lower vertices flags list');
-			fielddisplay(obj,'elementonbed','lower elements flags list');
-			fielddisplay(obj,'vertexonsurface','upper vertices flags list');
-			fielddisplay(obj,'elementonsurface','upper elements flags list');
-			fielddisplay(obj,'uppervertex','upper vertex list (NaN for vertex on the upper surface)');
-			fielddisplay(obj,'upperelements','upper element list (NaN for element on the upper layer)');
-			fielddisplay(obj,'lowervertex','lower vertex list (NaN for vertex on the lower surface)');
-			fielddisplay(obj,'lowerelements','lower element list (NaN for element on the lower layer');
-			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,'enum',MeshTypeEnum(),'data',StringToEnum(['Mesh' meshtype(obj)]),'format','Integer');
-			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','elements','format','DoubleMat','mattype',2);
-			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','elementonbed','format','BooleanMat','mattype',2);
-			WriteData(fid,'object',obj,'fieldname','elementonsurface','format','BooleanMat','mattype',2);
-			WriteData(fid,'object',obj,'fieldname','vertexonbed','format','BooleanMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','vertexonsurface','format','BooleanMat','mattype',1);
-			WriteData(fid,'object',obj,'fieldname','lowerelements','format','DoubleMat','mattype',2);
-			WriteData(fid,'object',obj,'fieldname','upperelements','format','DoubleMat','mattype',2);
-			WriteData(fid,'object',obj,'fieldname','average_vertex_connectivity','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','elements2d','format','DoubleMat','mattype',3);
-			WriteData(fid,'object',obj,'fieldname','numberofvertices2d','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','numberofelements2d','format','Integer');
-		end % }}}
-		function type = meshtype(obj) % {{{
-			if obj.dimension==2,
-				type = '2Dhorizontal';
-			else
-				type = '3D';
-			end
-		end % }}}
-	end
-end
Index: sm/trunk-jpl/src/m/classes/mesh.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh.py	(revision 17557)
+++ 	(revision )
@@ -1,191 +1,0 @@
-import numpy
-from fielddisplay import fielddisplay
-from EnumDefinitions import *
-from checkfield import checkfield
-import MatlabFuncs as m
-from WriteData import WriteData
-
-class mesh(object):
-	"""
-	MESH class definition
-
-	   Usage:
-	      mesh=mesh();
-	"""
-
-	def __init__(self): # {{{
-		self.x                           = float('NaN');
-		self.y                           = float('NaN');
-		self.z                           = float('NaN');
-		self.elements                    = float('NaN');
-		self.dimension                   = 0;
-		self.numberoflayers              = 0;
-		self.numberofelements            = 0;
-		self.numberofvertices            = 0;
-		self.numberofedges               = 0;
-		
-		self.lat                         = float('NaN');
-		self.long                        = float('NaN');
-		self.hemisphere                  = float('NaN');
-
-		self.elementonbed                = float('NaN');
-		self.elementonsurface            = float('NaN');
-		self.vertexonbed                 = float('NaN');
-		self.vertexonsurface             = float('NaN');
-		self.lowerelements               = float('NaN');
-		self.lowervertex                 = float('NaN');
-		self.upperelements               = float('NaN');
-		self.uppervertex                 = 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.x2d                         = float('NaN');
-		self.y2d                         = float('NaN');
-		self.elements2d                  = float('NaN');
-		self.numberofvertices2d          = 0;
-		self.numberofelements2d          = 0;
-
-		self.extractedvertices           = float('NaN');
-		self.extractedelements           = float('NaN');
-
-		#set defaults
-		self.setdefaultparameters()
-		#}}}
-	def __repr__(self): # {{{
-		string="   Mesh:" 
-
-
-		if self.dimension==3:
-			string="%s\n%s"%(string,"\n      Elements and vertices of the original 2d mesh:")
-			
-			string="%s\n%s"%(string,fielddisplay(self,"numberofelements2d","number of elements"))
-			string="%s\n%s"%(string,fielddisplay(self,"numberofvertices2d","number of vertices"))
-			string="%s\n%s"%(string,fielddisplay(self,"elements2d","vertex indices of the mesh elements"))
-			string="%s\n%s"%(string,fielddisplay(self,"x2d","vertices x coordinate [m]"))
-			string="%s\n%s"%(string,fielddisplay(self,"y2d","vertices y coordinate [m]"))
-
-			string="%s\n%s"%(string,"\n\n      Elements and vertices of the extruded 3d mesh:")
-		else:
-			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,"z","vertices z 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,"dimension","mesh dimension (2d or 3d)"))
-		string="%s\n%s"%(string,fielddisplay(self,"numberoflayers","number of extrusion layers"))
-		string="%s\n%s"%(string,fielddisplay(self,"vertexonbed","lower vertices flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"elementonbed","lower elements flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (-1 for vertex on the upper surface)"))
-		string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (-1 for element on the upper layer)"))
-		string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (-1 for vertex on the lower surface)"))
-		string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (-1 for element on the lower layer)"))
-		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,'fieldname','mesh.x','NaN',1,'size',[md.mesh.numberofvertices])
-		md = checkfield(md,'fieldname','mesh.y','NaN',1,'size',[md.mesh.numberofvertices])
-		md = checkfield(md,'fieldname','mesh.z','NaN',1,'size',[md.mesh.numberofvertices])
-		md = checkfield(md,'fieldname','mesh.elements','NaN',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1))
-		if md.mesh.dimension==2:
-			md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3])
-		else:
-			md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,6])
-		if numpy.any(numpy.logical_not(m.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,'fieldname','mesh.dimension','values',[2,3])
-		md = checkfield(md,'fieldname','mesh.numberoflayers','>=',0)
-		md = checkfield(md,'fieldname','mesh.numberofelements','>',0)
-		md = checkfield(md,'fieldname','mesh.numberofvertices','>',0)
-		#no checks for numberofedges lat long and hemisphere
-		md = checkfield(md,'fieldname','mesh.elementonbed','size',[md.mesh.numberofelements],'values',[0,1])
-		md = checkfield(md,'fieldname','mesh.elementonsurface','size',[md.mesh.numberofelements],'values',[0,1])
-		md = checkfield(md,'fieldname','mesh.vertexonbed','size',[md.mesh.numberofvertices],'values',[0,1])
-		md = checkfield(md,'fieldname','mesh.vertexonsurface','size',[md.mesh.numberofvertices],'values',[0,1])
-		if md.mesh.dimension==2:
-			md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',9,'message',"'mesh.average_vertex_connectivity' should be at least 9 in 2d")
-		else:
-			md = checkfield(md,'fieldname','mesh.average_vertex_connectivity','>=',24,'message',"'mesh.average_vertex_connectivity' should be at least 24 in 3d")
-
-		#Solution specific checks
-		if   solution==MasstransportSolutionEnum():
-			if md.masstransport.stabilization==3:
-				md = checkfield(md,'fieldname','mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes")
-		elif solution==BalancethicknessSolutionEnum():
-			if md.balancethickness.stabilization==3:
-				md = checkfield(md,'fieldname','mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes")
-		elif solution==TransientSolutionEnum():
-			if md.transient.ismasstransport and md.masstransport.stabilization==3:
-				md = checkfield(md,'fieldname','mesh.dimension','values',2,'message',"Discontinuous Galerkin only supported for 2d meshes")
-		elif solution==ThermalSolutionEnum():
-			md = checkfield(md,'fieldname','mesh.dimension','values',3,'message','thermal solution only supported on 3d meshes')
-
-		return md
-	# }}}
-	def meshtype(self): # {{{
-
-		if self.dimension==2:
-			return "2Dhorizontal"
-		else:
-			return "3D"
-	#}}}
-	def marshall(self,md,fid):    # {{{
-		WriteData(fid,'enum',MeshTypeEnum(),'data',StringToEnum("Mesh"+self.meshtype())[0],'format','Integer');
-		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','z','format','DoubleMat','mattype',1)
-		WriteData(fid,'object',self,'fieldname','elements','format','DoubleMat','mattype',2)
-		WriteData(fid,'object',self,'fieldname','numberoflayers','format','Integer')
-		WriteData(fid,'object',self,'fieldname','numberofelements','format','Integer')
-		WriteData(fid,'object',self,'fieldname','numberofvertices','format','Integer')
-		WriteData(fid,'object',self,'fieldname','elementonbed','format','BooleanMat','mattype',2)
-		WriteData(fid,'object',self,'fieldname','elementonsurface','format','BooleanMat','mattype',2)
-		WriteData(fid,'object',self,'fieldname','vertexonbed','format','BooleanMat','mattype',1)
-		WriteData(fid,'object',self,'fieldname','vertexonsurface','format','BooleanMat','mattype',1)
-		WriteData(fid,'object',self,'fieldname','lowerelements','format','DoubleMat','mattype',2)
-		WriteData(fid,'object',self,'fieldname','upperelements','format','DoubleMat','mattype',2)
-		WriteData(fid,'object',self,'fieldname','average_vertex_connectivity','format','Integer')
-		WriteData(fid,'object',self,'fieldname','elements2d','format','DoubleMat','mattype',3)
-		WriteData(fid,'object',self,'fieldname','numberofvertices2d','format','Integer')
-		WriteData(fid,'object',self,'fieldname','numberofelements2d','format','Integer')
-	# }}}
Index: /issm/trunk-jpl/src/m/classes/mesh2d.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh2d.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/classes/mesh2d.py	(revision 17558)
@@ -4,4 +4,5 @@
 from checkfield import checkfield
 import MatlabFuncs as m
+from WriteData import WriteData
 
 class mesh2d(object):
@@ -85,5 +86,4 @@
 		md = checkfield(md,'fieldname','mesh.x','NaN',1,'size',[md.mesh.numberofvertices])
 		md = checkfield(md,'fieldname','mesh.y','NaN',1,'size',[md.mesh.numberofvertices])
-		md = checkfield(md,'fieldname','mesh.z','NaN',1,'size',[md.mesh.numberofvertices])
 		md = checkfield(md,'fieldname','mesh.elements','NaN',1,'>',0,'values',numpy.arange(1,md.mesh.numberofvertices+1))
 		md = checkfield(md,'fieldname','mesh.elements','size',[md.mesh.numberofelements,3])
Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dprisms.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.m	(revision 17558)
@@ -46,4 +46,14 @@
 				case 0
 					obj=setdefaultparameters(obj);
+				case 1
+					self=mesh3dprisms();
+					object=varargin{1};
+					fields=fieldnames(object);
+					for i=1:length(fields)
+						field=fields{i};
+						if ismember(field,properties('mesh3dprisms')),
+							self.(field)=object.(field);
+						end
+					end
 				otherwise
 					error('constructor not supported');
Index: /issm/trunk-jpl/src/m/classes/mesh3dprisms.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/mesh3dprisms.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/classes/mesh3dprisms.py	(revision 17558)
@@ -4,4 +4,5 @@
 from checkfield import *
 from MatlabFuncs import *
+from WriteData import WriteData
 
 class mesh3dprisms(object):
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17558)
@@ -110,4 +110,12 @@
 				md.surfaceforcings=SMB();
 				md.surfaceforcings.mass_balance=mass_balance;
+			end
+			%2014 March 26th
+			if isa(md.mesh,'mesh'),
+				if(md.mesh.dimension==2),
+					md.mesh=mesh2d(md.mesh);
+				else
+					md.mesh=mesh3dprisms(md.mesh);
+				end
 			end
 		end% }}}
@@ -143,5 +151,5 @@
 
 			%Check that the model is really a 3d model
-			if ~md.mesh.dimension==3,
+			if ~strcmp(md.mesh.meshtype(),'3D'),
 				error('collapse error message: only 3d mesh can be collapsed')
 			end
@@ -179,10 +187,4 @@
 			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)
@@ -237,4 +239,5 @@
 
 			%Initialize with the 2d mesh
+			md.mesh=mesh2d();
 			md.mesh.x=md.mesh.x2d;
 			md.mesh.y=md.mesh.y2d;
@@ -258,6 +261,4 @@
 			md.mesh.numberoflayers=0;
 
-			%Update mesh type
-			md.mesh.dimension=2;
 		end % }}}
 		function md2 = extract(md,area) % {{{
@@ -327,5 +328,5 @@
 			elements_2(:,2)=Pnode(elements_2(:,2));
 			elements_2(:,3)=Pnode(elements_2(:,3));
-			if md1.mesh.dimension==3,
+			if isa(md1.mesh,'mesh3dprisms'),
 				elements_2(:,4)=Pnode(elements_2(:,4));
 				elements_2(:,5)=Pnode(elements_2(:,5));
@@ -383,5 +384,5 @@
 
 			%mesh.uppervertex mesh.lowervertex
-			if md1.mesh.dimension==3
+			if isa(md1.mesh,'mesh3dprisms'),
 				md2.mesh.uppervertex=md1.mesh.uppervertex(pos_node);
 				pos=find(~isnan(md2.mesh.uppervertex));
@@ -402,5 +403,5 @@
 
 			%Initial 2d mesh 
-			if md1.mesh.dimension==3
+			if isa(md1.mesh,'mesh3dprisms'),
 				flag_elem_2d=flag_elem(1:md1.mesh.numberofelements2d);
 				pos_elem_2d=find(flag_elem_2d);
@@ -461,5 +462,5 @@
 
 			%recreate segments
-			if md1.mesh.dimension==2
+			if isa(md1.mesh,'mesh2d'),
 				md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
 				md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
@@ -600,9 +601,29 @@
 				error('number of layers should be at least 2');
 			end
-			if md.mesh.dimension==3,
+			if strcmp(md.mesh.meshtype(),'3D')
 				error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
 			end
 
 			%Initialize with the 2d mesh
+			mesh2d = md.mesh;
+			md.mesh=mesh3dprisms();
+			md.mesh.x                           = mesh2d.x;
+			md.mesh.y                           = mesh2d.y;
+			md.mesh.elements                    = mesh2d.elements;
+			md.mesh.numberofelements            = mesh2d.numberofelements;
+			md.mesh.numberofvertices            = mesh2d.numberofvertices;
+
+			md.mesh.lat                         = mesh2d.lat;
+			md.mesh.long                        = mesh2d.long;
+			md.mesh.hemisphere                  = mesh2d.hemisphere;
+
+			md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
+			md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
+			md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
+			md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
+
+			md.mesh.extractedvertices           = mesh2d.extractedvertices;
+			md.mesh.extractedelements           = mesh2d.extractedelements;
+
 			x3d=[]; 
 			y3d=[];
@@ -649,7 +670,4 @@
 			md.mesh.numberofelements2d=md.mesh.numberofelements;
 			md.mesh.numberofvertices2d=md.mesh.numberofvertices;
-
-			%Update mesh type
-			md.mesh.dimension=3;
 
 			%Build global 3d mesh 
@@ -906,5 +924,4 @@
 			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
@@ -940,9 +957,4 @@
 			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;
@@ -1066,5 +1078,5 @@
 
 			%initialize subclasses
-			md.mesh             = mesh();
+			md.mesh             = mesh2d();
 			md.mask             = mask();
 			md.constants        = constants();
Index: /issm/trunk-jpl/src/m/classes/model.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/classes/model.py	(revision 17558)
@@ -4,4 +4,5 @@
 import sys
 from mesh import mesh
+from mesh3dprisms import mesh3dprisms
 from mask import mask
 from geometry import geometry
@@ -239,5 +240,5 @@
 		elements_2[:,1]=Pnode[elements_2[:,1]-1]
 		elements_2[:,2]=Pnode[elements_2[:,2]-1]
-		if md1.mesh.dimension==3:
+		if md1.mesh.__class__.__name__=='mesh3dprisms':
 			elements_2[:,3]=Pnode[elements_2[:,3]-1]
 			elements_2[:,4]=Pnode[elements_2[:,4]-1]
@@ -291,5 +292,5 @@
 
 		#mesh.uppervertex mesh.lowervertex
-		if md1.mesh.dimension==3:
+		if md1.mesh.__class__.__name__=='mesh3dprisms':
 			md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
 			pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]
@@ -309,5 +310,5 @@
 
 		#Initial 2d mesh 
-		if md1.mesh.dimension==3:
+		if md1.mesh.__class__.__name__=='mesh3dprisms':
 			flag_elem_2d=flag_elem[numpy.arange(0,md1.mesh.numberofelements2d)]
 			pos_elem_2d=numpy.nonzero(flag_elem_2d)[0]
@@ -362,5 +363,5 @@
 
 		#recreate segments
-		if md1.mesh.dimension==2:
+		if md1.mesh.__class__.__name__=='mesh2d':
 			[md2.mesh.vertexconnectivity]=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices)
 			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
@@ -498,8 +499,28 @@
 		if numlayers<2:
 			raise TypeError("number of layers should be at least 2")
-		if md.mesh.dimension==3:
+		if md.mesh.__class__.__name__=='mesh3dprisms':
 			raise TypeError("Cannot extrude a 3d mesh (extrude cannot be called more than once)")
 
 		#Initialize with the 2d mesh
+		mesh2d = md.mesh
+		md.mesh=mesh3dprisms()
+		md.mesh.x                           = mesh2d.x
+		md.mesh.y                           = mesh2d.y
+		md.mesh.elements                    = mesh2d.elements
+		md.mesh.numberofelements            = mesh2d.numberofelements
+		md.mesh.numberofvertices            = mesh2d.numberofvertices
+
+		md.mesh.lat                         = mesh2d.lat
+		md.mesh.long                        = mesh2d.long
+		md.mesh.hemisphere                  = mesh2d.hemisphere
+
+		md.mesh.vertexonboundary            = mesh2d.vertexonboundary
+		md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity
+		md.mesh.elementconnectivity         = mesh2d.elementconnectivity
+		md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity
+
+		md.mesh.extractedvertices           = mesh2d.extractedvertices
+		md.mesh.extractedelements           = mesh2d.extractedelements
+
 		x3d=numpy.empty((0))
 		y3d=numpy.empty((0))
@@ -544,7 +565,4 @@
 		md.mesh.numberofelements2d=md.mesh.numberofelements
 		md.mesh.numberofvertices2d=md.mesh.numberofvertices
-
-		#Update mesh type
-		md.mesh.dimension=3
 
 		#Build global 3d mesh 
Index: /issm/trunk-jpl/src/m/classes/oldclasses/mesh.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/oldclasses/mesh.m	(revision 17558)
+++ /issm/trunk-jpl/src/m/classes/oldclasses/mesh.m	(revision 17558)
@@ -0,0 +1,43 @@
+classdef mesh
+	properties (SetAccess=public) 
+		x                           = NaN;
+		y                           = NaN;
+		z                           = NaN
+		elements                    = NaN
+		dimension                   = 0;
+		numberoflayers              = 0;
+		numberofelements            = 0;
+		numberofvertices            = 0;
+		numberofedges               = 0;
+
+		lat                         = NaN
+		long                        = NaN
+		hemisphere                  = NaN
+
+		elementonbed                = NaN
+		elementonsurface            = NaN
+		vertexonbed                 = NaN
+		vertexonsurface             = NaN
+		lowerelements               = NaN
+		lowervertex                 = NaN
+		upperelements               = NaN
+		uppervertex                 = NaN
+		vertexonboundary            = NaN
+
+		edges                       = NaN
+		segments                    = NaN
+		segmentmarkers              = NaN
+		vertexconnectivity          = NaN
+		elementconnectivity         = NaN
+		average_vertex_connectivity = 0;
+
+		x2d                         = NaN
+		y2d                         = NaN
+		elements2d                  = NaN
+		numberofvertices2d          = 0;
+		numberofelements2d          = 0;
+
+		extractedvertices           = NaN
+		extractedelements           = NaN
+	end
+end
Index: /issm/trunk-jpl/src/m/mesh/argusmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/argusmesh.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/argusmesh.m	(revision 17558)
@@ -77,15 +77,10 @@
 
 %Finally, use model constructor to build a complete model: 
+md.mesh=mesh2d();
 md.mesh.elements=elements;
 md.mesh.x=x;
 md.mesh.y=y;
-md.z=z;
 md.mesh.numberofvertices=size(md.mesh.x,1);
 md.mesh.numberofelements=size(md.mesh.elements,1);
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
-md.mesh.dimension=2;
 md=addnote(md,notes);
 
Index: /issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 17558)
@@ -333,5 +333,5 @@
 
 if getfieldvalue(options,'vertical',0),
-	md.mesh=mesh2dvertical;
+	md.mesh=mesh2dvertical();
 	md.mesh.x=bamgmesh_out.Vertices(:,1);
 	md.mesh.y=bamgmesh_out.Vertices(:,2);
@@ -350,5 +350,5 @@
 	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 else
-	% plug results onto model
+	md.mesh=mesh2d();
 	md.mesh.x=bamgmesh_out.Vertices(:,1);
 	md.mesh.y=bamgmesh_out.Vertices(:,2);
@@ -359,15 +359,11 @@
 
 	%Fill in rest of fields:
-	md.mesh.dimension=2;
 	md.mesh.numberofelements=size(md.mesh.elements,1);
 	md.mesh.numberofvertices=length(md.mesh.x);
 	md.mesh.numberofedges=size(md.mesh.edges,1);
-	md.mesh.z=zeros(md.mesh.numberofvertices,1);
-	md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-	md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-	md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-	md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
 	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 end
+
+%Bamg private fields
 md.private.bamg=struct();
 md.private.bamg.mesh=bamgmesh(bamgmesh_out);
@@ -375,5 +371,4 @@
 md.mesh.elementconnectivity=md.private.bamg.mesh.ElementConnectivity;
 md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
-
 
 %Check for orphan
Index: /issm/trunk-jpl/src/m/mesh/bamg.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/bamg.py	(revision 17558)
@@ -1,4 +1,5 @@
 import os.path
 import numpy
+from mesh2d import mesh2d
 from collections import OrderedDict
 from pairoptions import pairoptions
@@ -320,4 +321,5 @@
 	md.private.bamg['mesh']=bamgmesh(bamgmesh_out)
 	md.private.bamg['geometry']=bamggeom(bamggeom_out)
+	md.mesh = mesh2d()
 	md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
 	md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
@@ -328,13 +330,7 @@
 
 	#Fill in rest of fields:
-	md.mesh.dimension=2
 	md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
 	md.mesh.numberofvertices=numpy.size(md.mesh.x)
 	md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
-	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
 	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
 	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
Index: /issm/trunk-jpl/src/m/mesh/meshconvert.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshconvert.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/meshconvert.m	(revision 17558)
@@ -36,12 +36,6 @@
 
 %Fill in rest of fields:
-md.mesh.dimension        = 2;
 md.mesh.numberofelements = size(md.mesh.elements,1);
 md.mesh.numberofvertices = length(md.mesh.x);
 md.mesh.numberofedges    = size(md.mesh.edges,1);
-md.mesh.z                = zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonbed      = ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface  = ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed     = ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface = ones(md.mesh.numberofelements,1);
 md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2)) = 1;
Index: /issm/trunk-jpl/src/m/mesh/meshconvert.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshconvert.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/meshconvert.py	(revision 17558)
@@ -2,4 +2,5 @@
 from collections import OrderedDict
 from BamgConvertMesh import BamgConvertMesh 
+from mesh2d   import mesh2d
 from bamgmesh import bamgmesh
 from bamggeom import bamggeom
@@ -33,4 +34,5 @@
 	md.private.bamg['mesh']     = bamgmesh(bamgmesh_out)
 	md.private.bamg['geometry'] = bamggeom(bamggeom_out)
+	md.mesh                     = mesh2d()
 	md.mesh.x                   = bamgmesh_out['Vertices'][:,0].copy()
 	md.mesh.y                   = bamgmesh_out['Vertices'][:,1].copy()
@@ -41,14 +43,7 @@
 
 	#Fill in rest of fields:
-	md.mesh.dimension          = 2
 	md.mesh.numberofelements   = numpy.size(md.mesh.elements,axis=0)
 	md.mesh.numberofvertices   = numpy.size(md.mesh.x)
 	md.mesh.numberofedges      = numpy.size(md.mesh.edges,axis=0)
-	md.mesh.z                  = numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonbed        = numpy.ones(md.mesh.numberofvertices,bool)
-	md.mask.vertexonwater      = numpy.zeros(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonsurface    = numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.elementonbed       = numpy.ones(md.mesh.numberofelements,bool)
-	md.mesh.elementonsurface   = numpy.ones(md.mesh.numberofelements,bool)
 	md.mesh.vertexonboundary   = numpy.zeros(md.mesh.numberofvertices,bool)
 	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshaddrifts.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshaddrifts.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshaddrifts.m	(revision 17558)
@@ -80,13 +80,6 @@
 %finish up "a la" mesh.h
 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
 
 %Now, build the connectivity tables for this mesh.
 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-
-%type of model
-md.mesh.dimension=2;
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m	(revision 17558)
@@ -82,11 +82,5 @@
 md.mesh.numberofelements=length(md.mesh.elements);
 md.mesh.numberofvertices=length(md.mesh.x);
-md.mesh.z=zeros(md.mesh.numberofvertices,1);
 md.mesh.vertexonboundary=zeros(length(md.mesh.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-md.flowequation.element_equation=3*ones(md.mesh.numberofelements,1);
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
 end
 
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py	(revision 17558)
@@ -80,13 +80,7 @@
 	md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
 	md.mesh.numberofvertices=numpy.size(md.mesh.x)
-	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
 	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
 	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
 	md.rifts.numrifts=length(md.rifts.riftstruct)
-	md.flowequation.element_equation=3*numpy.ones(md.mesh.numberofelements,int)
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
 
 	return md
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.m	(revision 17558)
@@ -34,10 +34,5 @@
 md.mesh.numberofelements=length(md.mesh.elements);
 md.mesh.numberofvertices=length(md.mesh.x);
-md.mesh.z=zeros(md.mesh.numberofvertices,1);
 md.mesh.vertexonboundary=zeros(length(md.mesh.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
 
 %get coordinates of rift tips
Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py	(revision 17558)
@@ -35,11 +35,6 @@
 	md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0)
 	md.mesh.numberofvertices=numpy.size(md.mesh.x)
-	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
 	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
 	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
 
 	#get coordinates of rift tips
Index: /issm/trunk-jpl/src/m/mesh/squaremesh.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/squaremesh.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/squaremesh.m	(revision 17558)
@@ -55,11 +55,9 @@
 
 %plug coordinates and nodes
+md.mesh=mesh2d();
 md.mesh.x=x;
 md.mesh.y=y;
-md.mesh.z=zeros(nods,1);
 md.mesh.numberofvertices=nods;
 md.mesh.vertexonboundary=zeros(nods,1);md.mesh.vertexonboundary(segments(:,1:2))=1;
-md.mesh.vertexonbed=ones(nods,1);
-md.mesh.vertexonsurface=ones(nods,1);
 
 %plug elements
@@ -67,11 +65,6 @@
 md.mesh.segments=segments;
 md.mesh.numberofelements=nel;
-md.mesh.elementonbed=ones(nel,1);
-md.mesh.elementonsurface=ones(nel,1);
 
 %Now, build the connectivity tables for this mesh.
 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-
-%plug other field
-md.mesh.dimension=2;
Index: /issm/trunk-jpl/src/m/mesh/squaremesh.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/squaremesh.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/squaremesh.py	(revision 17558)
@@ -57,12 +57,10 @@
 
 	#plug coordinates and nodes
+	md.mesh=mesh2d()
 	md.mesh.x=x
 	md.mesh.y=y
-	md.mesh.z=numpy.zeros((nods))
 	md.mesh.numberofvertices=nods
 	md.mesh.vertexonboundary=numpy.zeros((nods),bool)
 	md.mesh.vertexonboundary[segments[:,0:2]-1]=True
-	md.mesh.vertexonbed=numpy.ones((nods),bool)
-	md.mesh.vertexonsurface=numpy.ones((nods),bool)
 
 	#plug elements
@@ -70,6 +68,4 @@
 	md.mesh.segments=segments
 	md.mesh.numberofelements=nel
-	md.mesh.elementonbed=numpy.ones((nel),bool)
-	md.mesh.elementonsurface=numpy.ones((nel),bool)
 
 	#Now, build the connectivity tables for this mesh.
@@ -77,7 +73,3 @@
 	[md.mesh.elementconnectivity]=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity)
 
-	#plug other field
-	md.mesh.dimension=2
-
 	return md
-
Index: /issm/trunk-jpl/src/m/mesh/triangle.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/triangle.m	(revision 17558)
@@ -59,4 +59,5 @@
 
 %plug into md
+md.mesh=mesh2d();
 md.mesh.x=x;
 md.mesh.y=y;
@@ -68,15 +69,7 @@
 md.mesh.numberofelements=size(md.mesh.elements,1);
 md.mesh.numberofvertices=length(md.mesh.x);
-md.mesh.z=zeros(md.mesh.numberofvertices,1);
 md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
 
 %Now, build the connectivity tables for this mesh.
 md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-
-%type of model
-md.mesh.dimension=2;
Index: /issm/trunk-jpl/src/m/mesh/triangle.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/triangle.py	(revision 17558)
@@ -1,3 +1,4 @@
 import numpy
+from mesh2d import mesh2d
 from TriMesh import TriMesh
 from NodeConnectivity import NodeConnectivity
@@ -43,4 +44,5 @@
 
 	#Mesh using TriMesh
+	md.mesh=mesh2d()
 	[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area)
 	md.mesh.elements=md.mesh.elements.astype(int)
@@ -51,11 +53,6 @@
 	md.mesh.numberofelements = numpy.size(md.mesh.elements,axis=0)
 	md.mesh.numberofvertices = numpy.size(md.mesh.x)
-	md.mesh.z = numpy.zeros(md.mesh.numberofvertices)
 	md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool)
 	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
-	md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool)
-	md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool)
-	md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool)
 
 	#Now, build the connectivity tables for this mesh.
@@ -63,7 +60,3 @@
 	[md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
 
-	#type of model
-	md.mesh.dimension = 2.
-
 	return md
-
Index: /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 17558)
@@ -46,5 +46,5 @@
 
 %plug into md
-md.mesh=mesh2dvertical;
+md.mesh=mesh2dvertical();
 md.mesh.x=x;
 md.mesh.z=z;
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 17557)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 17558)
@@ -42,5 +42,5 @@
 %get nodes inside profile
 mesh.elementconnectivity=md.mesh.elementconnectivity;
-if md.mesh.dimension==2;
+if strcmp(md.mesh.meshtype(),'2Dhorizontal'),
 	mesh.elements=md.mesh.elements;
 	mesh.x=md.mesh.x;
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.py
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.py	(revision 17557)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py	(revision 17558)
@@ -46,5 +46,5 @@
 	#get nodes inside profile
 	mesh.elementconnectivity=copy.deepcopy(md.mesh.elementconnectivity)
-	if md.mesh.dimension==2:
+	if m.strcmp(md.mesh.meshtype(),'2Dhorizontal'):
 		mesh.elements=copy.deepcopy(md.mesh.elements)
 		mesh.x=copy.deepcopy(md.mesh.x)