Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshband.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshband.m	(revision 21368)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/meshband.m	(revision 21368)
@@ -0,0 +1,58 @@
+function band=meshband(mh,outerdomain,varargin)
+
+	%process options: 
+	options=pairoptions(varargin{:});
+
+	%some checks on the mesh: 
+	if (isempty(mh.x) | isempty(mh.y)  ) 
+		error('meshband error message: mesh has one of the following empty: ''x'',''y''');
+	end
+
+	%give ourselves a unique temporary directory: 
+	temproot=tempname; mkdir(temproot);
+
+	%align external segments on internal mesh: 
+	mh.segments=alignsegments(mh.segments);
+
+	%figure out domain outline: 
+	meshtodomain(mh,[temproot '/innerdomain.exp']);
+
+	%create domain outine from inner and outer domain
+	inner=expread([temproot '/innerdomain.exp']);
+	if inner.closed==0,
+		inner.nods=inner.nods+1;
+		inner.x(end+1)=inner.x(1);
+		inner.y(end+1)=inner.y(1);
+		inner.closed=1;
+		if getfieldvalue(options,'invert',0),
+			inner.x=flipud(inner.x);
+			inner.y=flipud(inner.y);
+		end
+	end
+	outer=expread(outerdomain);
+
+	domain(1)=outer;
+	domain(2)=inner;
+	
+	expwrite(domain,[temproot '/Band.exp']);
+	
+	if getfieldvalue(options,'plot',0),
+		figure(1),clf,hold on,axis image;
+		expdisp([temproot '/Band.exp'],'linestyle','k-*');
+	end
+
+	%mesh: 
+	md=bamg(model(),'domain',[temproot '/Band.exp'],'MaxCornerAngle',1e-15,'gradation',10000); band=md.mesh; clear md;
+
+	if getfieldvalue(options,'plot',0),
+		trisurf(band.elements,band.x,band.y,band.x),view(2),shading faceted;
+	end
+	
+	if length(band.segments) ~= (length(domain(1).x)+length(domain(2).x)-2),
+		error('band mesh not consistent');
+	end
+
+	%erase temporary directory: 
+	system(['rm -rf ' temproot]);
+
+end
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m	(revision 21368)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/mesh/modelmerge.m	(revision 21368)
@@ -0,0 +1,89 @@
+function md=modelmerge(md1,md2)
+%MODELMERGE  - merge two models by merging their meshes
+%
+%   Usage:
+%      md=modelmerge(md1,md2);
+
+	if nargin~=2 
+		help meshconvert
+		error('meshconvert error message: bad usage');
+	end
+	
+	md=model();
+
+	%first ,copy md1 mesh into md.mesh to initialize: 
+	md.mesh=md1.mesh;
+
+	%some initializatoin: 
+	elements1=md1.mesh.elements;
+	x1=md1.mesh.x;
+	y1=md1.mesh.y;
+	nods1=md1.mesh.numberofvertices;
+	nel1=md1.mesh.numberofelements;
+
+	elements2=md2.mesh.elements;
+	x2=md2.mesh.x;
+	y2=md2.mesh.y;
+	nods2=md2.mesh.numberofvertices;
+	nel2=md2.mesh.numberofelements;
+	segs2=md2.mesh.segments;
+
+	%offset elements2 by nods1: 
+	elements2=elements2+nods1;
+
+	tolerance=1e-5;
+	%go into the vertices on boundary of mesh 1, and figure out which ones are common to mesh2: 
+	verticesonboundary=find(md1.mesh.vertexonboundary); 
+	for i=1:length(verticesonboundary),
+		node1=verticesonboundary(i); xnode1=x1(node1); ynode1=y1(node1);
+		%is there another node with these coordinates in mesh2? 
+		ind=find(sqrt(((x2-xnode1).^2+(y2-ynode1).^2))<tolerance);
+		if ~isempty(ind),
+			x2(ind)=NaN;
+			y2(ind)=NaN;
+			pos=find(elements2==(ind+nods1)); elements2(pos)=node1;
+		end
+	end
+
+	%go through elements2 and drop counter on each vertex that is above the x2 and y2 vertices being dropped: 
+	while( ~isempty(find(isnan(x2)))),
+		for i=1:length(x2),
+			if isnan(x2(i)),
+				pos=find(elements2>(i+nods1));
+				elements2(pos)=elements2(pos)-1;
+				x2(i)=[];
+				y2(i)=[];
+				break;
+			end
+		end
+	end
+
+	%merge elements: 
+	elements=[elements1;elements2];
+
+	%merge vertices: 
+	x=[x1;x2]; 
+	y=[y1;y2];
+
+	%output: 
+	md.mesh.x=x;
+	md.mesh.y=y;
+	md.mesh.elements=elements;
+	md.mesh.numberofvertices=length(x);
+	md.mesh.numberofelements=size(elements,1);
+
+	%connectivities: 
+	md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+	md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+
+	%find segments: 
+	md.mesh.segments=findsegments(md);
+
+	%vertex on boundary: 
+	md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1);
+	md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+	%some checks: 
+	if max(md.mesh.elements)>md.mesh.numberofvertices, 
+		error('issue in modelmerge, one of the element ids is > number of vertices!');
+	end