source: issm/trunk/src/m/model/mesh/meshyams.m@ 8298

Last change on this file since 8298 was 8298, checked in by seroussi, 14 years ago

changed grid to node in matlab

File size: 7.1 KB
RevLine 
[1338]1function md=meshyams(md,varargin);
[1331]2%MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
[1327]3%
4% Usage:
[1338]5% md=meshyams(md,varargin);
6% where varargin is a lit of paired arguments.
7% arguments can be: 'domainoutline': Argus file containing the outline of the domain to be meshed
8% arguments can be: 'velocities': matlab file containing the velocities [m/yr]
9% optional arguments: 'groundeddomain': Argus file containing the outline of the grounded ice
10% this option is used to minimize the metric on water (no refinement)
11% optional arguments: 'resolution': initial mesh resolution [m]
12% optional arguments: 'nsteps': number of steps of mesh adaptation
13% optional arguments: 'epsilon': average interpolation error wished [m/yr]
14% optional arguments: 'hmin': minimum edge length
15% optional arguments: 'hmanx': maximum edge
[2749]16% optional arguments: 'riftoutline': if rifts are present, specifies rift outline file.
[1338]17%
18%
19% Examples:
20% md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat');
21% md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp');
22% md=meshyams(md,'domainoutline','Domain.exp','velocities','vel.mat','groundeddomain','ground.exp','nsteps',6,'epsilon',2,'hmin',500,'hmax',30000);
[1327]23
[1338]24%recover options
[2396]25options=pairoptions(varargin{:});
26options=deleteduplicates(options,1);
[1327]27
[1338]28%recover some fields
[1346]29disp('MeshYams Options:')
[5088]30domainoutline=getfieldvalue(options,'domainoutline');
[2396]31disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline));
[2688]32riftoutline=getfieldvalue(options,'riftoutline','N/A');
33disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline));
[2396]34groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
35disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain));
[5088]36velocities=getfieldvalue(options,'velocities');
[2396]37disp(sprintf(' %-15s: ''%s''','Velocities',velocities));
38resolution=getfieldvalue(options,'resolution',5000);
39disp(sprintf(' %-15s: %f','Resolution',resolution));
40nsteps=getfieldvalue(options,'nsteps',6);
41disp(sprintf(' %-15s: %i','nsteps',nsteps));
42gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
43disp(sprintf(' %-15s: %g','gradation',gradation(1)));
44epsilon=getfieldvalue(options,'epsilon',3);
45disp(sprintf(' %-15s: %f','epsilon',epsilon));
46hmin=getfieldvalue(options,'hmin',500);
47disp(sprintf(' %-15s: %f','hmin',hmin));
48hmax=getfieldvalue(options,'hmax',150*10^3);
49disp(sprintf(' %-15s: %f\n','hmax',hmax));
[1336]50
[1327]51%mesh with initial resolution
[1338]52disp('Initial mesh generation...');
[2688]53if strcmpi(riftoutline,'N/A');
54 md=mesh(md,domainoutline,resolution);
55else
56 md=mesh(md,domainoutline,riftoutline,resolution);
[3260]57 md=meshprocessrifts(md,domainoutline);
[2688]58end
[1346]59disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
[1327]60
61%load velocities
[1338]62disp('loading velocities...');
[2568]63Names=VelFindVarNames(velocities);
64Vel=load(velocities);
[1327]65
66%start mesh adaptation
67for i=1:nsteps,
68 disp(['Iteration #' num2str(i) '/' num2str(nsteps)]);
69
70 %interpolate velocities onto mesh
[1346]71 disp(' interpolating velocities...');
[8298]72 if strcmpi(Names.interp,'node'),
[2568]73 vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
74 vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
75 else
[2573]76 vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
77 vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
[2568]78 end
[1339]79 field=sqrt(vx_obs.^2+vy_obs.^2);
[1327]80
[8298]81 %set nodeonwater field
[1338]82 if ~strcmp(groundeddomain,'N/A'),
[8298]83 nodeground=ContourToMesh(md.elements,md.x,md.y,groundeddomain,'node',2);
84 md.nodeonwater=ones(md.numberofnodes,1);
85 md.nodeonwater(find(nodeground))=0;
[1336]86 else
[8298]87 md.nodeonwater=zeros(md.numberofnodes,1);
[1336]88 end
89
[1327]90 %adapt according to velocities
[1346]91 disp(' adapting...');
[1338]92 md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
[1346]93
[2737]94 %if we have rifts, we just messed them up, we need to recreate the segments that constitute those
95 %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement.
96 if md.numrifts,
[8298]97 md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
[2737]98 md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
99 md.segments=findsegments(md);
100 md=meshyamsrecreateriftsegments(md);
101 end
102
[1327]103end
104
105disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
[1336]106
[2078]107%Now, build the connectivity tables for this mesh.
[8298]108md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
[2078]109md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
110
[1346]111%recreate segments
112md.segments=findsegments(md);
[8298]113md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;
[1336]114
115%Fill in rest of fields:
[8298]116md.z=zeros(md.numberofnodes,1);
117md.nodeonbed=ones(md.numberofnodes,1);
118md.nodeonsurface=ones(md.numberofnodes,1);
[1336]119md.elementonbed=ones(md.numberofelements,1);
120md.elementonsurface=ones(md.numberofelements,1);
[1341]121if ~strcmp(groundeddomain,'N/A'),
[8298]122 nodeground=ContourToMesh(md.elements,md.x,md.y,groundeddomain,'node',2);
123 md.nodeonwater=ones(md.numberofnodes,1);
124 md.nodeonwater(find(nodeground))=0;
[1341]125else
[8298]126 md.nodeonwater=zeros(md.numberofnodes,1);
[1341]127end
[8298]128if strcmpi(Names.interp,'node'),
[2574]129 md.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
130 md.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
131else
132 md.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
133 md.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
134end
[1341]135md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
[2737]136
137%deal with rifts
138if md.numrifts,
139 %first, recreate rift segments
140 md=meshyamsrecreateriftsegments(md);
141
142 %using the segments, recreate the penaltypairs
143 for j=1:md.numrifts,
144 rift=md.rifts(j);
145
146 %build normals and lengths of segments:
147 lengths=sqrt((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 );
148 normalsx=cos(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
149 normalsy=sin(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1)))));
150
151 %ok, build penaltypairs:
152 numpenaltypairs=length(rift.segments)/2-1;
153 rift.penaltypairs=zeros(numpenaltypairs,7);
154
155 for i=1:numpenaltypairs,
156 rift.penaltypairs(i,1)=rift.segments(i,2);
157 rift.penaltypairs(i,2)=rift.segments(end-i,2);
158 rift.penaltypairs(i,3)=rift.segments(i,3);
159 rift.penaltypairs(i,4)=rift.segments(end-i,3);
160 rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1);
161 rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1);
162 rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2;
163 end
164 %renormalize norms:
165 norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2);
166 rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms;
167 rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms;
168
169 md.rifts(j)=rift;
170 end
171
172end
Note: See TracBrowser for help on using the repository browser.