Changeset 3199

Timestamp:

03/05/10 13:44:55 (15 years ago)

Author:

Eric.Larour

Message:

Added second thickness metric to bamg wrapper.

File:

• issm/trunk/src/m/classes/public/mesh/meshbamg.m (modified) (6 diffs)

issm/trunk/src/m/classes/public/mesh/meshbamg.m

@@ -6 +6 @@
 %      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]
+%      arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
+%      arguments can be: 'thickness': matlab file containing the thickness [m]
 %      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: 'epsilon': average interpolation error wished [m/yr [m]]
 %      optional arguments: 'hmin': minimum edge length
 %      optional arguments: 'hmanx': maximum edge
@@ -20 +21 @@
 %      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);
+%      md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','thickness','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
 
 %recover options
@@ -26 +28 @@
 
 %recover some fields
-disp('MeshYams Options:')
+disp('MeshBamg Options:')
 domainoutline=getfieldvalueerr(options,'domainoutline');
 disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
@@ -33 +35 @@
 velocities=getfieldvalueerr(options,'velocities');
 disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+thickness=getfieldvalue(options,'thickness','none');
+disp(sprintf('   %-15s: ''%s''','Thickness',thickness));
 resolution=getfieldvalue(options,'resolution',5000);
 disp(sprintf('   %-15s: %f','Resolution',resolution));
@@ -53 +57 @@
 %load velocities 
 disp('loading velocities...');
-Names=VelFindVarNames(velocities);
+NamesV=VelFindVarNames(velocities);
 Vel=load(velocities);
+
+%compute thickness?
+if ~strcmpi(thickness,'none'),
+        disp('loading thickness...');
+        NamesH=FieldFindVarNames(thickness);
+        Thi=load(thickness);
+        thicknesspresent=1;
+else
+        thicknesspresent=0;
+end
 
 %start mesh adaptation
@@ -62 +76 @@
         %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);
+        if strcmpi(NamesV.interp,'grid'),
+                vx_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.x,md.y,0);
+                vy_obs=InterpFromGridToMesh(Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.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);
+                vx_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vxname),md.x,md.y,0);
+                vy_obs=InterpFromMeshToMesh2d(Vel.(NamesV.indexname),Vel.(NamesV.xname),Vel.(NamesV.yname),Vel.(NamesV.vyname),md.x,md.y,0);
         end
         field=sqrt(vx_obs.^2+vy_obs.^2);
-
+        
+        if thicknesspresent,
+                disp('   interpolating thickness...');
+                if strcmpi(NamesH.interp,'grid'),
+                        h=InterpFromGridToMesh(Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.x,md.y,0);
+                else
+                        h=InterpFromMeshToMesh2d(Thi.(NamesH.indexname),Thi.(NamesH.xname),Thi.(NamesH.yname),Thi.(NamesH.dataname),md.x,md.y,0);
+                end
+                field=[field h];
+        end
+        
         %set gridonwater field
         if ~strcmp(groundeddomain,'N/A'),