[1338] | 1 | function 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] | 25 | options=pairoptions(varargin{:});
|
---|
| 26 | options=deleteduplicates(options,1);
|
---|
[1327] | 27 |
|
---|
[1338] | 28 | %recover some fields
|
---|
[1346] | 29 | disp('MeshYams Options:')
|
---|
[5088] | 30 | domainoutline=getfieldvalue(options,'domainoutline');
|
---|
[2396] | 31 | disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline));
|
---|
[2688] | 32 | riftoutline=getfieldvalue(options,'riftoutline','N/A');
|
---|
| 33 | disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline));
|
---|
[2396] | 34 | groundeddomain=getfieldvalue(options,'groundeddomain','N/A');
|
---|
| 35 | disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain));
|
---|
[5088] | 36 | velocities=getfieldvalue(options,'velocities');
|
---|
[2396] | 37 | disp(sprintf(' %-15s: ''%s''','Velocities',velocities));
|
---|
| 38 | resolution=getfieldvalue(options,'resolution',5000);
|
---|
| 39 | disp(sprintf(' %-15s: %f','Resolution',resolution));
|
---|
| 40 | nsteps=getfieldvalue(options,'nsteps',6);
|
---|
| 41 | disp(sprintf(' %-15s: %i','nsteps',nsteps));
|
---|
| 42 | gradation=getfieldvalue(options,'gradation',2*ones(nsteps,1));
|
---|
| 43 | disp(sprintf(' %-15s: %g','gradation',gradation(1)));
|
---|
| 44 | epsilon=getfieldvalue(options,'epsilon',3);
|
---|
| 45 | disp(sprintf(' %-15s: %f','epsilon',epsilon));
|
---|
| 46 | hmin=getfieldvalue(options,'hmin',500);
|
---|
| 47 | disp(sprintf(' %-15s: %f','hmin',hmin));
|
---|
| 48 | hmax=getfieldvalue(options,'hmax',150*10^3);
|
---|
| 49 | disp(sprintf(' %-15s: %f\n','hmax',hmax));
|
---|
[1336] | 50 |
|
---|
[1327] | 51 | %mesh with initial resolution
|
---|
[1338] | 52 | disp('Initial mesh generation...');
|
---|
[2688] | 53 | if strcmpi(riftoutline,'N/A');
|
---|
| 54 | md=mesh(md,domainoutline,resolution);
|
---|
| 55 | else
|
---|
| 56 | md=mesh(md,domainoutline,riftoutline,resolution);
|
---|
[3260] | 57 | md=meshprocessrifts(md,domainoutline);
|
---|
[2688] | 58 | end
|
---|
[1346] | 59 | disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
|
---|
[1327] | 60 |
|
---|
| 61 | %load velocities
|
---|
[1338] | 62 | disp('loading velocities...');
|
---|
[2568] | 63 | Names=VelFindVarNames(velocities);
|
---|
| 64 | Vel=load(velocities);
|
---|
[1327] | 65 |
|
---|
| 66 | %start mesh adaptation
|
---|
| 67 | for 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.
|
---|
[9619] | 96 | if md.rifts.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] | 103 | end
|
---|
| 104 |
|
---|
| 105 | disp(['Final mesh, number of elements: ' num2str(md.numberofelements)]);
|
---|
[1336] | 106 |
|
---|
[2078] | 107 | %Now, build the connectivity tables for this mesh.
|
---|
[8298] | 108 | md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
|
---|
[2078] | 109 | md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
|
---|
| 110 |
|
---|
[1346] | 111 | %recreate segments
|
---|
| 112 | md.segments=findsegments(md);
|
---|
[8298] | 113 | md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;
|
---|
[1336] | 114 |
|
---|
| 115 | %Fill in rest of fields:
|
---|
[8298] | 116 | md.z=zeros(md.numberofnodes,1);
|
---|
| 117 | md.nodeonbed=ones(md.numberofnodes,1);
|
---|
| 118 | md.nodeonsurface=ones(md.numberofnodes,1);
|
---|
[1336] | 119 | md.elementonbed=ones(md.numberofelements,1);
|
---|
| 120 | md.elementonsurface=ones(md.numberofelements,1);
|
---|
[1341] | 121 | if ~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] | 125 | else
|
---|
[8298] | 126 | md.nodeonwater=zeros(md.numberofnodes,1);
|
---|
[1341] | 127 | end
|
---|
[8298] | 128 | if strcmpi(Names.interp,'node'),
|
---|
[9681] | 129 | md.inversion.vx_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
|
---|
| 130 | md.inversion.vy_obs=InterpFromGridToMesh(Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
|
---|
[2574] | 131 | else
|
---|
[9681] | 132 | md.inversion.vx_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vxname),md.x,md.y,0);
|
---|
| 133 | md.inversion.vy_obs=InterpFromMeshToMesh2d(Vel.(Names.indexname),Vel.(Names.xname),Vel.(Names.yname),Vel.(Names.vyname),md.x,md.y,0);
|
---|
[2574] | 134 | end
|
---|
[9681] | 135 | md.inversion.vel_obs=sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
|
---|
[2737] | 136 |
|
---|
| 137 | %deal with rifts
|
---|
[9619] | 138 | if md.rifts.numrifts,
|
---|
[2737] | 139 | %first, recreate rift segments
|
---|
| 140 | md=meshyamsrecreateriftsegments(md);
|
---|
| 141 |
|
---|
| 142 | %using the segments, recreate the penaltypairs
|
---|
[9619] | 143 | for j=1:md.rifts.numrifts,
|
---|
| 144 | rift=md.rifts.riftstruct(j);
|
---|
[2737] | 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 |
|
---|
[9619] | 169 | md.rifts.riftstruct(j)=rift;
|
---|
[2737] | 170 | end
|
---|
| 171 |
|
---|
| 172 | end
|
---|