Index: /issm/trunk-jpl/src/m/model/contourenvelope.m
===================================================================
--- /issm/trunk-jpl/src/m/model/contourenvelope.m	(revision 12611)
+++ /issm/trunk-jpl/src/m/model/contourenvelope.m	(revision 12612)
@@ -35,8 +35,8 @@
 %Now, build the connectivity tables for this mesh.
 %Computing connectivity
-if size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices,
+if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d),
 	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
 end
-if size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements,
+if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d),
 	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
 end
@@ -44,9 +44,24 @@
 %get nodes inside profile
 mesh.elementconnectivity=md.mesh.elementconnectivity;
+if md.mesh.dimension==2;
+	mesh.elements=md.mesh.elements;
+	mesh.x=md.mesh.x;
+	mesh.y=md.mesh.y;
+	mesh.numberofvertices=md.mesh.numberofvertices;
+	mesh.numberofelements=md.mesh.numberofelements;
+else
+	mesh.elements=md.mesh.elements2d;
+	mesh.x=md.mesh.x2d;
+	mesh.y=md.mesh.y2d;
+	mesh.numberofvertices=md.mesh.numberofvertices2d;
+	mesh.numberofelements=md.mesh.numberofelements2d;
+end
+
 if nargin==2,
+
 	if isfile,
 		%get flag list of elements and nodes inside the contour
-		nodein=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,file,'node',1);
-		elemin=(sum(nodein(md.mesh.elements),2)==size(md.mesh.elements,2));
+		nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1);
+		elemin=(sum(nodein(mesh.elements),2)==size(mesh.elements,2));
 		%modify element connectivity
 		elemout=find(~elemin);
@@ -55,10 +70,10 @@
 	else
 		%get flag list of elements and nodes inside the contour
-		nodein=zeros(md.mesh.numberofvertices,1); 
-		elemin=zeros(md.mesh.numberofelements,1); 
+		nodein=zeros(mesh.numberofvertices,1); 
+		elemin=zeros(mesh.numberofelements,1); 
 		
 		pos=find(flags); 
 		elemin(pos)=1;
-		nodein(md.mesh.elements(pos,:))=1;
+		nodein(mesh.elements(pos,:))=1;
 
 		%modify element connectivity
@@ -87,11 +102,11 @@
 	els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
 	if length(els2)>1,
-		flag=intersect(md.mesh.elements(els2(1),:),md.mesh.elements(els2(2),:));
-		nods1=md.mesh.elements(el1,:);
+		flag=intersect(mesh.elements(els2(1),:),mesh.elements(els2(2),:));
+		nods1=mesh.elements(el1,:);
 		nods1(find(nods1==flag))=[];
 		segments(count,:)=[nods1 el1];
 
-		ord1=find(nods1(1)==md.mesh.elements(el1,:));
-		ord2=find(nods1(2)==md.mesh.elements(el1,:));
+		ord1=find(nods1(1)==mesh.elements(el1,:));
+		ord2=find(nods1(2)==mesh.elements(el1,:));
 
 		%swap segment nodes if necessary
@@ -104,12 +119,12 @@
 		count=count+1;
 	else
-		nods1=md.mesh.elements(el1,:);
-		flag=setdiff(nods1,md.mesh.elements(els2,:));
+		nods1=mesh.elements(el1,:);
+		flag=setdiff(nods1,mesh.elements(els2,:));
 		for j=1:3,
 			nods=nods1; nods(j)=[];
 			if any(ismember(flag,nods)),
 				segments(count,:)=[nods el1];
-				ord1=find(nods(1)==md.mesh.elements(el1,:));
-				ord2=find(nods(2)==md.mesh.elements(el1,:));
+				ord1=find(nods(1)==mesh.elements(el1,:));
+				ord2=find(nods(2)==mesh.elements(el1,:));
 				if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
 					temp=segments(count,1);
Index: /issm/trunk-jpl/src/m/model/modelextract.m
===================================================================
--- /issm/trunk-jpl/src/m/model/modelextract.m	(revision 12611)
+++ /issm/trunk-jpl/src/m/model/modelextract.m	(revision 12612)
@@ -210,4 +210,14 @@
 		md2.mesh.segments=contourenvelope(md2);
 		md2.mesh.vertexonboundary=zeros(numberofvertices2,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+	else
+		%First do the connectivity for the contourenvelope in 2d
+		md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements2d,md2.mesh.numberofvertices2d);
+		md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity);
+		md2.mesh.segments=contourenvelope(md2);
+		md2.mesh.vertexonboundary=zeros(numberofvertices2/md2.mesh.numberoflayers,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2))=1;
+		md2.mesh.vertexonboundary=repmat(md2.mesh.vertexonboundary,md2.mesh.numberoflayers,1);
+		%Then do it for 3d as usual
+		md2.mesh.vertexconnectivity=NodeConnectivity(md2.mesh.elements,md2.mesh.numberofvertices);
+		md2.mesh.elementconnectivity=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity);
 	end