Index: /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m	(revision 1327)
+++ /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m	(revision 1327)
@@ -0,0 +1,82 @@
+function md=BuildAdaptedMeshYams
+%BUILDADAPTEDMESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
+%
+%   Usage:
+%      md=BuildAdaptedMeshYams
+
+%PIG settings
+%{
+nsteps=4;
+resolution=10000;
+hmin=1500;
+hmax=10^7;
+gradation=3*ones(nsteps,1);
+epsilon=10^-0;
+domainoutline='DomainOutline.exp';
+%}
+
+%%{
+%Antarctica settings
+nsteps=6;
+resolution=10000;
+hmin=500;
+hmax=3*10^5;
+gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
+epsilon=5*10^-1;
+domainoutline='EnvDomainOutline200km.exp';
+%%}
+
+velpath='/u/wilkes-r1b/larour/ModelData/RignotAntarcticaVelMosaicRampErsAlos/RignotAntVelSmooth.mat'; scale=1;
+%velpath='/u/wilkes-r1b/larour/ModelData/RignotAntarcticaVel1km/RignotAntVel.mat'; scale=365*24*3600; %(from m/s to m/yr)
+%velpath='/u/wilkes-r1b/larour/ModelData/RignotPigVel1996/RignotPig1996_HoleFree.mat'; scale=1;
+%velpath='/proj/ice/larour/Glaciology/Model_Data/RignotAntarcticaVel1km/RignotAntVel';
+%velpath='/proj/ice/larour/Glaciology/Model_Data/BamberAntarcticaVel1km/velocity';
+%velpath='/proj/ice/larour/Glaciology/Model_Data/mosaicAnt1km_ramp_ers_alos/RignotAntVel';
+
+groundingline='GroundingLine.exp';
+icestreams='Icestream.exp';
+parameterfile='Basins.par';
+
+%clear data
+clear md
+
+%build new model
+md=model;
+
+%mesh with initial resolution
+md=mesh(md,domainoutline,resolution);
+disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);
+
+%load velocities 
+load(velpath);
+
+%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);
+
+	%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)]);
+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);
+return;
+
+%parameterize model
+disp('Parameterizing');
+md=geography(md,groundingline,icestreams);
+md=parameterize(md,parameterfile);
+
+%md.vx_obs=DataInterp(x_m,y_m,vx*md.yts,md.x,md.y);
+%md.vy_obs=DataInterp(x_m,y_m,vy*md.yts,md.x,md.y);
+%md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
+%return;
Index: /issm/trunk/src/m/utils/Mesh/YamsCall.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1327)
+++ /issm/trunk/src/m/utils/Mesh/YamsCall.m	(revision 1327)
@@ -0,0 +1,168 @@
+function md=YamsCall(md,field,hmin,hmax,gradation,epsilon),
+%YAMSCALL - call yams
+%
+%   build a metric using the Hessian of the given field
+%   call Yams and the output mesh is plugged onto the model
+%   -hmin = minimum edge length (m)
+%   -hmax = maximum edge length (m)
+%   -gradation = maximum edge length gradation between 2 elements
+%   -epsilon = average error on each element (m/yr)
+%
+%   Usage:
+%      md=YamsCall(md,field,hmin,hmax,epsilon)
+%
+%   Example:
+%      md=YamsCall(md,md.vel_obs,1500,10^8,1.3,0.9);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build metric %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+global ISSM_DIR
+
+%initialization
+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; 
+
+%build some usefull variables
+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.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=[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(:,1),3,1),numberofgrids,1)./weights ...
+	sparse(line,ones(linesize,1),repmat(areas.*hessian(:,2),3,1),numberofgrids,1)./weights ...
+	sparse(line,ones(linesize,1),repmat(areas.*hessian(:,3),3,1),numberofgrids,1)./weights ];
+
+test=1;
+
+tic
+	if test==1,
+		%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/hmin^2 0 1/hmin^2],length(pos1),1);
+		metric(pos2,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos2),1);
+
+	end
+	if test==2,
+		for i=1:md.numberofgrids
+			H=[hessian(i,1) hessian(i,2)
+			hessian(i,2) hessian(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(i,:)=[metricTria(1,1) metricTria(1,2) metricTria(2,2)];
+		end
+	end
+toc
+
+if any(isnan(metric)),
+	error('YamsCall error message: NaN in the metric...')
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build FILES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Tria=md.elements;
+Coor=[md.x md.y];
+save -ascii carre0.tria Tria
+save -ascii carre0.coor Coor
+save -ascii carre0.met  metric
+
+fid=fopen('carre0.mesh','w');
+
+%initialiation
+fprintf(fid,'\n%s\n%i\n','MeshVersionFormatted',1);
+
+%dimension
+fprintf(fid,'\n%s\n%i\n','Dimension',2);
+
+%Vertices
+fprintf(fid,'\n%s\n%i\n','Vertices',md.numberofgrids);
+for i=1:md.numberofgrids,
+	fprintf(fid,'\n%f %f %i',md.x(i),md.y(i),0);
+end
+
+%Triangles
+fprintf(fid,'\n\n%s\n%i\n','Triangles',md.numberofelements);
+for i=1:md.numberofelements,
+	fprintf(fid,'\n%i %i %i %i',md.elements(i,1),md.elements(i,2),md.elements(i,3),0);
+end
+
+%Corners
+pos=find(md.gridonboundary);
+fprintf(fid,'\n\n%s\n%i\n','Corners',length(pos));
+for i=1:length(pos),
+	fprintf(fid,'\n%i',pos(i));
+end
+
+%close
+fclose(fid);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% call Yams %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%call yams
+addpath([ISSM_DIR '/externalpackages/yams']);
+system(['yams2-linux -O 1 -v 10 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+rmpath([ISSM_DIR '/externalpackages/yams']);
+
+%plug new mesh
+Tria=load('carre1.tria');
+Coor=load('carre1.coor');
+md.x=Coor(:,1);
+md.y=Coor(:,2);
+md.z=zeros(size(Coor,1),1);
+md.elements=Tria;
+
+md.numberofgrids=size(Coor,1);
+md.numberofelements=size(Tria,1);