Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1337)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1337)
@@ -0,0 +1,112 @@
+function md=meshyams(md,velpath,domainoutline,varargin);
+%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
+%
+%   Usage:
+%      md=meshyams;
+
+%PIG settings
+%{
+nsteps=4;
+resolution=10000;
+hmin=1500;
+hmax=10^7;
+gradation=3*ones(nsteps,1);
+epsilon=10^-0;
+scale=1;
+%}
+
+%%{
+%Antarctica settings
+nsteps=6;
+resolution=5000;
+hmin=300;        %300m
+hmax=150*10^3;   %150km
+gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
+epsilon=2*10^-0; %3m/a interpolation error
+scale=1;
+%%}
+
+%check number of arguments
+waterflag=0;
+if nargin==4,
+	groundedoutline=varargin{1};
+	if exist(groundedoutline);
+		disp(['grounded ice domain found. Metric will be minimum on water']);
+		waterflag=1;
+	else
+		error(['MeshYams error message: file ' groundedoutline ' not found.']);
+	end
+end
+
+%mesh with initial resolution
+md=mesh(md,domainoutline,resolution);
+
+%load velocities 
+load(velpath);
+
+disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);
+
+%start mesh adaptation
+for i=1:nsteps,
+	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
+
+	%interpolate velocities onto mesh
+	disp('   interpolating velocities');
+	md.vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0)*scale;
+	md.vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0)*scale;
+	md.vel_obs=averaging(md,sqrt(md.vx_obs.^2+md.vy_obs.^2),2);
+
+	%set gridonwater field
+	if waterflag,
+		gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'node',2);
+		md.gridonwater=ones(md.numberofgrids,1);
+		md.gridonwater(find(gridground))=0;
+	else
+		md.gridonwater=zeros(md.numberofgrids,1);
+	end
+
+	%adapt according to velocities
+	disp('   adapting');
+	md=YamsCall(md,md.vel_obs,hmin,hmax,gradation(i),epsilon);
+end
+	
+disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
+
+%Now, build the connectivity tables for this mesh.
+md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
+md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
+
+%Recreate the segments
+elementconnectivity=md.elementconnectivity;
+elementonboundary=double(elementconnectivity(:,end)==0);
+pos=find(elementonboundary);
+num_segments=length(pos);
+segments=zeros(num_segments,3);
+for i=1:num_segments,
+	el1=pos(i);
+	els2=elementconnectivity(el1,find(elementconnectivity(el1,:)));
+	flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:));
+	nods1=md.elements(el1,:);
+	nods1(find(nods1==flag))=[];
+	segments(i,:)=[nods1 el1];
+
+	ord1=find(nods1(1)==md.elements(el1,:));
+	ord2=find(nods1(2)==md.elements(el1,:));
+
+	%swap segment grids if necessary
+	if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
+		temp=segments(i,1);
+		segments(i,1)=segments(i,2);
+		segments(i,2)=temp;
+	end
+	segments(i,1:2)=fliplr(segments(i,1:2));
+end
+md.segments=segments;
+
+%Fill in rest of fields:
+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);
Index: sm/trunk/src/m/utils/Mesh/MeshYams.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/MeshYams.m	(revision 1336)
+++ 	(revision )
@@ -1,112 +1,0 @@
-function md=MeshYams(md,velpath,domainoutline,varargin);
-%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
-%
-%   Usage:
-%      md=MeshYams;
-
-%PIG settings
-%{
-nsteps=4;
-resolution=10000;
-hmin=1500;
-hmax=10^7;
-gradation=3*ones(nsteps,1);
-epsilon=10^-0;
-scale=1;
-%}
-
-%%{
-%Antarctica settings
-nsteps=6;
-resolution=5000;
-hmin=300;        %300m
-hmax=150*10^3;   %150km
-gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=2*10^-0; %3m/a interpolation error
-scale=1;
-%%}
-
-%check number of arguments
-waterflag=0;
-if nargin==4,
-	groundedoutline=varargin{1};
-	if exist(groundedoutline);
-		disp(['grounded ice domain found. Metric will be minimum on water']);
-		waterflag=1;
-	else
-		error(['MeshYams error message: file ' groundedoutline ' not found.']);
-	end
-end
-
-%mesh with initial resolution
-md=mesh(md,domainoutline,resolution);
-
-%load velocities 
-load(velpath);
-
-disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);
-
-%start mesh adaptation
-for i=1:nsteps,
-	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
-
-	%interpolate velocities onto mesh
-	disp('   interpolating velocities');
-	md.vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0)*scale;
-	md.vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0)*scale;
-	md.vel_obs=averaging(md,sqrt(md.vx_obs.^2+md.vy_obs.^2),2);
-
-	%set gridonwater field
-	if waterflag,
-		gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'node',2);
-		md.gridonwater=ones(md.numberofgrids,1);
-		md.gridonwater(find(gridground))=0;
-	else
-		md.gridonwater=zeros(md.numberofgrids,1);
-	end
-
-	%adapt according to velocities
-	disp('   adapting');
-	md=YamsCall(md,md.vel_obs,hmin,hmax,gradation(i),epsilon);
-end
-	
-disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
-
-%Now, build the connectivity tables for this mesh.
-md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids);
-md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
-
-%Recreate the segments
-elementconnectivity=md.elementconnectivity;
-elementonboundary=double(elementconnectivity(:,end)==0);
-pos=find(elementonboundary);
-num_segments=length(pos);
-segments=zeros(num_segments,3);
-for i=1:num_segments,
-	el1=pos(i);
-	els2=elementconnectivity(el1,find(elementconnectivity(el1,:)));
-	flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:));
-	nods1=md.elements(el1,:);
-	nods1(find(nods1==flag))=[];
-	segments(i,:)=[nods1 el1];
-
-	ord1=find(nods1(1)==md.elements(el1,:));
-	ord2=find(nods1(2)==md.elements(el1,:));
-
-	%swap segment grids if necessary
-	if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
-		temp=segments(i,1);
-		segments(i,1)=segments(i,2);
-		segments(i,2)=temp;
-	end
-	segments(i,1:2)=fliplr(segments(i,1:2));
-end
-md.segments=segments;
-
-%Fill in rest of fields:
-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);