Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m	(revision 17687)
@@ -52,8 +52,10 @@
 
 %First find segments that are not completely on the front
-if ~strcmp(md.mesh.domaintype(),'3D'),
+if strcmp(elementtype(md.mesh),'Penta'),
+	numbernodesfront=4;
+elseif strcmp(elementtype(md.mesh),'Tria'),
 	numbernodesfront=2;
-else 
-	numbernodesfront=4;
+else
+	error('mesh type not supported yet');
 end
 segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0;
Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17686)
+++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 17687)
@@ -53,8 +53,10 @@
 
 	#First find segments that are not completely on the front
-	if not m.strcmp(md.mesh.domaintype(),'3D'):
+	if m.strcmp(md.mesh.elementtype(),'Penta'):
+		numbernodesfront=4
+	elif m.strcmp(md.mesh.elementtype(),'Tria'):
 		numbernodesfront=2
 	else:
-		numbernodesfront=4
+			raise StandardError("Mesh type not supported")
 	if any(md.mask.ice_levelset<=0):
 		values=md.mask.ice_levelset[md.mesh.segments[:,0:-1]-1]
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17687)
@@ -156,5 +156,5 @@
 
 			%Check that the model is really a 3d model
-			if ~strcmp(md.mesh.domaintype(),'3D'),
+			if ~strcmp(md.mesh.elementtype(),'Penta'),
 				error('collapse error message: only 3d mesh can be collapsed')
 			end
Index: /issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/contrib/ecco/PropagateFlagsUntilDistance.m	(revision 17687)
@@ -10,5 +10,5 @@
 
 %make 3d work in 2d: 
-if strcmp(domaintype(md.mesh),'3D'),
+if dimension(md.mesh)==3,
 	md.mesh.x=md.mesh.x2d;
 	md.mesh.y=md.mesh.y2d;
Index: /issm/trunk-jpl/src/m/extrusion/DepthAverage.m
===================================================================
--- /issm/trunk-jpl/src/m/extrusion/DepthAverage.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/extrusion/DepthAverage.m	(revision 17687)
@@ -9,5 +9,5 @@
 
 %check that the model given in input is 3d
-if ~strcmp(md.mesh.domaintype,'3D');
+if ~strcmp(md.mesh.elementtype(),'3D');
 	error('DepthAverage error message: the model given in input must be 3d')
 end
Index: /issm/trunk-jpl/src/m/parameterization/contourenvelope.m
===================================================================
--- /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.m	(revision 17687)
@@ -42,5 +42,5 @@
 %get nodes inside profile
 mesh.elementconnectivity=md.mesh.elementconnectivity;
-if strcmp(md.mesh.domaintype(),'2Dhorizontal'),
+if dimension(md.mesh)==2,
 	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 17686)
+++ /issm/trunk-jpl/src/m/parameterization/contourenvelope.py	(revision 17687)
@@ -47,5 +47,5 @@
 	#get nodes inside profile
 	elementconnectivity=copy.deepcopy(md.mesh.elementconnectivity)
-	if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
+	if md.mesh.dimension()==2:
 		elements=copy.deepcopy(md.mesh.elements)
 		x=copy.deepcopy(md.mesh.x)
Index: /issm/trunk-jpl/src/m/plot/applyoptions.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.m	(revision 17687)
@@ -49,5 +49,5 @@
 
 %view 
-if strcmp(domaintype(md.mesh),'3D') & ~exist(options,'layer'),
+if dimension(md.mesh)==3 & ~exist(options,'layer'),
 	view(getfieldvalue(options,'view',3));
 else
@@ -60,5 +60,5 @@
 	eval(['axis ' getfieldvalue(options,'axis')]);
 else
-	if strcmpi(domaintype(md.mesh),'2Dhorizontal') | exist(options,'layer'),
+	if dimension(md.mesh)==2 | exist(options,'layer'),
 		axis tight equal;
 	else
Index: /issm/trunk-jpl/src/m/plot/plot_edges.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_edges.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_edges.m	(revision 17687)
@@ -17,5 +17,5 @@
 end
 
-if (strcmp(domaintype(md.mesh),'2Dhorizontal')),
+if dimension(md.mesh)==2,
 	%plot mesh
 	A=elements(:,1); B=elements(:,2); C=elements(:,3); 
Index: /issm/trunk-jpl/src/m/plot/plot_icefront.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_icefront.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_icefront.m	(revision 17687)
@@ -19,5 +19,5 @@
 elementzeroice=sum(zeroice(md.mesh.elements),2);
 
-if (strcmp(domaintype(md.mesh),'2Dhorizontal')),
+if dimension(md.mesh)==2,
 	icefront=(elementice & elementnoice) & ~(elementice==2 & elementzeroice);
 
Index: /issm/trunk-jpl/src/m/plot/plot_penalties.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_penalties.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_penalties.m	(revision 17687)
@@ -21,5 +21,5 @@
 end
 
-if ~strcmp(domaintype(md.mesh),'3D'),
+if dimension(md.mesh)~=3,
 	error('no penalties to plot for ''2d'' model');
 elseif isempty(md.penalties),
Index: /issm/trunk-jpl/src/m/plot/plot_qmu_mass_flux_segments.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_qmu_mass_flux_segments.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_qmu_mass_flux_segments.m	(revision 17687)
@@ -13,5 +13,5 @@
 allsegments=md.qmu.mass_flux_segments;
 
-if (strcmp(domaintype(md.mesh),'2Dhorizontal')),
+if dimension(md.mesh)==2,
 
 	%recover segments
Index: /issm/trunk-jpl/src/m/plot/plot_referential.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_referential.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_referential.m	(revision 17687)
@@ -28,5 +28,5 @@
 Yhat=cross(Zhat,Xhat);
 
-if (strcmp(domaintype(md.mesh),'2Dhorizontal')),
+if dimension(md.mesh)==2,
 
 	%plot mesh
Index: /issm/trunk-jpl/src/m/plot/plot_segments.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_segments.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_segments.m	(revision 17687)
@@ -14,5 +14,5 @@
 segments=md.mesh.segments;
 
-if (strcmp(domaintype(md.mesh),'2Dhorizontal')),
+if dimension(md.mesh)==2,
 	%plot mesh
 	A=elements(:,1); B=elements(:,2); C=elements(:,3); 
Index: /issm/trunk-jpl/src/m/plot/plot_tensor_components.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_tensor_components.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_tensor_components.m	(revision 17687)
@@ -10,5 +10,5 @@
 upperplots=fix((i-1)/width);
 if upperplots==0, leftplots=i-1; else leftplots=i-width*upperplots-1; end
-if (strcmp(domaintype(md.mesh),'2Dhorizontal'))%3 components -> 3 indexes
+if dimension(md.mesh)==2 %3 components -> 3 indexes
 	index1=4*width*upperplots+2*leftplots+1;
 	index2=index1+1;
@@ -34,5 +34,5 @@
 end
 
-if ((strcmp(domaintype(md.mesh),'2Dhorizontal'))),
+if dimension(md.mesh)==2,
 	subplot(2*width,2*width,index1),
 	plot_unit(x,y,z,elements,tensor.xx,is2d,isplanet,datatype,options)
Index: /issm/trunk-jpl/src/m/plot/plot_transient_results.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_transient_results.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/plot_transient_results.m	(revision 17687)
@@ -41,5 +41,5 @@
 clear string;
 
-if strcmp(domaintype(md.mesh),'3D'),
+if dimension(md.mesh)==3,
 	string='plotmodel(md';
 	for i=1:length(md.results.transient),
@@ -67,5 +67,5 @@
 clear string;
 
-if strcmp(domaintype(md.mesh),'3D'),
+if dimension(md.mesh)==3,
 	string='plotmodel(md';
 	for i=2:length(md.results.transient),
Index: /issm/trunk-jpl/src/m/plot/processdata.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/processdata.m	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/processdata.m	(revision 17687)
@@ -65,5 +65,5 @@
 
 	%check number of columns, add zeros if necessary,
-	if (strcmp(md.mesh.domaintype(),'3D'))
+	if (dimension(md.mesh)==3)
 		if datasize(2)==2,
 			data=[data, zeros(datasize(1),1)];
@@ -85,5 +85,5 @@
 
 %treat the case datasize(1)=nodes2d
-if (strcmp(md.mesh.domaintype(),'3D') & datasize(1)==numberofvertices2d),
+if (dimension(md.mesh)==3 & datasize(1)==numberofvertices2d),
 	data=project3d(md,'vector',data,'type','node');
 	datasize(1)=md.mesh.numberofvertices;
@@ -92,5 +92,5 @@
 
 %treat the case datasize(1)=nodes2d
-if (strcmp(md.mesh.domaintype(),'3D') & datasize(1)==numberofelements2d),
+if (dimension(md.mesh)==3 & datasize(1)==numberofelements2d),
 	data=project3d(md,'vector',data,'type','element');
 	datasize(1)=md.mesh.numberofelements;
Index: /issm/trunk-jpl/src/m/plot/processmesh.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 17686)
+++ /issm/trunk-jpl/src/m/plot/processmesh.py	(revision 17687)
@@ -42,13 +42,11 @@
 
 		#is it a 2D plot?
-		if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
+		if md.mesh.dimension()==2:
 			is2d=1
-		elif m.strcmp(md.mesh.domaintype(),'3D'):
+		else:
 			if options.getfieldvalue('layer',0)>=1:
 				is2d=1
 			else:
 				is2d=0
-		else:
-			raise ValueError('processmesh error: dim = %d not supported' % md.mesh.dimension)
 
 		#layer projection?
@@ -64,5 +62,5 @@
 	else:
 		#Process mesh for plotting 
-		if m.strcmp(md.mesh.domaintype(),'2Dhorizontal'):
+		if md.mesh.dimension()==2:
 			is2d=1
 		else: