Index: /issm/trunk/src/m/classes/public/mesh/meshbamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/mesh/meshbamg.m	(revision 2615)
+++ /issm/trunk/src/m/classes/public/mesh/meshbamg.m	(revision 2615)
@@ -0,0 +1,118 @@
+function md=meshbamg(md,varargin);
+%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
+%
+%   Usage:
+%      md=meshyams(md,varargin);
+%      where varargin is a lit of paired arguments. 
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'velocities': matlab file containing the velocities [m/yr]
+%      optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
+%                          this option is used to minimize the metric on water (no refinement)
+%      optional arguments: 'resolution': initial mesh resolution [m]
+%      optional arguments: 'nsteps': number of steps of mesh adaptation
+%      optional arguments: 'epsilon': average interpolation error wished [m/yr]
+%      optional arguments: 'hmin': minimum edge length
+%      optional arguments: 'hmanx': maximum edge
+%      
+%
+%   Examples:
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
+
+%recover options
+options=pairoptions(varargin{:});
+options=deleteduplicates(options,1);
+
+%recover some fields
+disp('MeshYams Options:')
+domainoutline=getfieldvalueerr(options,'domainoutline');
+disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
+groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
+disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
+velocities=getfieldvalueerr(options,'velocities');
+disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+resolution=getfieldvalue(options,'resolution',5000);
+disp(sprintf('   %-15s: %f','Resolution',resolution));
+nsteps=getfieldvalue(options,'nsteps',6);
+disp(sprintf('   %-15s: %i','nsteps',nsteps));
+gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
+disp(sprintf('   %-15s: %g','gradation',gradation(1)));
+epsilon=getfieldvalue(options,'epsilon',3);
+disp(sprintf('   %-15s: %f','epsilon',epsilon));
+hmin=getfieldvalue(options,'hmin',500);
+disp(sprintf('   %-15s: %f','hmin',hmin));
+hmax=getfieldvalue(options,'hmax',150*10^3);
+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 
+disp('loading velocities...');
+Names=VelFindVarNames(velocities);
+Vel=load(velocities);
+
+%start mesh adaptation
+for i=1:nsteps,
+	disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
+
+	%interpolate velocities onto mesh
+	disp('   interpolating velocities...');
+	if strcmpi(Names.interp,'grid'),
+		vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
+		vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+	else
+		vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
+		vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+	end
+	field=sqrt(vx_obs.^2+vy_obs.^2);
+
+	%set gridonwater field
+	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
+
+	%adapt according to velocities
+	disp('   adapting...');
+	md=BamgCall(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 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.gridonbed=ones(md.numberofgrids,1);
+md.gridonsurface=ones(md.numberofgrids,1);
+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
+if strcmpi(Names.interp,'grid'),
+	md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
+	md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+else
+	md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
+	md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
+end
+md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);