Changeset 1342

Timestamp:

07/17/09 09:01:18 (16 years ago)

Author:

Mathieu Morlighem

Message:

Improved Yams and QualityMesh

Location:

issm/trunk/src/m/utils/Mesh

Files:

• ComputeHessian.m (added)
• ComputeMetric.m (added)
• MeshQuality.m (moved) (moved from issm/trunk/src/m/classes/public/mesh/meshquality.m ) (4 diffs)
• YamsCall.m (modified) (2 diffs)

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

@@ -1 @@
-function meshquality(md);
+function MeshQuality(md,epsilon,hmin,hmax);
 %MESHQUALITY - compute mesh quality
 %
 %   Usage:
-%      meshquality(md);
+%      MeshQuality(md,epsilon,hmin,hmax);
 
-t1=clock; fprintf('%s','      computing Hessian...');
-%initialization
-field=md.vel_obs;
+%Get some variables from the model
+index=md.elements;
 x=md.x;
 y=md.y;
-index=md.elements;
-numberofgrids=md.numberofgrids;
-numberofelements=md.numberofelements;
-gradx=zeros(numberofgrids,1);
-grady=zeros(numberofgrids,1);
-metric=zeros(numberofelements,3);
 
 %2d geometric parameter (do not change)
 scale=2/9; 
-hmin=500;
-hmax=10^6;
-epsilon=5;
 
-%build some usefull variables
-line=index(:);
-linesize=3*numberofelements;
+%Compute Hessian
+hessian=ComputeHessian(index,x,y,md.vel_obs,'node');
 
-%get areas and  nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
-[alpha beta]=GetNodalFunctionsCoeff(index,md.x,md.y);
-areas=GetAreas(index,md.x,md.y);
-
-%Compute gradient for each element
-grad_elx=sum(field(index).*alpha,2); 
-grad_ely=sum(field(index).*beta,2);
-
-%update weights that holds the volume of all the element holding the grid i
-weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
-
-%Compute gradient for each grid (average of the elements around)
-gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
-grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
-gradx=gradx./weights;
-grady=grady./weights;
-
-%Compute hessian for each element
-hessian_elem=[sum(gradx(index).*alpha,2) sum(grady(index).*alpha,2) sum(grady(index).*beta,2)];
-
-%Compute Hessian on the nodes (average of the elements around)
-hessian=[sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,1),3,1),numberofgrids,1)./weights ...
-        sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,2),3,1),numberofgrids,1)./weights ...
-        sparse(line,ones(linesize,1),repmat(areas.*hessian_elem(:,3),3,1),numberofgrids,1)./weights ];
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
-
-t1=clock; fprintf('%s','      computing metric...');
-%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,2); hessian(i,2)  hessian(i,3)]
-a=hessian(:,1); b=hessian(:,2); d=hessian(:,3);
-lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
-lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
-pos1=find(lambda1==0);
-pos2=find(lambda2==0);
-
-%Modify the eigen values to control the shape of the elements
-lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
-lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
-
-%compute eigen vectors
-norm1=sqrt(8*b.^2+2*(d-a).^2+2*(d-a).*sqrt((a-d).^2+4*b.^2));
-v1x=2*b./norm1;
-v1y=((d-a)+sqrt((a-d).^2+4*b.^2))./norm1;
-norm2=sqrt(8*b.^2+2*(d-a).^2-2*(d-a).*sqrt((a-d).^2+4*b.^2));
-v2x=2*b./norm2;
-v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2;
-
-%Compute new metric (for each node M=V*Lambda*V^-1)
-metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ...
-        (v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v1y.*v2y-lambda2.*v1y.*v2y) ...
-        (v1x.*v2y-v1y.*v2x).^(-1).*(-lambda1.*v2x.*v1y+lambda2.*v1x.*v2y)]);
-
-%some corrections for 0 eigen values
-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);
-if length(md.gridonwater)==numberofgrids;
+%Compute metric
+if length(md.gridonwater)==md.numberofgrids,
         pos=find(md.gridonwater);
-        metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+else
+        pos=[];
 end
-
-%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)),
-        error('YamsCall error message: NaN in the metric...')
-end
-t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
 
 %Get Areas
@@ -122 +49 @@
 
 %compute quality:
-Q=4*sqrt(3)*V./(L1+L2+L3);
-
-%display
-X=0:0.1:4;
-hist(Q,X);
-xlim([0 3]);
-title('mesh quality distribution');
+quality=4*sqrt(3)*V./(L1+L2+L3);
 
 %compute error
-%{
-%a=hessian_elem(:,1); b=hessian_elem(:,2); d=hessian_elem(:,3);
 a=hessian(:,1); b=hessian(:,2); d=hessian(:,3);
 a=a(index)*[1;1;1]/3;
@@ -141 +60 @@
 lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
 lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
-if length(md.elementonwater)==numberofelements;
+if length(md.gridonwater)==md.numberofgrids;
         pos=find(md.gridonwater);
-        lambda1(pos)=0;
-        lambda2(pos)=0;
-end
-%}
-if length(md.gridonwater)==numberofgrids;
         lambda1(pos)=0;
         lambda2(pos)=0;
@@ -159 +73 @@
 
 %display
-plotmodel(md,'data',epsilon,'title','Interpolation error','figure',2,'log',10)
-
-disp(['Average interpolation error = ' num2str(mean(epsilon))  ' m/yr']);
-disp(['Maximum interpolation error = ' num2str(max(epsilon))  ' m/yr']);
+X=0:0.1:4; hist(quality,X); xlim([0 3]); title('mesh quality distribution','FontSize',14);
+plotmodel(md,'data',epsilon,'title','Interpolation error','figure',2)
+disp(sprintf('\n%s','Mesh Quality'));
+disp(sprintf('   %s %g','Average Mesh quality:',mean(quality)));
+disp(sprintf('   %s %g','Maximum interpolation error:',max(quality)));
+disp(sprintf('\n%s','Interpolation Error'));
+disp(sprintf('   %s %g %s','Average interpolation error:',mean(epsilon),'m/yr'));
+disp(sprintf('   %s %g %s','Maximum interpolation error:',max(epsilon),'m/yr'));

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

@@ -19 +19 @@
 global ISSM_DIR
 
-t1=clock; fprintf('%s','      computing Hessian...');
-%initialization
-index=md.elements;
-numberofgrids=md.numberofgrids;
-numberofelements=md.numberofelements;
-gradx=zeros(numberofgrids,1);
-grady=zeros(numberofgrids,1);
-metric=zeros(numberofelements,3);
-
 %2d geometric parameter (do not change)
 scale=2/9; 
 
-%build some usefull variables
-line=index(:);
-summation=1/3*ones(3,1);
-linesize=3*numberofelements;
-
-%get areas and  nodal functions coefficients N(x,y)=alpha x + beta y + gamma 
-[alpha beta]=GetNodalFunctionsCoeff(index,md.x,md.y);
-areas=GetAreas(index,md.x,md.y);
-
-%Compute gradient for each element
-grad_elx=sum(field(index).*alpha,2); 
-grad_ely=sum(field(index).*beta,2);
-
-%update weights that holds the volume of all the element holding the grid i
-weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
-
-%Compute gradient for each grid (average of the elements around)
-gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
-grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
-gradx=gradx./weights;
-grady=grady./weights;
-
-%Compute hessian for each element
-hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*alpha,2) sum(grady(index).*beta,2)];
-
-%Compute Hessian on the nodes (average of the elements around)
-hessian=[sparse(line,ones(linesize,1),repmat(areas.*hessian(:,1),3,1),numberofgrids,1)./weights ...
-        sparse(line,ones(linesize,1),repmat(areas.*hessian(:,2),3,1),numberofgrids,1)./weights ...
-        sparse(line,ones(linesize,1),repmat(areas.*hessian(:,3),3,1),numberofgrids,1)./weights ];
+%Compute Hessian
+t1=clock; fprintf('%s','      computing Hessian...');
+hessian=ComputeHessian(md.elements,md.x,md.y,field,'node');
 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
+%Compute metric
 t1=clock; fprintf('%s','      computing metric...');
-%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,2); hessian(i,2)  hessian(i,3)]
-a=hessian(:,1); b=hessian(:,2); d=hessian(:,3);
-lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
-lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
-pos1=find(lambda1==0);
-pos2=find(lambda2==0);
-
-%Modify the eigen values to control the shape of the elements
-lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2);
-lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2);
-
-%compute eigen vectors
-norm1=sqrt(8*b.^2+2*(d-a).^2+2*(d-a).*sqrt((a-d).^2+4*b.^2));
-v1x=2*b./norm1;
-v1y=((d-a)+sqrt((a-d).^2+4*b.^2))./norm1;
-norm2=sqrt(8*b.^2+2*(d-a).^2-2*(d-a).*sqrt((a-d).^2+4*b.^2));
-v2x=2*b./norm2;
-v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2;
-
-%Compute new metric (for each node M=V*Lambda*V^-1)
-metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ...
-        (v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v1y.*v2y-lambda2.*v1y.*v2y) ...
-        (v1x.*v2y-v1y.*v2x).^(-1).*(-lambda1.*v2x.*v1y+lambda2.*v1x.*v2y)]);
-
-%some corrections for 0 eigen values
-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;
+if length(md.gridonwater)==md.numberofgrids,
         pos=find(md.gridonwater);
-        metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1);
+else
+        pos=[];
 end
-
-if any(isnan(metric)),
-        error('YamsCall error message: NaN in the metric...')
-end
+metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos);
 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build FILES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+%write files
 t1=clock; fprintf('%s','      writing initial mesh files...');
 save -ascii carre0.met  metric
@@ -132 +66 @@
 t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% call Yams %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 %call yams
 fprintf('%s\n','      call Yams...');
-system(['yams2-linux -O 1 -v 1 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+system(['yams2-linux -O 1 -v 2 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
 
 %plug new mesh