Changeset 1336

Timestamp:

07/16/09 12:41:01 (16 years ago)

Author:

Mathieu Morlighem

Message:

Added segments + metric tweaking on water

Location:

issm/trunk/src/m/utils/Mesh

Files:

• MeshYams.m (modified) (4 diffs)
• YamsCall.m (modified) (2 diffs)

issm/trunk/src/m/utils/Mesh/MeshYams.m

@@ -1 @@
-function md=MeshYams(md,velpath,domainoutline);
+function md=MeshYams(md,velpath,domainoutline,varargin);
 %MESHYAMS - Build model of Antarctica by refining according to observed velocity error estimator
 %
@@ -23 +23 @@
 hmax=150*10^3;   %150km
 gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=3*10^-0; %3m/a interpolation error
+epsilon=2*10^-0; %3m/a interpolation error
 scale=1;
 %%}
+
+%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
 
 %mesh with initial resolution
@@ -45 +57 @@
         md.vel_obs=averaging(md,sqrt(md.vx_obs.^2+md.vy_obs.^2),2);
 
+        %set gridonwater field
+        if waterflag,
+                gridground=ContourToMesh(md.elements,md.x,md.y,expread(groundedoutline,1),'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');
@@ -51 +72 @@
         
 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 the segments
+elementconnectivity=md.elementconnectivity;
+elementonboundary=double(elementconnectivity(:,end)==0);
+pos=find(elementonboundary);
+num_segments=length(pos);
+segments=zeros(num_segments,3);
+for i=1:num_segments,
+        el1=pos(i);
+        els2=elementconnectivity(el1,find(elementconnectivity(el1,:)));
+        flag=intersect(md.elements(els2(1),:),md.elements(els2(2),:));
+        nods1=md.elements(el1,:);
+        nods1(find(nods1==flag))=[];
+        segments(i,:)=[nods1 el1];
+
+        ord1=find(nods1(1)==md.elements(el1,:));
+        ord2=find(nods1(2)==md.elements(el1,:));
+
+        %swap segment grids if necessary
+        if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
+                temp=segments(i,1);
+                segments(i,1)=segments(i,2);
+                segments(i,2)=temp;
+        end
+        segments(i,1:2)=fliplr(segments(i,1:2));
+end
+md.segments=segments;
+
+%Fill in rest of fields:
+md.z=zeros(md.numberofgrids,1);
+md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=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);

issm/trunk/src/m/utils/Mesh/YamsCall.m

@@ -88 +88 @@
 
 %some corrections for 0 eigen values
-metric(pos1,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos1),1);
-metric(pos2,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos2),1);
+metric(pos1,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos1),1);
+metric(pos2,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos2),1);
+
+%take care of water elements
+if length(md.gridonwater)==numberofgrids;
+        pos=find(md.gridonwater);
+        metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+end
 
 if any(isnan(metric)),
@@ -143 +149 @@
 md.numberofgrids=size(Coor,1);
 md.numberofelements=size(Tria,1);
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+t2=clock;fprintf('%s\n\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
 %clean up: