Changeset 1346

Timestamp:

07/17/09 11:59:35 (16 years ago)

Author:

Mathieu Morlighem

Message:

improved yams 2...

Location:

issm/trunk/src/m

Files:

• classes/public/mesh/findsegments.m (added)
• classes/public/mesh/meshyams.m (modified) (3 diffs)
• utils/Mesh/ComputeMetric.m (modified) (3 diffs)
• utils/Mesh/YamsCall.m (modified) (3 diffs)

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

@@ -28 +28 @@
 
 %recover some fields
-domainoutline=yamsoptions.domainoutline;
-groundeddomain=yamsoptions.groundeddomain;
-velocities=yamsoptions.velocities;
-resolution=yamsoptions.resolution;
+disp('MeshYams Options:')
+domainoutline=yamsoptions.domainoutline;   disp(sprintf('   %-15s: ''%s''','DomainOutline',domainoutline));
+groundeddomain=yamsoptions.groundeddomain; disp(sprintf('   %-15s: ''%s''','GroundedDomain',groundeddomain));
+velocities=yamsoptions.velocities;         disp(sprintf('   %-15s: ''%s''','Velocities',velocities));
+resolution=yamsoptions.resolution;         disp(sprintf('   %-15s: %f','Resolution',resolution));
 gradation=yamsoptions.gradation;
-nsteps=yamsoptions.nsteps;
-epsilon=yamsoptions.epsilon;
-hmin=yamsoptions.hmin;
-hmax=yamsoptions.hmax;
+nsteps=yamsoptions.nsteps;                 disp(sprintf('   %-15s: %i','nsteps',nsteps));
+epsilon=yamsoptions.epsilon;               disp(sprintf('   %-15s: %f','epsilon',epsilon));
+hmin=yamsoptions.hmin;                     disp(sprintf('   %-15s: %f','hmin',hmin));
+hmax=yamsoptions.hmax;                     disp(sprintf('   %-15s: %f\n','hmax',hmax));
 
 %mesh with initial resolution
 disp('Initial mesh generation...');
 md=mesh(md,domainoutline,resolution);
+disp(['Initial mesh, number of elements: ' num2str(md.numberofelements)]);
 
 %load velocities 
@@ -51 +53 @@
 
         %interpolate velocities onto mesh
-        disp('   interpolating velocities');
+        disp('   interpolating velocities...');
         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);
@@ -66 +68 @@
 
         %adapt according to velocities
-        disp('   adapting');
+        disp('   adapting...');
         md=YamsCall(md,field,hmin,hmax,gradation(i),epsilon);
+
 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 the segments
-disp('Creating 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;
+%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.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
 md.gridonbed=ones(md.numberofgrids,1);
 md.gridonsurface=ones(md.numberofgrids,1);

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

@@ -15 +15 @@
 pos1=find(lambda1==0);
 pos2=find(lambda2==0);
+pos3=find(b==0 & lambda1==lambda2);
 
 %Modify the eigen values to control the shape of the elements
@@ -28 +29 @@
 v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2;
 
+v1x(pos3)=1; v1y(pos3)=0;
+v2x(pos3)=0; v2y(pos3)=1;
+
 %Compute new metric (for each node M=V*Lambda*V^-1)
 metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ...
@@ -40 +44 @@
 metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
 
+%take care of NaNs if any (use Matlab eig in a loop)
+[pos posj]=find(isnan(metric)); clear posj;
+if ~isempty(pos),
+        fprintf(' %i %s',length(pos),'NaN found in the metric. Use Matlab routine...');
+        for i=1:length(pos)
+                H=[hessian(pos(i),1) hessian(pos(i),2)
+                hessian(pos(i),2) hessian(pos(i),3)];
+                [u,v]=eig(H);
+                lambda1=v(1,1);
+                lambda2=v(2,2);
+                v(1,1)=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
+                v(2,2)=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
+
+                metricTria=u*v*u^(-1);
+                metric(pos(i),:)=[metricTria(1,1) metricTria(1,2) metricTria(2,2)];
+        end
+end
+
 if any(isnan(metric)),
-        error('ComputeMetric error message: NaN in the metric...')
+        error('ComputeMetric error message: NaN in the metric despite our efforts...')
 end

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

@@ -10 +10 @@
 %
 %   Usage:
-%      md=YamsCall(md,field,hmin,hmax,epsilon)
+%      md=YamsCall(md,field,hmin,hmax,gradation,epsilon);
 %
 %   Example:
 %      md=YamsCall(md,md.vel_obs,1500,10^8,1.3,0.9);
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build metric %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 global ISSM_DIR
@@ -57 +55 @@
 fprintf(fid,'%i %i %i %i\n',[md.elements zeros(md.numberofelements,1)]');
 
-%Corners
-pos=find(md.gridonboundary);
-fprintf(fid,'\n\n%s\n%i\n\n','Corners',length(pos));
-fprintf(fid,'%i\n',pos');
-
 %close
 fclose(fid);
@@ -68 +61 @@
 %call yams
 fprintf('%s\n','      call Yams...');
-system(['yams2-linux -O 1 -v 2 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+system(['yams2-linux -O 1 -v -0 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
 
 %plug new mesh