Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 26830)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 26831)
@@ -928,13 +928,31 @@
 			x_edges = 0.5*(md.mesh.x(edges(:,1)) + md.mesh.x(edges(:,2)));
 			y_edges = 0.5*(md.mesh.y(edges(:,1)) + md.mesh.y(edges(:,2)));
-			md2.mesh.x = [md2.mesh.x;x_edges];
-			md2.mesh.y = [md2.mesh.y;y_edges];
-			md2.mesh.elements = [...
+			xnew = [md2.mesh.x;x_edges];
+			ynew = [md2.mesh.y;y_edges];
+			indexnew = [...
 				index(:,1)          nbv+edges_tria(:,3) nbv+edges_tria(:,2);...
 				nbv+edges_tria(:,2) nbv+edges_tria(:,3) nbv+edges_tria(:,1);...
 				nbv+edges_tria(:,2) nbv+edges_tria(:,1) index(:,3);...
 				nbv+edges_tria(:,3) index(:,2)          nbv+edges_tria(:,1)];
-			md2.mesh.numberofelements = 4*nbe;
-			md2.mesh.numberofvertices = nbv + size(edges,1);
+			%md2.mesh.numberofelements = 4*nbe;
+			%md2.mesh.numberofvertices = nbv + size(edges,1);
+
+			%Call Bamg to update other mesh properties
+			[bamgmesh_out bamggeom_out]=BamgConvertMesh(indexnew,xnew,ynew);
+			md2.mesh.x              = bamgmesh_out.Vertices(:,1);
+			md2.mesh.y              = bamgmesh_out.Vertices(:,2);
+			md2.mesh.elements       = bamgmesh_out.Triangles(:,1:3);
+			md2.mesh.edges          = bamgmesh_out.IssmEdges;
+			md2.mesh.segments       = bamgmesh_out.IssmSegments(:,1:3);
+			md2.mesh.segmentmarkers = bamgmesh_out.IssmSegments(:,4);
+			md2.mesh.numberofelements = size(md2.mesh.elements,1);
+			md2.mesh.numberofvertices = length(md2.mesh.x);
+			md2.mesh.numberofedges    = size(md2.mesh.edges,1);
+			md2.mesh.vertexonboundary = zeros(md2.mesh.numberofvertices,1); md2.mesh.vertexonboundary(md2.mesh.segments(:,1:2)) = 1;
+
+			%Deal with boudary
+			md2.mesh.vertexonboundary = [md.mesh.vertexonboundary;sum(md.mesh.vertexonboundary(edges),2)==2];
+			md2.mesh.elementconnectivity=bamgmesh_out.ElementConnectivity;
+			md2.mesh.elementconnectivity(find(isnan(md2.mesh.elementconnectivity)))=0;
 			disp(['   Old number of elements: ' num2str(nbe)]);
 			disp(['   New number of elements: ' num2str(4*nbe)]);