Changeset 1338

Timestamp:

07/16/09 14:55:38 (16 years ago)

Author:

Mathieu Morlighem

Message:

new way of calling meshyams

Location:

issm/trunk/src/m

Files:

• classes/public/BasinConstrain.m (modified) (1 diff)
• classes/public/extrude.m (modified) (1 diff)
• classes/public/mesh/meshyams.m (modified) (4 diffs)
• classes/public/mesh/process_mesh_options.m (added)
• classes/public/mesh/recover_mesh_options.m (added)
• classes/public/recover_areas.m (modified) (1 diff)
• classes/public/setelementstype.m (modified) (1 diff)
• classes/public/solve.m (modified) (1 diff)
• utils/Miscellaneous/recover_options.m (added)

issm/trunk/src/m/classes/public/BasinConstrain.m

@@ -61 +61 @@
 md.dirichletvalues_diag(grids,2)=md.vy_obs(grids);
 
-
-%make sure any grid with NaN velocity is spc'd:
-pos=find(isnan(md.vel_obs_raw));
-md.gridondirichlet_diag(pos)=1;
-%we spc to the smoothed value, so that control methods don't go berserk trying to figure out what reaction force to apply for the spc to stand.
-md.dirichletvalues_diag(pos,1)=md.vx_obs(pos); 
-md.dirichletvalues_diag(pos,2)=md.vy_obs(pos); 
-
 %make sure icefronts that are completely spc'd are taken out:
 free_segments=find(sum(md.gridondirichlet_diag(md.segmentonneumann_diag(:,1:2)),2)~=2);

issm/trunk/src/m/classes/public/extrude.m

@@ -36 +36 @@
 end
 
-if extrusionexponent<0,
+if extrusionexponent<=0,
         help extrude;
-        error('extrusionexponent must be >0');
+        error('extrusionexponent must be >=0');
 end
 

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

@@ -1 @@
-function md=meshyams(md,velpath,domainoutline,varargin);
+function md=meshyams(md,varargin);
 %MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
 %
 %   Usage:
-%      md=meshyams;
+%      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);
 
-%PIG settings
-%{
-nsteps=4;
-resolution=10000;
-hmin=1500;
-hmax=10^7;
-gradation=3*ones(nsteps,1);
-epsilon=10^-0;
-scale=1;
-%}
+%recover options
+yamsoptions=recover_options(varargin{:});
 
-%%{
-%Antarctica settings
-nsteps=6;
-resolution=5000;
-hmin=300;        %300m
-hmax=150*10^3;   %150km
-gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=2*10^-0; %3m/a interpolation error
-scale=1;
-%%}
+%add default options
+yamsoptions=process_mesh_options(yamsoptions);
 
-%check number of arguments
-waterflag=0;
-if nargin==4,
-        groundedoutline=varargin{1};
-        if exist(groundedoutline);
-                disp(['grounded ice domain found. Metric will be minimum on water']);
-                waterflag=1;
-        else
-                error(['MeshYams error message: file ' groundedoutline ' not found.']);
-        end
-end
+%recover some fields
+domainoutline=yamsoptions.domainoutline;
+groundeddomain=yamsoptions.groundeddomain;
+velocities=yamsoptions.velocities;
+resolution=yamsoptions.resolution;
+gradation=yamsoptions.gradation;
+nsteps=yamsoptions.nsteps;
+epsilon=yamsoptions.epsilon;
+hmin=yamsoptions.hmin;
+hmax=yamsoptions.hmax;
 
 %mesh with initial resolution
+disp('Initial mesh generation...');
 md=mesh(md,domainoutline,resolution);
 
 %load velocities 
-load(velpath);
-
-disp(['First mesh, number of elements: ' num2str(md.numberofelements)]);
+disp('loading velocities...');
+load(velocities);
 
 %start mesh adaptation
@@ -53 +52 @@
         %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);
+        vx_obs=InterpFromGrid(x_m,y_m,vx,md.x,md.y,0);
+        vy_obs=InterpFromGrid(x_m,y_m,vy,md.x,md.y,0);
+        field=averaging(md,sqrt(vx_obs.^2+vy_obs.^2),2);
 
         %set gridonwater field
-        if waterflag,
-                gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'node',2);
+        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;
@@ -68 +67 @@
         %adapt according to velocities
         disp('   adapting');
-        md=YamsCall(md,md.vel_obs,hmin,hmax,gradation(i),epsilon);
+        md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
 end
         
@@ -78 +77 @@
 
 %Recreate the segments
+disp('Creating segments...');
 elementconnectivity=md.elementconnectivity;
 elementonboundary=double(elementconnectivity(:,end)==0);

issm/trunk/src/m/classes/public/recover_areas.m

@@ -15 +15 @@
 
         %make sure length(varargin) is even
-        varargin=varargin{1};
         if mod(length(varargin),2),
                 error('recover_elementstype_options error message: an even number of options is necessary');

issm/trunk/src/m/classes/public/setelementstype.m

@@ -26 +26 @@
 end
 
-[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);
+[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin{:});
 
 %Flag the elements that has not been flagged as filltype

issm/trunk/src/m/classes/public/solve.m

@@ -20 +20 @@
 
 %recover options
-options=recover_solve_options(md,varargin{:});
+options=recover_options(varargin{:});
 
 %add default options