Index: /issm/trunk/src/m/classes/public/mesh/findsegments.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/findsegments.m	(revision 1346)
+++ /issm/trunk/src/m/classes/public/mesh/findsegments.m	(revision 1346)
@@ -0,0 +1,35 @@
+function segments=findsegments(md)
+%FINDSEGMENTS - build segments model field
+%
+%   Usage:
+%      segments=findsegments(md);
+
+%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
Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1345)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1346)
@@ -28,17 +28,19 @@
 
 %recover some fields
-domainoutline=yamsoptions.domainoutline;
-groundeddomain=yamsoptions.groundeddomain;
-velocities=yamsoptions.velocities;
-resolution=yamsoptions.resolution;
+disp('MeshYams Options:')
+domainoutline=yamsoptions.domainoutline;   disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
+groundeddomain=yamsoptions.groundeddomain; disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
+velocities=yamsoptions.velocities;         disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+resolution=yamsoptions.resolution;         disp(sprintf('   %-15s: %f','Resolution',resolution));
 gradation=yamsoptions.gradation;
-nsteps=yamsoptions.nsteps;
-epsilon=yamsoptions.epsilon;
-hmin=yamsoptions.hmin;
-hmax=yamsoptions.hmax;
+nsteps=yamsoptions.nsteps;                 disp(sprintf('   %-15s: %i','nsteps',nsteps));
+epsilon=yamsoptions.epsilon;               disp(sprintf('   %-15s: %f','epsilon',epsilon));
+hmin=yamsoptions.hmin;                     disp(sprintf('   %-15s: %f','hmin',hmin));
+hmax=yamsoptions.hmax;                     disp(sprintf('   %-15s: %f\n','hmax',hmax));
 
 %mesh with initial resolution
 disp('Initial mesh generation...');
 md=mesh(md,domainoutline,resolution);
+disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
 
 %load velocities 
@@ -51,5 +53,5 @@
 
 	%interpolate velocities onto mesh
-	disp('   interpolating velocities');
+	disp('   interpolating velocities...');
 	vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0);
 	vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0);
@@ -66,45 +68,17 @@
 
 	%adapt according to velocities
-	disp('   adapting');
+	disp('   adapting...');
 	md=YamsCall(md,field,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
-disp('Creating 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;
+%recreate segments
+md.segments=findsegments(md);
+md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
 
 %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);
Index: /issm/trunk/src/m/utils/Mesh/ComputeMetric.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/ComputeMetric.m	(revision 1345)
+++ /issm/trunk/src/m/utils/Mesh/ComputeMetric.m	(revision 1346)
@@ -15,4 +15,5 @@
 pos1=find(lambda1==0);
 pos2=find(lambda2==0);
+pos3=find(b==0 & lambda1==lambda2);
 
 %Modify the eigen values to control the shape of the elements
@@ -28,4 +29,7 @@
 v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2;
 
+v1x(pos3)=1; v1y(pos3)=0;
+v2x(pos3)=0; v2y(pos3)=1;
+
 %Compute new metric (for each node M=V*Lambda*V^-1)
 metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ...
@@ -40,5 +44,23 @@
 metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
 
+%take care of NaNs if any (use Matlab eig in a loop)
+[pos posj]=find(isnan(metric)); clear posj;
+if ~isempty(pos),
+	fprintf(' %i %s',length(pos),'NaN found in the metric. Use Matlab routine...');
+	for i=1:length(pos)
+		H=[hessian(pos(i),1) hessian(pos(i),2)
+		hessian(pos(i),2) hessian(pos(i),3)];
+		[u,v]=eig(H);
+		lambda1=v(1,1);
+		lambda2=v(2,2);
+		v(1,1)=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
+		v(2,2)=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
+
+		metricTria=u*v*u^(-1);
+		metric(pos(i),:)=[metricTria(1,1) metricTria(1,2) metricTria(2,2)];
+	end
+end
+
 if any(isnan(metric)),
-	error('ComputeMetric error message: NaN in the metric...')
+	error('ComputeMetric error message: NaN in the metric despite our efforts...')
 end
Index: /issm/trunk/src/m/utils/Mesh/YamsCall.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1345)
+++ /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1346)
@@ -10,10 +10,8 @@
 %
 %   Usage:
-%      md=YamsCall(md,field,hmin,hmax,epsilon)
+%      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
 %
 %   Example:
 %      md=YamsCall(md,md.vel_obs,1500,10^8,1.3,0.9);
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build metric %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 global ISSM_DIR
@@ -57,9 +55,4 @@
 fprintf(fid,'%i %i %i %i\n',[md.elements zeros(md.numberofelements,1)]');
 
-%Corners
-pos=find(md.gridonboundary);
-fprintf(fid,'\n\n%s\n%i\n\n','Corners',length(pos));
-fprintf(fid,'%i\n',pos');
-
 %close
 fclose(fid);
@@ -68,5 +61,5 @@
 %call yams
 fprintf('%s\n','      call Yams...');
-system(['yams2-linux -O 1 -v 2 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
 
 %plug new mesh