function md=meshyams(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 % optional arguments: 'riftoutline': if rifts are present, specifies rift outline file. % % % 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=getfieldvalue(options,'domainoutline'); disp(sprintf(' %-15s: ''%s''','DomainOutline',domainoutline)); riftoutline=getfieldvalue(options,'riftoutline','N/A'); disp(sprintf(' %-15s: ''%s''','riftoutline',riftoutline)); groundeddomain=getfieldvalue(options,'groundeddomain','N/A'); disp(sprintf(' %-15s: ''%s''','GroundedDomain',groundeddomain)); velocities=getfieldvalue(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...'); if strcmpi(riftoutline,'N/A'); md=mesh(md,domainoutline,resolution); else md=mesh(md,domainoutline,riftoutline,resolution); md=meshprocessrifts(md,domainoutline); end 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,groundeddomain,'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=YamsCall(md,field,hmin,hmax,gradation(i),epsilon); %if we have rifts, we just messed them up, we need to recreate the segments that constitute those %rifts, because the segments are used in YamsCall to freeze the rifts elements during refinement. if md.numrifts, md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofgrids); md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity); md.segments=findsegments(md); md=meshyamsrecreateriftsegments(md); end 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,groundeddomain,'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); %deal with rifts if md.numrifts, %first, recreate rift segments md=meshyamsrecreateriftsegments(md); %using the segments, recreate the penaltypairs for j=1:md.numrifts, rift=md.rifts(j); %build normals and lengths of segments: lengths=sqrt((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))).^2 + (md.y(rift.segments(:,1))-md.y(rift.segments(:,2))).^2 ); normalsx=cos(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1))))); normalsy=sin(atan2((md.x(rift.segments(:,1))-md.x(rift.segments(:,2))) , (md.y(rift.segments(:,2))-md.y(rift.segments(:,1))))); %ok, build penaltypairs: numpenaltypairs=length(rift.segments)/2-1; rift.penaltypairs=zeros(numpenaltypairs,7); for i=1:numpenaltypairs, rift.penaltypairs(i,1)=rift.segments(i,2); rift.penaltypairs(i,2)=rift.segments(end-i,2); rift.penaltypairs(i,3)=rift.segments(i,3); rift.penaltypairs(i,4)=rift.segments(end-i,3); rift.penaltypairs(i,5)=normalsx(i)+normalsx(i+1); rift.penaltypairs(i,6)=normalsy(i)+normalsy(i+1); rift.penaltypairs(i,7)=(lengths(i)+lengths(i+1))/2; end %renormalize norms: norms=sqrt(rift.penaltypairs(:,5).^2+rift.penaltypairs(:,6).^2); rift.penaltypairs(:,5)=rift.penaltypairs(:,5)./norms; rift.penaltypairs(:,6)=rift.penaltypairs(:,6)./norms; md.rifts(j)=rift; end end