Index: /issm/trunk-jpl/src/m/mesh/augment2dmesh.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/augment2dmesh.m	(revision 19954)
+++ /issm/trunk-jpl/src/m/mesh/augment2dmesh.m	(revision 19954)
@@ -0,0 +1,75 @@
+function mh=augment2dmesh(mh,mhband,varargin)
+%AUGMENT2DMESH - augment mh mesh with a band around it (provided by mhband)
+%
+%   Usage:
+%      mh=augment2dmesh(mh,mhband);
+%
+%   Example: 
+%      md.mesh=augment2dmesh(md.mesh,md2.mesh);
+%
+
+%First process options
+options=pairoptions(varargin{:});
+
+%Offset the mesh2 elements: 
+mhband.elements=mhband.elements+mh.numberofvertices;
+mhband.segments(:,1:2)=mhband.segments(:,1:2)+mh.numberofvertices;
+mhband.segments(:,3)=mhband.segments(:,3)+mh.numberofelements;
+
+%The innner segment of md2 and the outer segments of md1 are identical. Go into  the elements of 
+%md2 and set them to their md1 equivalent: 
+flag=0;
+if flag,
+	for i=1:length(mhband.segments),
+		node2=mhband.segments(i,1);
+		%this node2 has an equivalent on the segments  of md1: 
+		for j=1:length(mh.segments),
+			node1=mh.segments(j,1);
+			if mhband.x(node2-mh.numberofvertices) == mh.x(node1) &&  mhband.y(node2-mh.numberofvertices) == mh.y(node1),
+				%go into the mesh of md2, and replace by node1.
+				pos=find(mhband.elements==node2); mhband.elements(pos)=node1;
+				segs=mhband.segments(:,1:2); pos=find(segs==node2); segs(pos)=node1; mhband.segments(:,1:2)=segs;
+				break;
+			end
+		end
+	end
+end	
+
+%Do the merge: 
+mh.elements=[mh.elements;mhband.elements];
+mh.x=[mh.x;mhband.x];
+mh.y=[mh.y;mhband.y];
+mh.lat=[mh.lat;mhband.lat];
+mh.long=[mh.long;mhband.long];
+mh.segments=[mh.segments;mhband.segments];
+
+%Remove orphans:
+x=mh.x; y=mh.y; lat=mh.lat; long=mh.long; 
+elements=mh.elements; segments=mh.segments;
+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)];
+	y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)];
+	lat=[lat(1:orphan(i)-(i-1)-1); lat(orphan(i)-(i-1)+1:end)];
+	long=[long(1:orphan(i)-(i-1)-1); long(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
+
+mh.elements=elements;
+mh.x=x;
+mh.y=y;
+mh.lat=lat;
+mh.long=long;
+mh.segments=segments;
+mh.numberofelements=length(mh.elements);
+mh.numberofvertices=length(mh.x);
Index: /issm/trunk-jpl/src/m/mesh/mesh3dsurfaceplug2d.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/mesh3dsurfaceplug2d.m	(revision 19954)
+++ /issm/trunk-jpl/src/m/mesh/mesh3dsurfaceplug2d.m	(revision 19954)
@@ -0,0 +1,77 @@
+function mh=mh3dsurfaceplug2d(mh,mh2,flags,segments,xsegs,ysegs,varargin)
+%MESH3DSURFACEPLUG2D - plug 2d mesh into a 3D surface mesh
+%
+%   Usage:
+%      mh=mesh3dsurfaceplug2d(mh,mh2);
+%
+
+	%First process options
+	options=pairoptions(varargin{:});
+
+	%Remove the elements that are flagged: 
+	pos=find(flags); 
+	mh.elements(pos,:)=[];
+	mh.numberofelements=size(mh.elements,1);
+
+	%Offset mh2.elements by number of vertices in the 3D structure: 
+	mh2.elements=mh2.elements+mh.numberofvertices;
+	mh2.segments(:,1:2)=mh2.segments(:,1:2)+mh.numberofvertices;
+	mh2.segments(:,3)=mh2.segments(:,3)+mh.numberofelements;
+
+	%The segments of md2 and the outer segments of md are identical. Go into  the elements of 
+	%dmd2 and set them to their md equivalent: 
+	for i=1:length(mh2.segments),
+		node2=mh2.segments(i,1);
+		%this node2 has an equivalent on the segments  of md: 
+		for j=1:length(segments),
+			node1=segments(j,1);
+			if mh2.x(node2-mh.numberofvertices) == xsegs(j) &&  mh2.y(node2-mh.numberofvertices) == ysegs(j),
+				%go into the mesh of md2, and replace by node1.
+				pos=find(mh2.elements==node2); mh2.elements(pos)=node1;
+				segs=mh2.segments(:,1:2); pos=find(segs==node2); segs(pos)=node1; mh2.segments(:,1:2)=segs;
+				break;
+			end
+		end
+	end
+
+	%Do the merge: 
+	mh.elements=[mh.elements;mh2.elements];
+	mh.lat=[mh.lat;mh2.lat];
+	mh.long=[mh.long;mh2.long];
+	mh.segments=[mh.segments;mh2.segments];
+	mh.numberofvertices=length(mh.lat);
+	mh.numberofelements=size(mh.elements,1);
+	
+	%Remove orphans:
+	lat=mh.lat; long=mh.long; 
+	elements=mh.elements; segments=mh.segments;
+	orphan=find(~ismember([1:length(lat)],sort(unique(elements(:)))));
+	for i=1:length(orphan),
+		%disp('WARNING: removing orphans');
+		%get rid of the orphan node i
+		%update lat and long
+		lat=[lat(1:orphan(i)-(i-1)-1); lat(orphan(i)-(i-1)+1:end)];
+		long=[long(1:orphan(i)-(i-1)-1); long(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
+	
+	mh.elements=elements;
+	mh.lat=lat;
+	mh.long=long;
+	mh.segments=segments;
+	mh.numberofelements=length(mh.elements);
+	mh.numberofvertices=length(mh.lat);
+
+	%reconstruct x,y and z:
+	R=mh.r(1);
+	mh.x = R .* cosd(mh.lat) .* cosd(mh.long);
+	mh.y = R .* cosd(mh.lat) .* sind(mh.long);
+	mh.z = R .* sind(mh.lat);
+	mh.r=sqrt(mh.x.^2+mh.y.^2+mh.z.^2);
Index: /issm/trunk-jpl/src/m/mesh/patchglobe.m
===================================================================
--- /issm/trunk-jpl/src/m/mesh/patchglobe.m	(revision 19954)
+++ /issm/trunk-jpl/src/m/mesh/patchglobe.m	(revision 19954)
@@ -0,0 +1,82 @@
+function mh=patchglobe(mh,mh2d,varargin)
+
+	%process options: 
+	options=pairoptions(varargin{:});
+
+	%recover basic options:
+	hem=getfieldvalue(options,'hem');
+	bandwidth=getfieldvalue(options,'bandwidth',100000);
+
+	%give ourselves a unique temporary directory: 
+	temproot=tempname; mkdir(temproot);
+
+	%figure out domain outline: 
+	meshtodomain(mh2d,[temproot '/Patch.exp']);
+
+	%broaden this domain outline: 
+	expcoarsen([temproot '/PatchBroad.exp'],[temproot '/Patch.exp'],200000);
+	expcontract([temproot '/PatchBroad.exp'],[temproot '/PatchBroad.exp'],-bandwidth);
+
+	%now flag vertices that are within the contour:
+	[x,y]=ll2xy(mh.lat,mh.long,hem);
+	if hem==1,
+		flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat>=0;
+	else
+		flagsnods=ContourToNodes(x,y,[temproot '/PatchBroad.exp'],1) & mh.lat<=0;
+	end
+
+	%expand flags to any element that touches the contour: 
+	pos=find(sum(flagsnods(mh.elements),2));
+	flags=zeros(mh.numberofelements,1); flags(pos)=1;
+
+	%need to find the segment enveloppe of these elements:
+	mh.segments=contourenvelope(mh,flags);
+
+	%segments need to be ordered in line: 
+	mh.segments=alignsegments(mh.segments);
+
+	%x,y for segments: 
+	[xsegs,ysegs]=ll2xy(mh.lat(mh.segments(:,1)),mh.long(mh.segments(:,1)),hem);
+
+	%create contour out of these segments:
+	meshtodomain(mh,[temproot '/PatchEnveloppe.exp'],'latlong','on');
+	expll2xy([temproot '/PatchEnveloppe.exp'],hem);
+
+	%now, create domain outine from broad enveloppe and initial mesh 
+	env=expread([temproot '/PatchEnveloppe.exp']); 
+	dom=expread([temproot '/Patch.exp']);
+	
+	%close the contours:
+	env(1).x=[env(1).x;env(1).x(1)];
+	env(1).y=[env(1).y;env(1).y(1)];
+	dom(1).x=[dom(1).x;dom(1).x(1)];
+	dom(1).y=[dom(1).y;dom(1).y(1)];
+
+	%flip inner hole: 
+	dom(1).x=flipud(dom(1).x);
+	dom(1).y=flipud(dom(1).y);
+
+	domain(1)=env; 
+	domain(2)=dom;
+	expwrite(domain,[temproot '/PatchBand.exp']);
+
+	%mesh: 
+	mdb=bamg(model(),'domain',[temproot '/PatchBand.exp'],'MaxCornerAngle',1e-15,'gradation',10000); 
+	mhb=mdb.mesh; clear mdb;
+
+	%double check: 
+	if length(mhb.segments) ~= (length(dom(1).x)+length(env(1).x)-2),
+		%error('band mesh not consistent');
+	end
+
+	%augment patch with band
+	[mhb.lat,mhb.long]=xy2ll(mhb.x,mhb.y,hem);
+	mh2db=augment2dmesh(mh2d,mhb);
+
+	%merge inner band and earth: 
+	mh=mesh3dsurfaceplug2d(mh,mh2db,flags,mh.segments,xsegs,ysegs);
+
+	%erase temporary directory: 
+	system(['rm -rf ' temproot]);
+
+end