Index: /issm/trunk-jpl/src/m/mesh/bamg.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 16323)
+++ /issm/trunk-jpl/src/m/mesh/bamg.m	(revision 16324)
@@ -316,28 +316,48 @@
 [bamgmesh_out bamggeom_out]=BamgMesher(bamg_mesh,bamg_geometry,bamg_options);
 
-% plug results onto model
+if getfieldvalue(options,'vertical',0),
+	md.mesh=mesh2dvertical;
+	md.mesh.x=bamgmesh_out.Vertices(:,1);
+	md.mesh.z=bamgmesh_out.Vertices(:,2);
+	md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
+	md.mesh.edges=bamgmesh_out.IssmEdges;
+	md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
+	md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
+
+	%Fill in rest of fields:
+	md.mesh.numberofelements=size(md.mesh.elements,1);
+	md.mesh.numberofvertices=length(md.mesh.x);
+	md.mesh.numberofedges=size(md.mesh.edges,1);
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+	%Now, build the connectivity tables for this mesh.
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+else
+	% plug results onto model
+	md.mesh.x=bamgmesh_out.Vertices(:,1);
+	md.mesh.y=bamgmesh_out.Vertices(:,2);
+	md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
+	md.mesh.edges=bamgmesh_out.IssmEdges;
+	md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
+	md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
+
+	%Fill in rest of fields:
+	md.mesh.dimension=2;
+	md.mesh.numberofelements=size(md.mesh.elements,1);
+	md.mesh.numberofvertices=length(md.mesh.x);
+	md.mesh.numberofedges=size(md.mesh.edges,1);
+	md.mesh.z=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
+	md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
+	md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
+	md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+end
 md.private.bamg=struct();
 md.private.bamg.mesh=bamgmesh(bamgmesh_out);
 md.private.bamg.geometry=bamggeom(bamggeom_out);
-md.mesh.x=bamgmesh_out.Vertices(:,1);
-md.mesh.y=bamgmesh_out.Vertices(:,2);
-md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
-md.mesh.edges=bamgmesh_out.IssmEdges;
-md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
-md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
-
-%Fill in rest of fields:
-md.mesh.dimension=2;
-md.mesh.numberofelements=size(md.mesh.elements,1);
-md.mesh.numberofvertices=length(md.mesh.x);
-md.mesh.numberofedges=size(md.mesh.edges,1);
-md.mesh.z=zeros(md.mesh.numberofvertices,1);
-md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
-md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
-md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
-md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);
-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
 md.mesh.elementconnectivity=md.private.bamg.mesh.ElementConnectivity;
 md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
+
 
 %Check for orphan
Index: /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 16324)
+++ /issm/trunk-jpl/src/m/mesh/triangle2dvertical.m	(revision 16324)
@@ -0,0 +1,62 @@
+function md=triangle(md,domainname,resolution)
+%TRIANGLE - create model mesh using the triangle package
+%
+%   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
+%   where md is a @model object, domainname is the name of an Argus domain outline file, 
+%   and resolution is a characteristic length for the mesh (same unit as the domain outline
+%   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
+%
+%   Usage:
+%      md=triangle(md,domainname,resolution)
+%
+%   Examples:
+%      md=triangle(md,'DomainOutline.exp',1000);
+
+%Check that mesh was not already run, and warn user: 
+if md.mesh.numberofelements~=0,
+	choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s');
+	if ~strcmp(choice,'y')
+		disp('no meshing done ... exiting');
+		return
+	end
+end
+
+area=resolution^2;
+
+%Mesh using TriMesh
+[elements,x,z,segments,segmentmarkers]=TriMesh(domainname,'',area);
+
+%check that all the created nodes belong to at least one element
+orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
+for i=1:length(orphan),
+	disp('WARNING: removing orphans');
+	%get rid of the orphan node i
+	%update x and y
+	x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)];
+	z=[z(1:orphan(i)-(i-1)-1); z(orphan(i)-(i-1)+1:end)];
+	%update elements
+	pos=find(elements>orphan(i)-(i-1));
+	elements(pos)=elements(pos)-1;
+	%update segments
+	pos1=find(segments(:,1)>orphan(i)-(i-1));
+	pos2=find(segments(:,2)>orphan(i)-(i-1));
+	segments(pos1,1)=segments(pos1,1)-1;
+	segments(pos2,2)=segments(pos2,2)-1;
+end
+
+%plug into md
+md.mesh=mesh2dvertical;
+md.mesh.x=x;
+md.mesh.z=z;
+md.mesh.elements=elements;
+md.mesh.segments=segments;
+md.mesh.segmentmarkers=segmentmarkers;
+
+%Fill in rest of fields:
+md.mesh.numberofelements=length(md.mesh.elements);
+md.mesh.numberofvertices=length(md.mesh.x);
+md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+%Now, build the connectivity tables for this mesh.
+md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);