Index: /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m	(revision 17718)
+++ /issm/trunk-jpl/src/m/contrib/dassflow/importgmsh.m	(revision 17719)
@@ -1,3 +1,3 @@
-function md = importgmsh(filename)
+function md = importgmsh(filename,dim)
 
 %some checks
@@ -30,4 +30,6 @@
 x = A(2,:)';
 y = A(3,:)';
+z = A(4,:)';
+
 A=fscanf(fid,'%s',1);
 if ~strcmp(A,'$EndNodes'), 
@@ -43,26 +45,60 @@
 disp(['Number of elements: ' num2str(nbt) ]);
 counter = 0;
-index   = zeros(0,3);
-segments       = zeros(0,2);
-segmentmarkers = zeros(0,1);
+if (dim==2),
+	index   = zeros(0,3);
+	segments       = zeros(0,2);
+	segmentmarkers = zeros(0,1);
+elseif (dim==3),
+	index   = zeros(0,4);
+	segments       = zeros(0,3);
+	segmentmarkers = zeros(0,1);
+else
+	error('not supported');
+end
+
 while(counter<nbt);
 	id = fscanf(fid,'%i',1);
 	ty = fscanf(fid,'%i',1);
-	if(ty>2) error('Type not supported'); end
 	nbf = fscanf(fid,'%i',1);
 	flags = fscanf(fid,'%i',nbf);
 
-	if ty==1,
-		A=fscanf(fid,'%i %i',2);
-		segments(end+1,:)=A;
-		if(flags(1)==5 & flags(2)==3), segmentmarkers(end+1)=3; 
-		elseif(flags(1)==1 & flags(2)==4), segmentmarkers(end+1)=4;
-		elseif(flags(1)==2 & flags(2)==1), segmentmarkers(end+1)=1;
-		elseif(flags(1)==4 & flags(2)==2), segmentmarkers(end+1)=2;
-		else error(['flags ' num2str(flags') ' not supported']);
-		end
-	else
-		A=fscanf(fid,'%i %i %i',3);
-		index(end+1,:)=A;
+	switch(ty)
+		case 1, %segments
+			A=fscanf(fid,'%i %i',2);
+			if (dim==2),  %Actual element
+				segments(end+1,:)=A;
+				if    (flags(1)==5 & flags(2)==3), segmentmarkers(end+1)=3; 
+				elseif(flags(1)==1 & flags(2)==4), segmentmarkers(end+1)=4;
+				elseif(flags(1)==2 & flags(2)==1), segmentmarkers(end+1)=1;
+				elseif(flags(1)==4 & flags(2)==2), segmentmarkers(end+1)=2;
+				else error(['flags ' num2str(flags') ' not supported']);
+				end
+			else
+				error('not supported');
+			end
+		case 2, %tria
+			A=fscanf(fid,'%i %i %i',3);
+			if (dim==2), %Actual element
+				index(end+1,:)=A;
+			else         %Boundary element
+				segments(end+1,:)=A;
+				if    (flags(1)==3), segmentmarkers(end+1)=3; 
+				elseif(flags(1)==4), segmentmarkers(end+1)=4;
+				elseif(flags(1)==1), segmentmarkers(end+1)=1;
+				elseif(flags(1)==2), segmentmarkers(end+1)=2;
+				elseif(flags(1)==5), segmentmarkers(end+1)=5;
+				elseif(flags(1)==6), segmentmarkers(end+1)=6;
+				else error(['flags ' num2str(flags') ' not supported']);
+				end
+			end
+		case 4, %tetra
+			A=fscanf(fid,'%i %i %i %i',4);
+			if (dim==3), %Actual element
+				index(end+1,:)=A;
+			else
+				error('not supported');
+			end
+		otherwise,
+			error(['Type ' num2str(ty) ' not supported']);
 	end
 	counter = counter + 1;
@@ -70,9 +106,18 @@
 
 %recreate segments
-nbs = size(segments,1);
-segments = [segments zeros(nbs,1)];
-for i=1:nbs,
-	E = find(sum(ismember(index,segments(i,:)),2)>1);
-	segments(i,3)=E;
+if dim==2,
+	nbs = size(segments,1);
+	segments = [segments zeros(nbs,1)];
+	for i=1:nbs,
+		E = find(sum(ismember(index,segments(i,:)),2)>1);
+		segments(i,3)=E;
+	end
+else
+	nbs = size(segments,1);
+	segments = [segments zeros(nbs,1)];
+	for i=1:nbs,
+		E = find(sum(ismember(index,segments(i,:)),2)>2);
+		segments(i,4)=E;
+	end
 end
 
@@ -81,10 +126,25 @@
 
 %Create model
-md=meshconvert(model,index,x,y);
-md.mesh=mesh2dvertical(md.mesh);
-md.mesh.segmentmarkers=segmentmarkers;
-md.mesh.segments=segments;
-md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonbase(find(vertexflags(md.mesh,1)))=1;
-md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface(find(vertexflags(md.mesh,3)))=1;
+if 0, %2d triangles
+	md=meshconvert(model,index,x,y);
+	md.mesh=mesh2dvertical(md.mesh);
+	md.mesh.segmentmarkers=segmentmarkers;
+	md.mesh.segments=segments;
+	md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonbase(find(vertexflags(md.mesh,1)))=1;
+	md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonsurface(find(vertexflags(md.mesh,3)))=1;
+else
+	md=model();
+	md.mesh=mesh3dprisms();
+	md.mesh.x = x;
+	md.mesh.y = y;
+	md.mesh.z = z;
+	md.mesh.elements = index;
+	md.mesh.numberofelements=size(md.mesh.elements,1);
+	md.mesh.numberofvertices=length(md.mesh.x);
+	md.mesh.vertexonbase=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonbase(find(segments(find(segmentmarkers==1),1:3)))=1;
+	md.mesh.vertexonsurface=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonsurface(find(segments(find(segmentmarkers==3),1:3)))=1;
+end