Index: sm/trunk-jpl/src/m/model/mesh/meshadaptation.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/meshadaptation.m	(revision 11014)
+++ 	(revision )
@@ -1,94 +1,0 @@
-function md2=meshadaptation(md,field,scale,epsilon,threshold)
-%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
-%   threshold: value under which no refinement is done.
-%
-%   Usage:
-%      md2=meshadaptation(md,field,scale,epsilon,threshold)
-%
-%   Example:
-%      md2=meshadaptation(md,md.initialization.vel,1.2,0.5,1);
-
-%some checks
-if ~(md.mesh.dimension==2)
-	error('meshadaptation error message: adaptation for 3d meshes not implemented yet')
-end
-if length(field)~=md.mesh.numberofvertices
-	error('meshadaptation error message: input field length shoud be numberofnodes')
-end
-
-disp(sprintf('      metric computation') )
-
-%initialization
-index=md.mesh.elements;
-numberofnodes=md.mesh.numberofvertices;
-numberofelements=md.mesh.numberofelements;
-gradx=zeros(numberofnodes,1);
-grady=zeros(numberofnodes,1);
-metric=zeros(numberofelements,1);
-
-%take the threshold, a characteristic length, and get the area equivalent threshold:
-threshold=threshold^2/2;
-
-%build some usefull variables
-field=field(:);
-line=index(:);
-summation=1/3*ones(3,1);
-linesize=3*numberofelements;
-
-%get areas and  nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
-[alpha beta]=GetNodalFunctionsCoeff(index,md.mesh.x,md.mesh.y);
-areas=GetAreas(index,md.mesh.x,md.mesh.y);
-
-%Compute gradient for each element
-grad_elx=sum(field(index).*alpha,2); 
-grad_ely=sum(field(index).*beta,2);
-
-%update weights that holds the volume of all the element holding the node i
-weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofnodes,1);
-
-%Compute gradient for each node (average of the elements around)
-gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofnodes,1);
-grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofnodes,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<threshold))=threshold;
-
-%Remesh with this new metric
-disp(sprintf('      initial number of element: %i', md.mesh.numberofelements))
-md2=meshrefine(md,full(metric));
-disp(sprintf('      new number of elements:    %i', md2.mesh.numberofelements))
Index: sm/trunk-jpl/src/m/model/mesh/meshbamg.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/meshbamg.m	(revision 11014)
+++ 	(revision )
@@ -1,112 +1,0 @@
-function md=meshbamg(md,varargin);
-%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
-%
-%   Usage:
-%      md=meshyams(md,varargin);
-%      where varargin is a lit of paired arguments. 
-%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
-%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
-%      arguments can be: 'thickness': matlab file containing the thickness [m]
-%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
-%                          this option is used to minimize the metric on water (no refinement)
-%      optional arguments: 'resolution': initial mesh resolution [m]
-%      optional arguments: 'nsteps': number of steps of mesh adaptation
-%      optional arguments: 'epsilon': average interpolation error wished [m/yr [m]]
-%      optional arguments: 'hmin': minimum edge length
-%      optional arguments: 'hmanx': maximum edge
-%      
-%
-%   Examples:
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','thickness','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
-
-%recover options
-options=pairoptions(varargin{:});
-options=deleteduplicates(options,1);
-
-%recover some fields
-disp('MeshBamg Options:')
-domainoutline=getfieldvalue(options,'domainoutline');
-disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
-groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
-disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
-velocities=getfieldvalue(options,'velocities');
-disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
-thickness=getfieldvalue(options,'thickness','none');
-disp(sprintf('   %-15s: ''%s''','Thickness',thickness));
-resolution=getfieldvalue(options,'resolution',5000);
-disp(sprintf('   %-15s: %f','Resolution',resolution));
-nsteps=getfieldvalue(options,'nsteps',6);
-disp(sprintf('   %-15s: %i','nsteps',nsteps));
-gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
-disp(sprintf('   %-15s: %g','gradation',gradation(1)));
-epsilon=getfieldvalue(options,'epsilon',3);
-disp(sprintf('   %-15s: %f','epsilon',epsilon));
-hmin=getfieldvalue(options,'hmin',500);
-disp(sprintf('   %-15s: %f','hmin',hmin));
-hmax=getfieldvalue(options,'hmax',150*10^3);
-disp(sprintf('   %-15s: %f\n','hmax',hmax));
-
-%mesh with initial resolution
-disp('Initial mesh generation...');
-md=setmesh(md,domainoutline,resolution);
-disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
-
-%load velocities 
-disp('loading velocities...');
-NamesV=VelFindVarNames(velocities);
-Vel=load(velocities);
-
-%compute thickness?
-if ~strcmpi(thickness,'none'),
-	disp('loading thickness...');
-	NamesH=FieldFindVarNames(thickness);
-	Thi=load(thickness);
-	thicknesspresent=1;
-else
-	thicknesspresent=0;
-end
-
-%start mesh adaptation
-for i=1:nsteps,
-	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
-
-	%interpolate velocities onto mesh
-	disp('   interpolating velocities...');
-	if strcmpi(NamesV.interp,'node'),
-		vx_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.mesh.x,md.mesh.y,0);
-	else
-		vx_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.mesh.x,md.mesh.y,0);
-	end
-	field=sqrt(vx_obs.^2+vy_obs.^2);
-	
-	if thicknesspresent,
-		disp('   interpolating thickness...');
-		if strcmpi(NamesH.interp,'node'),
-			h=InterpFromGridToMesh(Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.mesh.x,md.mesh.y,0);
-		else
-			h=InterpFromMeshToMesh2d(Thi.(NamesH.indexname),Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.mesh.x,md.mesh.y,0);
-		end
-		field=[field h];
-	end
-	
-	%set nodeonwater field
-	if ~strcmp(groundeddomain,'N/A'),
-		nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
-		md.nodeonwater=ones(md.mesh.numberofvertices,1);
-		md.nodeonwater(find(nodeground))=0;
-	else
-		md.nodeonwater=zeros(md.mesh.numberofvertices,1);
-	end
-
-	%adapt according to velocities
-	disp('   adapting...');
-	md=bamg(md,'field',field,'hmin',hmin,'hmax',hmax,'gradation',gradation(i),'err',epsilon);
-
-end
-	
-disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
Index: sm/trunk-jpl/src/m/model/mesh/meshexprefine.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/meshexprefine.m	(revision 11014)
+++ 	(revision )
@@ -1,55 +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
-
-%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 domainame file into a matlab array (x,y):
-refinearea=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,domainname,'element',1);
-aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y);
-
-%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-jpl/src/m/model/mesh/meshyams.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/meshyams.m	(revision 11014)
+++ 	(revision )
@@ -1,172 +1,0 @@
-function md=meshyams(md,varargin);
-%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
-%
-%   Usage:
-%      md=meshyams(md,varargin);
-%      where varargin is a lit of paired arguments. 
-%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
-%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
-%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
-%                          this option is used to minimize the metric on water (no refinement)
-%      optional arguments: 'resolution': initial mesh resolution [m]
-%      optional arguments: 'nsteps': number of steps of mesh adaptation
-%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
-%      optional arguments: 'hmin': minimum edge length
-%      optional arguments: 'hmanx': maximum edge
-%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
-%      
-%
-%   Examples:
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
-%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
-
-%recover options
-options=pairoptions(varargin{:});
-options=deleteduplicates(options,1);
-
-%recover some fields
-disp('MeshYams Options:')
-domainoutline=getfieldvalue(options,'domainoutline');
-disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
-riftoutline=getfieldvalue(options,'riftoutline','N/A');
-disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
-groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
-disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
-velocities=getfieldvalue(options,'velocities');
-disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
-resolution=getfieldvalue(options,'resolution',5000);
-disp(sprintf('   %-15s: %f','Resolution',resolution));
-nsteps=getfieldvalue(options,'nsteps',6);
-disp(sprintf('   %-15s: %i','nsteps',nsteps));
-gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
-disp(sprintf('   %-15s: %g','gradation',gradation(1)));
-epsilon=getfieldvalue(options,'epsilon',3);
-disp(sprintf('   %-15s: %f','epsilon',epsilon));
-hmin=getfieldvalue(options,'hmin',500);
-disp(sprintf('   %-15s: %f','hmin',hmin));
-hmax=getfieldvalue(options,'hmax',150*10^3);
-disp(sprintf('   %-15s: %f\n','hmax',hmax));
-
-%mesh with initial resolution
-disp('Initial mesh generation...');
-if strcmpi(riftoutline,'N/A');
-	md=setmesh(md,domainoutline,resolution);
-else
-	md=setmesh(md,domainoutline,riftoutline,resolution);
-	md=meshprocessrifts(md,domainoutline);
-end
-disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
-
-%load velocities 
-disp('loading velocities...');
-Names=VelFindVarNames(velocities);
-Vel=load(velocities);
-
-%start mesh adaptation
-for i=1:nsteps,
-	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
-
-	%interpolate velocities onto mesh
-	disp('   interpolating velocities...');
-	if strcmpi(Names.interp,'node'),
-		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-	else
-		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-	end
-	field=sqrt(vx_obs.^2+vy_obs.^2);
-
-	%set mask.vertexonwater  field
-	if ~strcmp(groundeddomain,'N/A'),
-		nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
-		md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
-		md.mask.vertexonwater(find(nodeground))=0;
-	else
-		md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
-	end
-
-	%adapt according to velocities
-	disp('   adapting...');
-	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.rifts.numrifts, 
-		md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
-		md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
-		md.mesh.segments=findsegments(md);
-		md=meshyamsrecreateriftsegments(md);
-	end
-
-end
-	
-disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
-
-%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);
-
-%recreate segments
-md.mesh.segments=findsegments(md);
-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
-
-%Fill in rest of fields:
-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);
-if ~strcmp(groundeddomain,'N/A'),
-	nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
-	md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
-	md.mask.vertexonwater(find(nodeground))=0;
-else
-	md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
-end
-if strcmpi(Names.interp,'node'),
-	md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-	md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-else
-	md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
-	md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
-end
-md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
-
-%deal with rifts 
-if md.rifts.numrifts,
-	%first, recreate rift segments
-	md=meshyamsrecreateriftsegments(md);
-
-	%using the segments, recreate the penaltypairs
-	for j=1:md.rifts.numrifts,
-		rift=md.rifts.riftstruct(j);
-
-		%build normals and lengths of segments:
-		lengths=sqrt((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))).^2 + (md.mesh.y(rift.segments(:,1))-md.mesh.y(rift.segments(:,2))).^2 );
-		normalsx=cos(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
-		normalsy=sin(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.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.riftstruct(j)=rift;
-	end
-
-end
Index: sm/trunk-jpl/src/m/model/mesh/reorder.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/reorder.m	(revision 11014)
+++ 	(revision )
@@ -1,27 +1,0 @@
-function md=reorder(md)
-%REORDER - reorder nodes and elements of a model
-%
-%   Usage: 
-%      md=reorder(md)
-
-%some checks
-if md.mesh.dimension==3,
-	error('reorder error message: 3d models not supported yet, exiting...')
-end
-disp('reorder warning: only the mesh fields are reorder. The model needs to be reparameterized');
-
-%reorder elements
-newelements=randperm(md.mesh.numberofelements)';
-tnewelements=zeros(md.mesh.numberofelements,1);tnewelements(newelements)=[1:md.mesh.numberofelements]';
-
-%reorder nodes
-newnodes=randperm(md.mesh.numberofvertices)';
-tnewnodes=zeros(md.mesh.numberofvertices,1);tnewnodes(newnodes)=[1:md.mesh.numberofvertices]';
-
-%update all fields
-md.mesh.elements=tnewnodes(md.mesh.elements(newelements,:));
-md.mesh.segments=[tnewnodes(md.mesh.segments(:,1)) tnewnodes(md.mesh.segments(:,2)) tnewelements(md.mesh.segments(:,3))];
-md.mesh.x=md.mesh.x(newnodes);
-md.mesh.y=md.mesh.y(newnodes);
-md.mesh.z=md.mesh.z(newnodes);
-md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
Index: /issm/trunk-jpl/src/m/model/mesh/yams.m
===================================================================
--- /issm/trunk-jpl/src/m/model/mesh/yams.m	(revision 11015)
+++ /issm/trunk-jpl/src/m/model/mesh/yams.m	(revision 11015)
@@ -0,0 +1,171 @@
+function md=yams(md,varargin);
+%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
+%
+%   Usage:
+%      md=yams(md,varargin);
+%      where varargin is a lit of paired arguments. 
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
+%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
+%                          this option is used to minimize the metric on water (no refinement)
+%      optional arguments: 'resolution': initial mesh resolution [m]
+%      optional arguments: 'nsteps': number of steps of mesh adaptation
+%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
+%      optional arguments: 'hmin': minimum edge length
+%      optional arguments: 'hmanx': maximum edge
+%      optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
+%      
+%
+%   Examples:
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat');
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
+%      md=yams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
+
+%recover options
+options=pairoptions(varargin{:});
+options=deleteduplicates(options,1);
+
+%recover some fields
+disp('MeshYams Options:')
+domainoutline=getfieldvalue(options,'domainoutline');
+disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
+riftoutline=getfieldvalue(options,'riftoutline','N/A');
+disp(sprintf('   %-15s: ''%s''','riftoutline',riftoutline));
+groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
+disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
+velocities=getfieldvalue(options,'velocities');
+disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+resolution=getfieldvalue(options,'resolution',5000);
+disp(sprintf('   %-15s: %f','Resolution',resolution));
+nsteps=getfieldvalue(options,'nsteps',6);
+disp(sprintf('   %-15s: %i','nsteps',nsteps));
+gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
+disp(sprintf('   %-15s: %g','gradation',gradation(1)));
+epsilon=getfieldvalue(options,'epsilon',3);
+disp(sprintf('   %-15s: %f','epsilon',epsilon));
+hmin=getfieldvalue(options,'hmin',500);
+disp(sprintf('   %-15s: %f','hmin',hmin));
+hmax=getfieldvalue(options,'hmax',150*10^3);
+disp(sprintf('   %-15s: %f\n','hmax',hmax));
+
+%mesh with initial resolution
+disp('Initial mesh generation...');
+if strcmpi(riftoutline,'N/A');
+	md=setmesh(md,domainoutline,resolution);
+else
+	md=setmesh(md,domainoutline,riftoutline,resolution);
+	md=meshprocessrifts(md,domainoutline);
+end
+disp(['Initial mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
+
+%load velocities 
+disp('loading velocities...');
+Names=VelFindVarNames(velocities);
+Vel=load(velocities);
+
+%start mesh adaptation
+for i=1:nsteps,
+	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
+
+	%interpolate velocities onto mesh
+	disp('   interpolating velocities...');
+	if strcmpi(Names.interp,'node'),
+		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+	else
+		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+	end
+	field=sqrt(vx_obs.^2+vy_obs.^2);
+
+	%set mask.vertexonwater  field
+	if ~strcmp(groundeddomain,'N/A'),
+		nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
+		md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
+		md.mask.vertexonwater(find(nodeground))=0;
+	else
+		md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
+	end
+
+	%adapt according to velocities
+	disp('   adapting...');
+	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.rifts.numrifts, 
+		md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+		md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+		md.mesh.segments=findsegments(md);
+		md=yamsrecreateriftsegments(md);
+	end
+
+end
+	
+disp(['Final mesh, number of elements: ' num2str(md.mesh.numberofelements)]);
+
+%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);
+
+%recreate segments
+md.mesh.segments=findsegments(md);
+md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+
+%Fill in rest of fields:
+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);
+if ~strcmp(groundeddomain,'N/A'),
+	nodeground=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,groundeddomain,'node',2);
+	md.mask.vertexonwater=ones(md.mesh.numberofvertices,1);
+	md.mask.vertexonwater(find(nodeground))=0;
+else
+	md.mask.vertexonwater=zeros(md.mesh.numberofvertices,1);
+end
+if strcmpi(Names.interp,'node'),
+	md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+	md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+else
+	md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.mesh.x,md.mesh.y,0);
+	md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.mesh.x,md.mesh.y,0);
+end
+md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
+
+%deal with rifts 
+if md.rifts.numrifts,
+	%first, recreate rift segments
+	md=meshyamsrecreateriftsegments(md);
+
+	%using the segments, recreate the penaltypairs
+	for j=1:md.rifts.numrifts,
+		rift=md.rifts.riftstruct(j);
+
+		%build normals and lengths of segments:
+		lengths=sqrt((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))).^2 + (md.mesh.y(rift.segments(:,1))-md.mesh.y(rift.segments(:,2))).^2 );
+		normalsx=cos(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.y(rift.segments(:,1)))));
+		normalsy=sin(atan2((md.mesh.x(rift.segments(:,1))-md.mesh.x(rift.segments(:,2))) , (md.mesh.y(rift.segments(:,2))-md.mesh.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.riftstruct(j)=rift;
+	end
+end