Index: /issm/trunk/src/m/classes/public/mesh/meshquality.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshquality.m	(revision 1341)
+++ /issm/trunk/src/m/classes/public/mesh/meshquality.m	(revision 1341)
@@ -0,0 +1,164 @@
+function meshquality(md);
+%MESHQUALITY - compute mesh quality
+%
+%   Usage:
+%      meshquality(md);
+
+t1=clock; fprintf('%s','      computing Hessian...');
+%initialization
+field=md.vel_obs;
+x=md.x;
+y=md.y;
+index=md.elements;
+numberofgrids=md.numberofgrids;
+numberofelements=md.numberofelements;
+gradx=zeros(numberofgrids,1);
+grady=zeros(numberofgrids,1);
+metric=zeros(numberofelements,3);
+
+%2d geometric parameter (do not change)
+scale=2/9; 
+hmin=500;
+hmax=10^6;
+epsilon=5;
+
+%build some usefull variables
+line=index(:);
+linesize=3*numberofelements;
+
+%get areas and  nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
+[alpha beta]=GetNodalFunctionsCoeff(index,md.x,md.y);
+areas=GetAreas(index,md.x,md.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 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_elem=[sum(gradx(index).*alpha,2) sum(grady(index).*alpha,2) sum(grady(index).*beta,2)];
+
+%Compute Hessian on the nodes (average of the elements around)
+hessian=[sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,1),3,1),numberofgrids,1)./weights ...
+	sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,2),3,1),numberofgrids,1)./weights ...
+	sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,3),3,1),numberofgrids,1)./weights ];
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+t1=clock; fprintf('%s','      computing metric...');
+%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,2); hessian(i,2)  hessian(i,3)]
+a=hessian(:,1); b=hessian(:,2); d=hessian(:,3);
+lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
+lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
+pos1=find(lambda1==0);
+pos2=find(lambda2==0);
+
+%Modify the eigen values to control the shape of the elements
+lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
+lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
+
+%compute eigen vectors
+norm1=sqrt(8*b.^2+2*(d-a).^2+2*(d-a).*sqrt((a-d).^2+4*b.^2));
+v1x=2*b./norm1;
+v1y=((d-a)+sqrt((a-d).^2+4*b.^2))./norm1;
+norm2=sqrt(8*b.^2+2*(d-a).^2-2*(d-a).*sqrt((a-d).^2+4*b.^2));
+v2x=2*b./norm2;
+v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2;
+
+%Compute new metric (for each node M=V*Lambda*V^-1)
+metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ...
+	(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v1y.*v2y-lambda2.*v1y.*v2y) ...
+	(v1x.*v2y-v1y.*v2x).^(-1).*(-lambda1.*v2x.*v1y+lambda2.*v1x.*v2y)]);
+
+%some corrections for 0 eigen values
+metric(pos1,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos1),1);
+metric(pos2,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos2),1);
+if length(md.gridonwater)==numberofgrids;
+	pos=find(md.gridonwater);
+	metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+end
+
+%take care of water elements
+if length(md.gridonwater)==numberofgrids;
+	pos=find(md.gridonwater);
+	metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+end
+
+if any(isnan(metric)),
+	error('YamsCall error message: NaN in the metric...')
+end
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%Get Areas
+areas=GetAreas(index,x,y);
+
+%length edges vectors
+e1x=[x(index(:,2))-x(index(:,1))];
+e1y=[y(index(:,2))-y(index(:,1))];
+e2x=[x(index(:,3))-x(index(:,2))];
+e2y=[y(index(:,3))-y(index(:,2))];
+e3x=[x(index(:,1))-x(index(:,3))];
+e3y=[y(index(:,1))-y(index(:,3))];
+
+%metric of each the 3 grids for each element
+M1=metric(index(:,1),:);
+M2=metric(index(:,2),:);
+M3=metric(index(:,3),:);
+
+%Get edge length in the metric
+L1=1/2*(sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y))+sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y)));
+L2=1/2*(sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y))+sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y)));
+L3=1/2*(sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y))+sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y)));
+
+%area in the metric
+V=1/3*areas.*(sqrt(M1(:,1).*M1(:,3)-M1(:,2).^2)+sqrt(M2(:,1).*M2(:,3)-M2(:,2).^2)+sqrt(M3(:,1).*M3(:,3)-M3(:,2).^2));
+
+%compute quality:
+Q=4*sqrt(3)*V./(L1+L2+L3);
+
+%display
+X=0:0.1:4;
+hist(Q,X);
+xlim([0 3]);
+title('mesh quality distribution');
+
+%compute error
+%{
+%a=hessian_elem(:,1); b=hessian_elem(:,2); d=hessian_elem(:,3);
+a=hessian(:,1); b=hessian(:,2); d=hessian(:,3);
+a=a(index)*[1;1;1]/3;
+b=b(index)*[1;1;1]/3;
+d=d(index)*[1;1;1]/3;
+lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
+lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
+lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
+lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
+if length(md.elementonwater)==numberofelements;
+	pos=find(md.gridonwater);
+	lambda1(pos)=0;
+	lambda2(pos)=0;
+end
+%}
+if length(md.gridonwater)==numberofgrids;
+	lambda1(pos)=0;
+	lambda2(pos)=0;
+end
+lambda1=lambda1(index)*[1;1;1]/3;
+lambda2=lambda2(index)*[1;1;1]/3;
+
+lambdamax=max(lambda1,lambda2);
+hmax=max(max(sqrt(e1x.^2+e1y.^2),sqrt(e2x.^2+e2y.^2)),sqrt(e3x.^2+e3y.^2));
+epsilon=scale*hmax.^2.*lambdamax;
+
+%display
+plotmodel(md,'data',epsilon,'title','Interpolation error','figure',2,'log',10)
+
+disp(['Average interpolation error = ' num2str(mean(epsilon))  ' m/yr']);
+disp(['Maximum interpolation error = ' num2str(max(epsilon))  ' m/yr']);
Index: /issm/trunk/src/m/classes/public/mesh/meshyams.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1340)
+++ /issm/trunk/src/m/classes/public/mesh/meshyams.m	(revision 1341)
@@ -111,2 +111,12 @@
 md.elementonbed=ones(md.numberofelements,1);
 md.elementonsurface=ones(md.numberofelements,1);
+if ~strcmp(groundeddomain,'N/A'),
+	gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundeddomain,1),'node',2);
+	md.gridonwater=ones(md.numberofgrids,1);
+	md.gridonwater(find(gridground))=0;
+else
+	md.gridonwater=zeros(md.numberofgrids,1);
+end
+md.vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0);
+md.vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0);
+md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);