Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 17753)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 17754)
@@ -1144,9 +1144,78 @@
 			md.mesh=mesh3dtetras(md.mesh);
 
-			%Split in 3 tetras
-			subelement1 = [1 2 3 5];
-			subelement2 = [4 6 5 1];
-			subelement3 = [5 6 3 1];
-			elements=[md.mesh.elements(:,subelement1);md.mesh.elements(:,subelement2);md.mesh.elements(:,subelement3)];
+			%Subdivision from Philipp Furnstahl (http://studierstube.icg.tugraz.at/thesis/fuernstahl_thesis.pdf)
+			steiner  = 0;
+			nbv      = md.mesh.numberofvertices;
+			nbt      = 3*md.mesh.numberofelements;
+			elements = zeros(nbt,4);
+			for i=1:md.mesh.numberofelements
+				v1=md.mesh.elements(i,1); v2=md.mesh.elements(i,2); v3=md.mesh.elements(i,3);
+				v4=md.mesh.elements(i,4); v5=md.mesh.elements(i,5); v6=md.mesh.elements(i,6);
+				if(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
+					steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
+					md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
+					md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
+					md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
+					elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
+					elements(3*(i-1)+2,:) = [v1 v2 v4 v7];
+					elements(3*(i-1)+3,:) = [v2 v4 v5 v7];
+					elements(end+1,:) = [v2 v3 v5 v7];
+					elements(end+1,:) = [v3 v5 v6 v7];
+					elements(end+1,:) = [v1 v3 v6 v7];
+					elements(end+1,:) = [v1 v4 v6 v7];
+					elements(end+1,:) = [v4 v5 v6 v7];
+				elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v2 v4 v6];
+					elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
+					elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
+				elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
+					elements(3*(i-1)+2,:) = [v2 v3 v4 v5];
+					elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
+				elseif(min(v2,v4)<min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v2 v3 v4];
+					elements(3*(i-1)+2,:) = [v2 v4 v5 v6];
+					elements(3*(i-1)+3,:) = [v2 v3 v4 v6];
+				elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
+					elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+					elements(3*(i-1)+3,:) = [v1 v3 v5 v6];
+				elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)>min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v4 v5 v6];
+					elements(3*(i-1)+2,:) = [v1 v2 v5 v6];
+					elements(3*(i-1)+3,:) = [v1 v2 v3 v6];
+				elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)<min(v2,v6)),
+					elements(3*(i-1)+1,:) = [v1 v3 v4 v5];
+					elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+					elements(3*(i-1)+3,:) = [v3 v4 v5 v6];
+				elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)<min(v3,v4) & min(v3,v5)<min(v2,v6)),
+					%NEW
+					elements(3*(i-1)+1,:) = [v1 v5 v6 v4];
+					elements(3*(i-1)+2,:) = [v1 v2 v3 v5];
+					elements(3*(i-1)+3,:) = [v5 v6 v3 v1];
+				elseif(min(v2,v4)>min(v1,v5) & min(v1,v6)>min(v3,v4) & min(v3,v5)>min(v2,v6)),
+					steiner = steiner+1; nbv = nbv+1; nbt = nbt+5; v7 = nbv;
+					md.mesh.x=[md.mesh.x; mean(md.mesh.x(md.mesh.elements(i,:)))];
+					md.mesh.y=[md.mesh.y; mean(md.mesh.y(md.mesh.elements(i,:)))];
+					md.mesh.z=[md.mesh.z; mean(md.mesh.z(md.mesh.elements(i,:)))];
+					elements(3*(i-1)+1,:) = [v1 v2 v3 v7];
+					elements(3*(i-1)+2,:) = [v1 v4 v5 v7];
+					elements(3*(i-1)+3,:) = [v1 v2 v5 v7];
+					elements(end+1,:) = [v2 v5 v6 v7];
+					elements(end+1,:) = [v2 v3 v6 v7];
+					elements(end+1,:) = [v3 v4 v6 v7];
+					elements(end+1,:) = [v1 v3 v4 v7];
+					elements(end+1,:) = [v4 v5 v6 v7];
+				else
+					error('Case not supported'); %not supposed to happen!
+				end
+			end
+			%%Split in 3 tetras
+			%subelement1 = [1 2 3 5];
+			%subelement2 = [4 6 5 1];
+			%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'])
+
 			pos_elements = repmat([1:md.mesh.numberofelements]',3,1);