Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17773)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17774)
@@ -1210,4 +1210,13 @@
 					error('Case not supported'); %not supposed to happen!
 				end
+				%Reorder elements to make sure they are direct
+				for j=1:3
+					element = elements(3*(i-1)+j,:);
+					matrix = [md.mesh.x(element), md.mesh.y(element), md.mesh.z(element), ones(4,1)];
+					if det(matrix)>0,
+						elements(3*(i-1)+j,1)=element(2);
+						elements(3*(i-1)+j,2)=element(1);
+					end
+				end
 			end
 			%%Split in 3 tetras
@@ -1216,5 +1225,10 @@
 			%subelement3 = [5 6 3 1];
 			%elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
-			disp([num2str(steiner) ' Steiner points had to be included'])
+			if steiner==0,
+				disp('No Steiner point required to split prismatic mesh into tets');
+			else
+				disp([num2str(steiner) ' Steiner points had to be included'])
+				error('Steiner point not supported yet');
+			end
 
 			pos_elements = repmat([1:md.mesh.numberofelements]',3,1);
@@ -1224,6 +1238,8 @@
 
 			%p and q (same deal, except for element that are on the bedrock: )
-			md.friction.p=md.friction.p(pos_elements);
-			md.friction.q=md.friction.q(pos_elements);
+			if ~isnan(md.friction.p),
+				md.friction.p=md.friction.p(pos_elements);
+				md.friction.q=md.friction.q(pos_elements);
+			end
 
 			%elementstype
@@ -1237,5 +1253,7 @@
 
 			%materials
-			md.materials.rheology_n=md.materials.rheology_n(pos_elements);
+			if ~isnan(md.materials.rheology_n),
+				md.materials.rheology_n=md.materials.rheology_n(pos_elements);
+			end
 
 			%increase connectivity if less than 25: