Index: /issm/trunk/src/m/utils/Mesh/BamgCall.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/BamgCall.m	(revision 2614)
+++ /issm/trunk/src/m/utils/Mesh/BamgCall.m	(revision 2614)
@@ -0,0 +1,86 @@
+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,gradation,epsilon);
+%
+%   Example:
+%      md=YamsCall(md,md.vel_obs,1500,10^8,1.3,0.9);
+
+global ISSM_DIR
+
+%2d geometric parameter (do not change)
+scale=2/9; 
+
+%Compute Hessian
+t1=clock; fprintf('%s','      computing Hessian...');
+hessian=ComputeHessian(md.elements,md.x,md.y,field,'node');
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%Compute metric
+t1=clock; fprintf('%s','      computing metric...');
+if length(md.gridonwater)==md.numberofgrids,
+	pos=find(md.gridonwater);
+else
+	pos=[];
+end
+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%write files
+t1=clock; fprintf('%s','      writing initial mesh files...');
+fid=fopen('carre0.met','w');
+fprintf(fid,'%i %i\n',md.numberofgrids,3);
+fprintf(fid,'%i %i %i\n',metric');
+fclose(fid);
+
+fid=fopen('carre0.mesh','w');
+
+%initialiation
+fprintf(fid,'%s %i\n','MeshVersionFormatted',0);
+
+%dimension
+fprintf(fid,'\n%s\n%i\n','Dimension',2);
+
+%Vertices
+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.numberofgrids);
+fprintf(fid,'%8g %8g %i\n',[md.x md.y ones(md.numberofgrids,1)]');
+
+%Triangles
+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.numberofelements);
+fprintf(fid,'%i %i %i %i\n',[md.elements ones(md.numberofelements,1)]');
+numberofelements1=md.numberofelements;
+
+%close
+fclose(fid);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%call yams
+fprintf('%s\n','      call Bamg...');
+system(['bamg -ratio ' num2str(gradation) ' -splitpbedge -nbv 1000000 -M carre0.met -b carre0.mesh -o carre1.mesh']);
+
+%plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
+A=meshread('carre1.mesh');
+md.x=A.x;
+md.y=A.y;
+md.z=zeros(A.nods,1);
+md.elements=A.index;
+md.numberofgrids=A.nods;
+md.numberofelements=A.nels;
+numberofelements2=md.numberofelements;
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
+%display number of elements
+fprintf('\n%s %i','      inital number of elements:',numberofelements1);
+fprintf('\n%s %i\n\n','      new    number of elements:',numberofelements2);
+
+%clean up:
+system('rm carre0.mesh carre0.met carre1.mesh');
Index: /issm/trunk/src/m/utils/Mesh/meshread.m
===================================================================
--- /issm/trunk/src/m/utils/Mesh/meshread.m	(revision 2614)
+++ /issm/trunk/src/m/utils/Mesh/meshread.m	(revision 2614)
@@ -0,0 +1,37 @@
+function Struct=meshread(filename);
+
+%some checks
+if ~exist(filename),
+	error(['meshread error message: file ' filename ' not found!']);
+end
+
+fid=fopen(filename,'r');
+
+while (~feof(fid)),
+
+	A=fscanf(fid,'%s',1);
+
+	if strcmp(A,'MeshVersionFormatted');
+		Struct.Version=fscanf(fid,'%s',1);
+
+	elseif strcmp(A,'Dimension'),
+		Struct.Dimension=fscanf(fid,'%i',1);
+
+	elseif strcmp(A,'Vertices'),
+		Struct.nods=fscanf(fid,'%i',1);
+		A=fscanf(fid,'%f %f %f',[3 Struct.nods]);
+		Struct.x=A(1,:)';
+		Struct.y=A(2,:)';
+
+	elseif strcmp(A,'Triangles'),
+		Struct.nels=fscanf(fid,'%i',1);
+		A=fscanf(fid,'%i %i %i',[4 Struct.nels]);
+		Struct.index=A(1:3,:)';
+
+	else
+		%do nothing
+
+	end
+end
+
+fclose(fid);