Index: /issm/trunk/src/m/classes/public/mesh/meshnodensity.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshnodensity.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/meshnodensity.m	(revision 2737)
@@ -0,0 +1,82 @@
+function md=meshnodensity(md,domainname,varargin)
+%MESH - create model mesh
+%
+%   This routine creates a model mesh using TriMeshNoDensity and a domain outline
+%   where md is a @model object, domainname is the name of an Argus domain outline file, 
+%   Riftname is an optional argument (Argus domain outline) describing rifts.
+%   The  difference with mesh.m is that the resolution of the mesh follows that of the domain 
+%   outline and the riftoutline
+%
+%   Usage:
+%      md=mesh(md,domainname)
+%   or md=mesh(md,domainname,riftname);
+%
+%   Examples:
+%      md=mesh(md,'DomainOutline.exp');
+%      md=mesh(md,'DomainOutline.exp','Rifts.exp');
+
+if (nargin==2),
+	riftname='';
+end
+if (nargin==3),
+	riftname=varargin{1};
+end
+
+%Mesh using TriMeshNoDensity
+if strcmp(riftname,''),
+	[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMeshNoDensity(domainname);
+else
+	[elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname);
+
+	%check that all the created grids belong to at least one element
+	orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
+	for i=1:length(orphan),
+		%get rid of the orphan grid 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)];
+		%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.x=x;
+	md.y=y;
+	md.elements=elements;
+	md.segments=segments;
+	md.segmentmarkers=segmentmarkers;
+end
+
+%Fill in rest of fields:
+md.numberofelements=length(md.elements);
+md.numberofgrids=length(md.x);
+md.z=zeros(md.numberofgrids,1);
+md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
+md.gridonbed=ones(md.numberofgrids,1);
+md.gridonsurface=ones(md.numberofgrids,1);
+md.elementonbed=ones(md.numberofelements,1);
+md.elementonsurface=ones(md.numberofelements,1);
+
+%Now, build the connectivity tables for this mesh.
+md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+
+%outline names
+md.domainoutline=char(textread(domainname,'%s','delimiter','\n'));
+if strcmp(riftname,''),
+	md.riftoutline='';
+else
+	md.riftoutline=readfile(riftname);
+end
+
+%type of model
+md.type='2d';
+	
+%augment counter  keeping track of what has been done to this model
+md.counter=1;
Index: sm/trunk/src/m/classes/public/mesh/meshprocessrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshprocessrifts.m	(revision 2736)
+++ 	(revision )
@@ -1,28 +1,0 @@
-function md=meshprocessrifts(md)
-%MESHPROCESSRIFTS - process mesh when rifts are present
-%
-%   split rifts inside mesh (rifts are defined by presence of
-%   segments inside the domain outline)
-%
-%   Usage:
-%      md=meshprocessrifts(md)
-
-%some checks on list of arguments
-if ((nargin~=1) | (nargout~=1)),
-	help meshprocessrifts
-	error('meshprocessrifts error messagei: bad usage');
-end
-
-[md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);
-
-%Fill in rest of fields:
-md.numberofelements=length(md.elements);
-md.numberofgrids=length(md.x);
-md.z=zeros(md.numberofgrids,1);
-md.gridonboundary=zeros(length(md.x),1); md.gridonboundary(md.segments(:,1:2))=1;
-md.numrifts=length(md.rifts);
-md.elements_type=3*ones(md.numberofelements,1);
-md.gridonbed=ones(md.numberofgrids,1);
-md.gridonsurface=ones(md.numberofgrids,1);
-md.elementonbed=ones(md.numberofelements,1);
-md.elementonsurface=ones(md.numberofelements,1);
Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 2736)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 2737)
@@ -54,4 +54,5 @@
 else
 	md=mesh(md,domainoutline,riftoutline,resolution);
+	md=meshprocessrifts(md);
 end
 disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
@@ -90,4 +91,13 @@
 	md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
 
+	%if we have rifts, we just messed them up, we need to recreate the segments that constitute those 
+	%rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
+	if md.numrifts, 
+		md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+		md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+		md.segments=findsegments(md);
+		md=meshyamsrecreateriftsegments(md);
+	end
+
 end
 	
@@ -123,2 +133,44 @@
 end
 md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
+
+%deal with rifts 
+if md.numrifts,
+	%first, recreate rift segments
+	md=meshyamsrecreateriftsegments(md);
+
+	%using the segments, recreate the penaltypairs
+	for j=1:md.numrifts,
+		rift=md.rifts(j);
+
+		%build normals and lengths of segments:
+		lengths=sqrt((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 );
+		normalsx=cos(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
+		normalsy=sin(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
+
+		%ok, build penaltypairs: 
+		numpenaltypairs=length(rift.segments)/2-1;
+		rift.penaltypairs=zeros(numpenaltypairs,7);
+
+		for i=1:numpenaltypairs,
+			rift.penaltypairs(i,1)=rift.segments(i,2);
+			rift.penaltypairs(i,2)=rift.segments(end-i,2);
+			rift.penaltypairs(i,3)=rift.segments(i,3);
+			rift.penaltypairs(i,4)=rift.segments(end-i,3);
+			rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
+			rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
+			rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
+		end
+		%renormalize norms: 
+		norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
+		rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
+		rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
+
+		md.rifts(j)=rift;
+	end
+
+end
+
+
+
+
+
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshaddrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshaddrifts.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshaddrifts.m	(revision 2737)
@@ -0,0 +1,98 @@
+function md=meshaddrifts(md,riftname);
+%MESHADDRIFTS - add rifts to a preloaded mesh (typically, an argus mesh)
+%
+%   Usage:
+%      md=meshaddrifts(md,riftname);
+%
+%        where md is a model with a preexisting mesh, and riftname is the name of an .exp file.
+%        The format of the riftname file is as follows: a list of pairs of closed and open contours. 
+%        The closed contour defines the envelopped of the rift, the open contour that follows in the 
+%        file defines the rift. The density of the rift should be chosen carefully in the file, as it 
+%        will be used to defined the rift contour density of the mesh. The open contour density will 
+%        be preserved. There can be as many pairs of closed contour and rift contour as wished.
+
+
+
+
+%read rift: 
+domains=expread(riftname,1);
+contours=domains(1:2:end);
+rifts=domains(2:2:end);
+
+%now loop over rifts: 
+for rift_i=1:length(rifts),
+	
+	%refine rift to desired resolution: 
+	contour=contours(rift_i);
+	rift=rifts(rift_i);
+	
+	delete('Meshaddrifts.Rift.exp');
+	expwrite(rift,'Meshaddrifts.Rift.Coarse.exp');
+	expcoarsen('Meshaddrifts.Rift.exp','Meshaddrifts.Rift.Coarse.exp',rift.density);
+	delete('Meshaddrifts.Rift.Coarse.exp');
+	
+	%extract model:
+	expwrite(contour,'Meshaddrifts.Contour.exp');
+	md2=modelextract(md,'Meshaddrifts.Contour.exp');
+	
+	%create domain of md2 model: 
+	md2.segments=contourenvelope(md2,'Meshaddrifts.Contour.exp');
+	domain_index=md2.segments(1,1:2);
+	while (domain_index(end)~=domain_index(1)),
+		pos=find(md2.segments(:,1)==domain_index(end));
+		domain_index(end+1)=md2.segments(pos,2);
+	end
+	
+	domain.x=md2.x(domain_index);
+	domain.y=md2.y(domain_index);
+	domain.name='Meshaddrifts.Domain.exp';
+	domain.density=1;
+	expwrite(domain,'Meshaddrifts.Domain.exp');
+	
+	%unloop domain index: used for later.
+	domain_index=domain_index(1:end-1);
+	
+	%remesh md2 using new domain outline, and rift profile: 
+	md2=meshnodensity(md2,'Meshaddrifts.Domain.exp','Meshaddrifts.Rift.exp');
+	md2=meshprocessrifts(md2);
+	
+	%plug md2 mesh into md mesh: 
+	[md.elements,md.x,md.y,md.z,md.numberofelements,md.numberofgrids,elconv,nodeconv,elconv2,nodeconv2]=meshplug(md.elements,md.x,md.y,md.z,...
+								md2.elements,md2.x,md2.y,md2.z,md2.extractedgrids,md2.extractedelements,domain_index);
+
+	%update md2 rifts using elconv and nodeconv, and plug them into md: 
+	md2.rifts=updateriftindexing(md2.rifts,elconv2,nodeconv2);
+	
+	for i=1:md.numrifts,
+		md.rifts(i)=updateriftindexing(md.rifts(i),elconv,nodeconv);
+	end
+	
+	if md.numrifts==0,
+		md.rifts=md2.rifts;
+		md.numrifts=1;
+	else
+		md.rifts(end+1,1)=md2.rifts;
+		md.numrifts=md.numrifts+1;
+	end
+	
+	md.segments(:,1:2)=nodeconv(md.segments(:,1:2));
+	md.segments(:,3)=elconv(md.segments(:,3));
+
+end
+
+%finish up "a la" mesh.h
+md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
+md.gridonbed=ones(md.numberofgrids,1);
+md.gridonsurface=ones(md.numberofgrids,1);
+md.elementonbed=ones(md.numberofelements,1);
+md.elementonsurface=ones(md.numberofelements,1);
+
+%Now, build the connectivity tables for this mesh.
+md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+
+%type of model
+md.type='2d';
+	
+%augment counter  keeping track of what has been done to this model
+md.counter=1;
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshplug.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshplug.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshplug.m	(revision 2737)
@@ -0,0 +1,89 @@
+function [elements,x,y,z,numberofelements,numberofgrids,elconv,nodeconv,elconv2,nodeconv2]=meshplug(elements,x,y,z,elements2,x2,y2,z2,extractedgrids,extractedelements,domain);
+%MESHPLUG - embed mesh into another one
+%     See also meshaddrifts
+
+
+%initialize elconv,nodeconv conversion tables from md mesh to new md mesh
+elconv=1:size(elements,1); elconv=elconv';
+nodeconv=1:size(x,1); nodeconv=nodeconv';
+
+%take away old elements in area of interest: 
+elements(extractedelements,:)=[];
+element_offset=size(elements,1);
+
+%update elconv after having extracted the area of interest elements
+temp_elconv=elconv; temp_elconv(extractedelements)=[];
+temp_elconvnum=1:length(temp_elconv);
+elconv(temp_elconv)=temp_elconvnum;
+elconv(extractedelements)=NaN;
+
+%initialize elconv2 and nodeconv2, conversion tables from md2 mesh to new md mesh
+elconv2=1:size(elements2,1);elconv2=elconv2'+element_offset;
+nodeconv2=(size(x,1)+1):(size(x,1)+size(x2,1)); nodeconv2=nodeconv2';
+
+extractedgrids_minusborder=extractedgrids;
+extractedgrids_minusborder(domain)=[];
+
+x(extractedgrids_minusborder)=NaN;
+y(extractedgrids_minusborder)=NaN;
+
+%now, plug md2 mesh: 
+
+%first, offset all ids of md2 mesh
+elements2=elements2+length(x);
+
+%NaN border grids in second mesh
+x2(1:length(domain))=NaN;
+y2(1:length(domain))=NaN;
+
+%redirect border grids in elements2  to elements
+for i=1:length(domain),
+	pos=find(elements2==(i+length(x)));
+	elements2(pos)=extractedgrids(domain(i));
+end
+
+%same deal for nodeconv2:
+for i=1:length(domain),
+	nodeconv2(i)=extractedgrids(domain(i));
+end
+
+
+%plug elements
+elements=[elements;elements2];
+
+
+%now, increase number of grids
+x=[x; x2];
+y=[y; y2];
+z=[z; z2];
+
+%now, get rid of NaN in x:
+while  ~isempty(find(isnan(x))),
+
+	pos=find(isnan(x));
+	grid=pos(1);
+
+	%collapse grid
+	x(grid)=[];
+	y(grid)=[];
+	z(grid)=[];
+
+	%renumber all grids > grid in elements
+	pos=find(elements>grid);
+	elements(pos)=elements(pos)-1;
+
+	%same deal for nodeconv2: 
+	pos=find(nodeconv2>grid);
+	nodeconv2(pos)=nodeconv2(pos)-1;
+
+end
+
+numberofgrids=length(x);
+numberofelements=length(elements);
+
+%finish nodeconv: 
+temp_nodeconv=nodeconv;  temp_nodeconv(extractedgrids_minusborder)=[];
+temp_nodeconvnum=1:length(temp_nodeconv);
+nodeconv(temp_nodeconv)=temp_nodeconvnum;
+nodeconv(extractedgrids_minusborder)=NaN;
+
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshprocessrifts.m	(revision 2737)
@@ -0,0 +1,34 @@
+function md=meshprocessrifts(md)
+%MESHPROCESSRIFTS - process mesh when rifts are present
+%
+%   split rifts inside mesh (rifts are defined by presence of
+%   segments inside the domain outline)
+%
+%   Usage:
+%      md=meshprocessrifts(md)
+
+%some checks on list of arguments
+if ((nargin~=1) | (nargout~=1)),
+	help meshprocessrifts
+	error('meshprocessrifts error messagei: bad usage');
+end
+
+[md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);
+
+%Fill in rest of fields:
+md.numberofelements=length(md.elements);
+md.numberofgrids=length(md.x);
+md.z=zeros(md.numberofgrids,1);
+md.gridonboundary=zeros(length(md.x),1); md.gridonboundary(md.segments(:,1:2))=1;
+md.numrifts=length(md.rifts);
+md.elements_type=3*ones(md.numberofelements,1);
+md.gridonbed=ones(md.numberofgrids,1);
+md.gridonsurface=ones(md.numberofgrids,1);
+md.elementonbed=ones(md.numberofelements,1);
+md.elementonsurface=ones(md.numberofelements,1);
+
+%get coordinates of rift tips
+for i=1:md.numrifts,
+	md.rifts(i).tip1coordinates=[md.x(md.rifts(i).tips(1)) md.y(md.rifts(i).tips(1))];
+	md.rifts(i).tip2coordinates=[md.x(md.rifts(i).tips(2)) md.y(md.rifts(i).tips(2))];
+end
Index: /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/meshyamsrecreateriftsegments.m	(revision 2737)
@@ -0,0 +1,34 @@
+function md=meshyamsrecreateriftsegments(md)
+
+	%recreate rift segments: just used for yams. temporaroy routine.
+	pos_record=[];
+	if md.numrifts,
+		for i=1:md.numrifts,
+			rift=md.rifts(i);
+
+			%find tip1 and tip2 for this rift, in the new mesh created by yams.
+			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+			tip1=md.segments(pos,1);
+			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip2coordinates(1),rift.tip2coordinates(2));
+			tip2=md.segments(pos,1);
+
+			%start from tip1, and build segments of this rift. 
+			pos=find_point(md.x(md.segments(:,1)),md.y(md.segments(:,1)),rift.tip1coordinates(1),rift.tip1coordinates(2));
+			pos_record=[pos_record; pos];
+			riftsegs=md.segments(pos,:);
+			while 1,
+				A=riftsegs(end,1); B=riftsegs(end,2); el=riftsegs(end,3);
+				%find other segment that holds B.
+				pos=find(md.segments(:,1)==B);
+				pos_record=[pos_record; pos];
+				riftsegs=[riftsegs; md.segments(pos,:)];
+				if riftsegs(end,2)==tip1, 
+					break;
+				end
+			end
+			md.rifts(i).segments=riftsegs;
+			md.rifts(i).tips=[tip1 tip2];
+		end
+	end
+	%take out rift segments from segments
+	md.segments(pos_record,:)=[];
Index: /issm/trunk/src/m/classes/public/mesh/rifts/rifttipsrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/rifttipsrefine.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/rifttipsrefine.m	(revision 2737)
@@ -0,0 +1,26 @@
+function md=rifttipsrefine(md,filename,resolution,circleradius);
+%RIFTTIPSREFINE - refine mesh near rift tips
+%
+%   Usage:
+%      md=rifttipsrefine(md,resolution,circleradius);
+
+numberofgrids=50;
+
+%take rifts, and create refinement circles around tips
+rifts=expread(filename,1);
+
+!echo -n "" > Circles.exp
+for i=1:length(rifts),
+	tip1=[rifts(i).x(1) rifts(i).y(1)];
+	tip2=[rifts(i).x(end) rifts(i).y(end)];
+	%create circle around tip
+	expcreatecircle('Circle1.exp',tip1(1),tip1(2),circleradius,numberofgrids);
+	expcreatecircle('Circle2.exp',tip2(1),tip2(2),circleradius,numberofgrids);
+	!cat Circles.exp Circle1.exp Circle2.exp > Circles2.exp
+	!mv Circles2.exp Circles.exp
+	!rm -rf Circle1.exp Circle2.exp
+end
+
+md=meshexprefine(md,'Circles.exp',resolution);
+
+system('rm -rf Circles.exp');
Index: /issm/trunk/src/m/classes/public/mesh/rifts/updateriftindexing.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifts/updateriftindexing.m	(revision 2737)
+++ /issm/trunk/src/m/classes/public/mesh/rifts/updateriftindexing.m	(revision 2737)
@@ -0,0 +1,11 @@
+function rift=updateriftindexing(rift,elconv,nodeconv)
+%UPDATERIFTINDEXING - update rift indexing, using mesh to new mesh conversion tables
+%     See also meshaddrift
+
+rift.segments(:,1:2)=nodeconv(rift.segments(:,1:2));
+rift.segments(:,3)=elconv(rift.segments(:,3));
+rift.pairs=elconv(rift.pairs);
+rift.tips=nodeconv(rift.tips);
+
+rift.penaltypairs(:,1:2)=nodeconv(rift.penaltypairs(:,1:2));
+rift.penaltypairs(:,3:4)=elconv(rift.penaltypairs(:,3:4));
Index: sm/trunk/src/m/classes/public/mesh/rifttipsrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifttipsrefine.m	(revision 2736)
+++ 	(revision )
@@ -1,33 +1,0 @@
-function md=rifttipsrefine(md,resolution,circleradius);
-%RIFTTIPSREFINE - refine mesh near rift tips
-%
-%   Usage:
-%      md=rifttipsrefine(md,resolution,circleradius);
-
-numberofgrids=50;
-
-%First check rifts are present!
-if strcmpi(md.riftoutline,''),
-	error('No rifts found in the outlines!');
-end
-
-%Ok, dumps riftoutline into a temporary file
-writefile('Rifts.exp.temp',md.riftoutline);
-
-%take rifts, and create refinement circles around tips
-rifts=expread('Rifts.exp.temp',1);
-!echo -n "" > Circles.exp
-for i=1:length(rifts),
-	tip1=[rifts(i).x(1) rifts(i).y(1)];
-	tip2=[rifts(i).x(end) rifts(i).y(end)];
-	%create circle around tip
-	expcreatecircle('Circle1.exp',tip1(1),tip1(2),circleradius,numberofgrids);
-	expcreatecircle('Circle2.exp',tip2(1),tip2(2),circleradius,numberofgrids);
-	!cat Circles.exp Circle1.exp Circle2.exp > Circles2.exp
-	!mv Circles2.exp Circles.exp
-	!rm -rf Circle1.exp Circle2.exp
-end
-
-md=meshexprefine(md,'Circles.exp',resolution);
-
-system('rm -rf Circles.exp');