Index: sm/trunk/src/m/classes/public/mesh.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh.m	(revision 1213)
+++ 	(revision )
@@ -1,95 +1,0 @@
-function md=mesh(md,domainname,varargin)
-%MESH - create model mesh
-%
-%   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=mesh(md,domainname,resolution)
-%   or md=mesh(md,domainname,riftname, resolution)
-%
-%   Examples:
-%      md=mesh(md,'DomainOutline.exp',1000);
-%      md=mesh(md,'DomainOutline.exp','Rifts.exp',1500);
-
-%Figure out a characteristic area. Resolution is a grid oriented concept (ex a 1000m  resolution grid would 
-%be made of 1000*1000 area squares). 
-if (nargin==3),
-	resolution=varargin{1};
-	riftname='';
-end
-if (nargin==4),
-	riftname=varargin{1};
-	resolution=varargin{2};
-end
-
-%Check that mesh was not already run, and warn user: 
-if subsref(md,struct('type','.','subs','counter'))>=1,
-	choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s');
-	if ~strcmp(choice,'y')
-		error('no meshing done ... exiting');
-	end
-end
-
-area=resolution^2;
-
-%Mesh using TriMesh
-if strcmp(riftname,''),
-	[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes');
-else
-	[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes');
-
-	%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=readfile(domainname);
-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: /issm/trunk/src/m/classes/public/mesh/mesh.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/mesh.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/mesh.m	(revision 1214)
@@ -0,0 +1,95 @@
+function md=mesh(md,domainname,varargin)
+%MESH - create model mesh
+%
+%   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=mesh(md,domainname,resolution)
+%   or md=mesh(md,domainname,riftname, resolution)
+%
+%   Examples:
+%      md=mesh(md,'DomainOutline.exp',1000);
+%      md=mesh(md,'DomainOutline.exp','Rifts.exp',1500);
+
+%Figure out a characteristic area. Resolution is a grid oriented concept (ex a 1000m  resolution grid would 
+%be made of 1000*1000 area squares). 
+if (nargin==3),
+	resolution=varargin{1};
+	riftname='';
+end
+if (nargin==4),
+	riftname=varargin{1};
+	resolution=varargin{2};
+end
+
+%Check that mesh was not already run, and warn user: 
+if subsref(md,struct('type','.','subs','counter'))>=1,
+	choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s');
+	if ~strcmp(choice,'y')
+		error('no meshing done ... exiting');
+	end
+end
+
+area=resolution^2;
+
+%Mesh using TriMesh
+if strcmp(riftname,''),
+	[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes');
+else
+	[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes');
+
+	%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=readfile(domainname);
+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: /issm/trunk/src/m/classes/public/mesh/meshadaptation.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshadaptation.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/meshadaptation.m	(revision 1214)
@@ -0,0 +1,100 @@
+function md2=meshadaptation(md,field,scale,epsilon)
+%MESHADAPTATION - remesh a model for a given field
+%
+%   remesh a model for a given field, use the hessian computation to spread the error on the element
+%   scale: 1.2 large scale leads to higher refining.
+%   epsilon: .5  large epsilon will refine fewer elements
+%
+%   Usage:
+%      md2=meshadaptation(md,field,scale,epsilon)
+%
+%   Example:
+%      md2=meshadaptation(md,md.vel,1.2,0.5);
+
+%some checks
+if ~strcmpi(md.type,'2d')
+	error('meshadaptation error message: adaptation for 3d meshes not implemented yet')
+end
+if length(field)~=md.numberofgrids
+	error('meshadaptation error message: input field length shoud be numberofgrids')
+end
+
+disp(sprintf('      metric computation') )
+
+%load some variables (it is much faster if the variab;es are loaded from md once for all) 
+numberofelements=md.numberofelements;
+numberofgrids=md.numberofgrids;
+index=md.elements;
+x=md.x; y=md.y; z=md.z;
+
+%initialization
+alpha=zeros(md.numberofelements,3);
+beta=zeros(md.numberofelements,3);
+gradx=zeros(md.numberofgrids,1);
+grady=zeros(md.numberofgrids,1);
+metric=zeros(md.numberofelements,1);
+
+%build some usefull variables
+field=field(:);
+line=index(:);
+summation=1/3*ones(3,1);
+linesize=3*numberofelements;
+x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
+
+%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
+invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
+alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
+beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
+
+%Compute gradient for each element
+grad_elx=sum(field(index).*alpha,2); 
+grad_ely=sum(field(index).*beta,2);
+
+%compute the volume of each element
+areas=area(md);
+
+%update weights that holds the volume of all the element holding the grid i
+weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
+
+%Compute gradient for each grid (average of the elements around)
+gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
+grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
+gradx=gradx./weights;
+grady=grady./weights;
+
+%Compute hessian for each element
+hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*beta,2) sum(grady(index).*alpha,2)];
+
+%Compute the new metric
+hmin=min(sqrt(areas))/scale;
+hmax=max(sqrt(areas));
+
+%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,3); hessian(i,3)  hessian(i,2)]
+a=hessian(:,1); b=hessian(:,3); d=hessian(:,2);
+lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
+lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
+
+%Modify the eigen values to control the shape of the elements
+lambda1=min(max(abs(lambda1)/epsilon,1/hmax^2),1/hmin^2);
+lambda2=min(max(abs(lambda2)/epsilon,1/hmax^2),1/hmin^2);
+
+%find the corresponding eigen vectors of the initial eigen values
+%NOT USED FOR NOW
+%div1=sqrt(1+(b./(lambda1-d)).^2);
+%div2=sqrt(1+(b./(lambda2-d)).^2);
+%vector1=[1./div1 (b./(lambda1-d))./div1];
+%vector2=[1./div2 (b./(lambda2-d))./div2];
+
+%compute the metric
+%R=[vector1(:,1) vector2(:,1) vector1(:,2) vector2(:,2)];
+%det=vector1(:,1).*vector2(:,2)-vector1(:,2).*vector1(:,2);
+%Rinv=[vector2(:,2)./det   -vector2(:,1)./det   -vector1(:,2)./det   vector1(:,1)./det];
+metric=20*1./sqrt(lambda1.^2+lambda2.^2);
+
+%ok, metric can't go under 1km resolution.
+metric(find(metric<10^6))=10^6;
+
+%Remesh with this new metric
+disp(sprintf('      initial number of element: %i', md.numberofelements))
+md2=meshrefine(md,full(metric));
+disp(sprintf('      new number of elements:    %i', md2.numberofelements))
Index: /issm/trunk/src/m/classes/public/mesh/meshexprefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshexprefine.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/meshexprefine.m	(revision 1214)
@@ -0,0 +1,59 @@
+function md=meshexprefine(md,domainname,newresolution)
+%MESHEXPREFINE - refine mesh from a model
+%
+%   This routine refines a mesh for a given area using an Argus domain outline 
+%   and a new resolution for this domain.
+%   Appending ~ at the beginning of domainname will make this routine refine
+%   the elements which are not in the domain outline. 
+%   newresolution is the resolution to which the flagged elements must be refined.
+%
+%   Usage:
+%      md=meshexprefine(md,domainname,newresolution)
+%
+%   Example:
+%      md=meshexprefine(md,'RefineAreas.exp',1000)
+
+%some checks on list of arguments
+if ((nargin~=3) | (nargout~=1)),
+	meshexprefineusage();
+	error('meshexprefine error message');
+end
+
+if ~ischar(domainname) 
+	meshexprefineusage();
+	error('meshexprefine error message');
+end
+
+if (md.counter<1)
+	error('meshexprefine error message: you need to run mesh.m first on this model');
+end
+
+%Check the first letter of domainname, if it is ~, then we need to strip it away from 
+%domainname. The subsequent refinement will be done on the mask of the domain outline.
+if strcmpi(domainname(1),'~'),
+	mask=1;
+	domainname=domainname(2:length(domainname));
+else
+	mask=0;
+end
+
+%Read domainname file into a matlab array (x,y):
+refinearea=ContourToMesh(md.elements,md.x,md.y,expread(domainname,1),'element',1);
+aires=area(md);
+
+%flags areas within the domain
+if mask==1,
+	pos=find(~refinearea);
+else
+	pos=find(refinearea);
+end
+aires(pos)=newresolution*newresolution/2; %triangle area
+
+pos=find(refinearea);
+aires(pos)=2*aires(pos); %be sure this does not get refined.
+
+%refine using the new area vector
+md=meshrefine(md,aires);
+
+%return model
+end
Index: /issm/trunk/src/m/classes/public/mesh/meshprocessrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshprocessrifts.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/meshprocessrifts.m	(revision 1214)
@@ -0,0 +1,31 @@
+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)),
+	meshprocessriftsusage();
+	error('meshprocessrifts error message');
+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);
+
+function meshprocessriftsusage(),
+disp('usage: md=meshprocessrifts(md)');
Index: /issm/trunk/src/m/classes/public/mesh/meshrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshrefine.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/meshrefine.m	(revision 1214)
@@ -0,0 +1,40 @@
+function md=meshrefine(md,areas)
+%MESHREFINE:  refined the mesh from a model, according to an area metric.
+%
+%   Usage:
+%      md=meshrefine(md,metric)
+
+%some checks on list of arguments
+if ((nargin~=2) | (nargout~=1)),
+	meshrefineusage();
+	error('meshrefine error message');
+end
+if ( (isempty(areas)) |  (length(areas)~=md.numberofelements) | (length(find(isnan(areas))))),
+	meshrefineusage();
+	error('meshrefine error message');
+end
+
+%Refine using TriMeshRefine
+%[md.elements,md.x,md.y,md.segments,md.rifts]=TriMeshRefine(md.elements,md.x,md.y,md.segments,md.rifts,areas,'yes');
+[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMeshRefine(md.elements,md.x,md.y,md.segments,md.segmentmarkers,areas,'yes');
+
+%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);
+
+%type of model
+md.type='2d';
+	
+%augment counter  keeping track of what has been done to this model
+md.counter=1;
+end
+
+function meshrefineusage(),
+disp('usage: md=meshrefine(md,areas)');
+end
Index: /issm/trunk/src/m/classes/public/mesh/rifttipsrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/rifttipsrefine.m	(revision 1214)
+++ /issm/trunk/src/m/classes/public/mesh/rifttipsrefine.m	(revision 1214)
@@ -0,0 +1,33 @@
+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');
Index: sm/trunk/src/m/classes/public/meshadaptation.m
===================================================================
--- /issm/trunk/src/m/classes/public/meshadaptation.m	(revision 1213)
+++ 	(revision )
@@ -1,100 +1,0 @@
-function md2=meshadaptation(md,field,scale,epsilon)
-%MESHADAPTATION - remesh a model for a given field
-%
-%   remesh a model for a given field, use the hessian computation to spread the error on the element
-%   scale: 1.2 large scale leads to higher refining.
-%   epsilon: .5  large epsilon will refine fewer elements
-%
-%   Usage:
-%      md2=meshadaptation(md,field,scale,epsilon)
-%
-%   Example:
-%      md2=meshadaptation(md,md.vel,1.2,0.5);
-
-%some checks
-if ~strcmpi(md.type,'2d')
-	error('meshadaptation error message: adaptation for 3d meshes not implemented yet')
-end
-if length(field)~=md.numberofgrids
-	error('meshadaptation error message: input field length shoud be numberofgrids')
-end
-
-disp(sprintf('      metric computation') )
-
-%load some variables (it is much faster if the variab;es are loaded from md once for all) 
-numberofelements=md.numberofelements;
-numberofgrids=md.numberofgrids;
-index=md.elements;
-x=md.x; y=md.y; z=md.z;
-
-%initialization
-alpha=zeros(md.numberofelements,3);
-beta=zeros(md.numberofelements,3);
-gradx=zeros(md.numberofgrids,1);
-grady=zeros(md.numberofgrids,1);
-metric=zeros(md.numberofelements,1);
-
-%build some usefull variables
-field=field(:);
-line=index(:);
-summation=1/3*ones(3,1);
-linesize=3*numberofelements;
-x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
-
-%compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
-invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
-alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
-beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
-
-%Compute gradient for each element
-grad_elx=sum(field(index).*alpha,2); 
-grad_ely=sum(field(index).*beta,2);
-
-%compute the volume of each element
-areas=area(md);
-
-%update weights that holds the volume of all the element holding the grid i
-weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
-
-%Compute gradient for each grid (average of the elements around)
-gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
-grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
-gradx=gradx./weights;
-grady=grady./weights;
-
-%Compute hessian for each element
-hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*beta,2) sum(grady(index).*alpha,2)];
-
-%Compute the new metric
-hmin=min(sqrt(areas))/scale;
-hmax=max(sqrt(areas));
-
-%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,3); hessian(i,3)  hessian(i,2)]
-a=hessian(:,1); b=hessian(:,3); d=hessian(:,2);
-lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
-lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
-
-%Modify the eigen values to control the shape of the elements
-lambda1=min(max(abs(lambda1)/epsilon,1/hmax^2),1/hmin^2);
-lambda2=min(max(abs(lambda2)/epsilon,1/hmax^2),1/hmin^2);
-
-%find the corresponding eigen vectors of the initial eigen values
-%NOT USED FOR NOW
-%div1=sqrt(1+(b./(lambda1-d)).^2);
-%div2=sqrt(1+(b./(lambda2-d)).^2);
-%vector1=[1./div1 (b./(lambda1-d))./div1];
-%vector2=[1./div2 (b./(lambda2-d))./div2];
-
-%compute the metric
-%R=[vector1(:,1) vector2(:,1) vector1(:,2) vector2(:,2)];
-%det=vector1(:,1).*vector2(:,2)-vector1(:,2).*vector1(:,2);
-%Rinv=[vector2(:,2)./det   -vector2(:,1)./det   -vector1(:,2)./det   vector1(:,1)./det];
-metric=20*1./sqrt(lambda1.^2+lambda2.^2);
-
-%ok, metric can't go under 1km resolution.
-metric(find(metric<10^6))=10^6;
-
-%Remesh with this new metric
-disp(sprintf('      initial number of element: %i', md.numberofelements))
-md2=meshrefine(md,full(metric));
-disp(sprintf('      new number of elements:    %i', md2.numberofelements))
Index: sm/trunk/src/m/classes/public/meshexprefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/meshexprefine.m	(revision 1213)
+++ 	(revision )
@@ -1,59 +1,0 @@
-function md=meshexprefine(md,domainname,newresolution)
-%MESHEXPREFINE - refine mesh from a model
-%
-%   This routine refines a mesh for a given area using an Argus domain outline 
-%   and a new resolution for this domain.
-%   Appending ~ at the beginning of domainname will make this routine refine
-%   the elements which are not in the domain outline. 
-%   newresolution is the resolution to which the flagged elements must be refined.
-%
-%   Usage:
-%      md=meshexprefine(md,domainname,newresolution)
-%
-%   Example:
-%      md=meshexprefine(md,'RefineAreas.exp',1000)
-
-%some checks on list of arguments
-if ((nargin~=3) | (nargout~=1)),
-	meshexprefineusage();
-	error('meshexprefine error message');
-end
-
-if ~ischar(domainname) 
-	meshexprefineusage();
-	error('meshexprefine error message');
-end
-
-if (md.counter<1)
-	error('meshexprefine error message: you need to run mesh.m first on this model');
-end
-
-%Check the first letter of domainname, if it is ~, then we need to strip it away from 
-%domainname. The subsequent refinement will be done on the mask of the domain outline.
-if strcmpi(domainname(1),'~'),
-	mask=1;
-	domainname=domainname(2:length(domainname));
-else
-	mask=0;
-end
-
-%Read domainname file into a matlab array (x,y):
-refinearea=ContourToMesh(md.elements,md.x,md.y,expread(domainname,1),'element',1);
-aires=area(md);
-
-%flags areas within the domain
-if mask==1,
-	pos=find(~refinearea);
-else
-	pos=find(refinearea);
-end
-aires(pos)=newresolution*newresolution/2; %triangle area
-
-pos=find(refinearea);
-aires(pos)=2*aires(pos); %be sure this does not get refined.
-
-%refine using the new area vector
-md=meshrefine(md,aires);
-
-%return model
-end
Index: sm/trunk/src/m/classes/public/meshprocessrifts.m
===================================================================
--- /issm/trunk/src/m/classes/public/meshprocessrifts.m	(revision 1213)
+++ 	(revision )
@@ -1,31 +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)),
-	meshprocessriftsusage();
-	error('meshprocessrifts error message');
-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);
-
-function meshprocessriftsusage(),
-disp('usage: md=meshprocessrifts(md)');
Index: sm/trunk/src/m/classes/public/meshrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/meshrefine.m	(revision 1213)
+++ 	(revision )
@@ -1,40 +1,0 @@
-function md=meshrefine(md,areas)
-%MESHREFINE:  refined the mesh from a model, according to an area metric.
-%
-%   Usage:
-%      md=meshrefine(md,metric)
-
-%some checks on list of arguments
-if ((nargin~=2) | (nargout~=1)),
-	meshrefineusage();
-	error('meshrefine error message');
-end
-if ( (isempty(areas)) |  (length(areas)~=md.numberofelements) | (length(find(isnan(areas))))),
-	meshrefineusage();
-	error('meshrefine error message');
-end
-
-%Refine using TriMeshRefine
-%[md.elements,md.x,md.y,md.segments,md.rifts]=TriMeshRefine(md.elements,md.x,md.y,md.segments,md.rifts,areas,'yes');
-[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMeshRefine(md.elements,md.x,md.y,md.segments,md.segmentmarkers,areas,'yes');
-
-%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);
-
-%type of model
-md.type='2d';
-	
-%augment counter  keeping track of what has been done to this model
-md.counter=1;
-end
-
-function meshrefineusage(),
-disp('usage: md=meshrefine(md,areas)');
-end
Index: sm/trunk/src/m/classes/public/rifttipsrefine.m
===================================================================
--- /issm/trunk/src/m/classes/public/rifttipsrefine.m	(revision 1213)
+++ 	(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');