Changeset 1329

Timestamp:

07/16/09 08:38:55 (16 years ago)

Author:

Mathieu Morlighem

Message:

speeded up Yams call (no more loop in fprintf... we should do the same for expwrite vim YamsCall.m )

Location:

issm/trunk/src/m/utils/Mesh

Files:

• BuildAdaptedMeshYams.m (modified) (1 diff)
• YamsCall.m (modified) (4 diffs)

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

@@ -19 +19 @@
 %Antarctica settings
 nsteps=6;
-resolution=10000;
-hmin=500;
+resolution=8000;
+hmin=50;
 hmax=3*10^5;
 gradation=[1.5*ones(2,1);3*ones(nsteps-2,1)];
-epsilon=5*10^-1;
+epsilon=2*10^-0;
 domainoutline='EnvDomainOutline200km.exp';
 %%}

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

@@ -19 +19 @@
 global ISSM_DIR
 
+t1=clock; fprintf('%s','      computing Hessian...');
 %initialization
 index=md.elements;
@@ -59 +60 @@
         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 ];
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
-test=1;
+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);
 
-tic
-        if test==1,
-                %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);
 
-                %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 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)]);
 
-                %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/hmin^2 0 1/hmin^2],length(pos1),1);
-                metric(pos2,:)=repmat([1/hmin^2 0 1/hmin^2],length(pos2),1);
-
-        end
-        if test==2,
-                for i=1:md.numberofgrids
-                        H=[hessian(i,1) hessian(i,2)
-                        hessian(i,2) hessian(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(i,:)=[metricTria(1,1) metricTria(1,2) metricTria(2,2)];
-                end
-        end
-toc
+%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);
 
 if any(isnan(metric)),
         error('YamsCall error message: NaN in the metric...')
 end
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Build FILES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-Tria=md.elements;
-Coor=[md.x md.y];
-save -ascii carre0.tria Tria
-save -ascii carre0.coor Coor
+t1=clock; fprintf('%s','      writing initial mesh files...');
 save -ascii carre0.met  metric
 
@@ -129 +110 @@
 
 %Vertices
-fprintf(fid,'\n%s\n%i\n','Vertices',md.numberofgrids);
-for i=1:md.numberofgrids,
-        fprintf(fid,'\n%f %f %i',md.x(i),md.y(i),0);
-end
+fprintf(fid,'\n%s\n%i\n\n','Vertices',md.numberofgrids);
+fprintf(fid,'%8g %8g %i\n',[md.x md.y zeros(md.numberofgrids,1)]');
 
 %Triangles
-fprintf(fid,'\n\n%s\n%i\n','Triangles',md.numberofelements);
-for i=1:md.numberofelements,
-        fprintf(fid,'\n%i %i %i %i',md.elements(i,1),md.elements(i,2),md.elements(i,3),0);
-end
+fprintf(fid,'\n\n%s\n%i\n\n','Triangles',md.numberofelements);
+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','Corners',length(pos));
-for i=1:length(pos),
-        fprintf(fid,'\n%i',pos(i));
-end
+fprintf(fid,'\n\n%s\n%i\n\n','Corners',length(pos));
+fprintf(fid,'%i\n',pos');
 
 %close
 fclose(fid);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% call Yams %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %call yams
+fprintf('%s\n','      call Yams...');
 addpath([ISSM_DIR '/externalpackages/yams']);
-system(['yams2-linux -O 1 -v 10 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
+system(['yams2-linux -O 1 -v 1 -ecp -hgrad ' num2str(gradation)  ' carre0 carre1']);
 rmpath([ISSM_DIR '/externalpackages/yams']);
+t2=clock;fprintf('%s\n',[' done (' num2str(etime(t2,t1)) ' seconds)']);
 
 %plug new mesh
+t1=clock; fprintf('\n%s','      reading final mesh files...');
 Tria=load('carre1.tria');
 Coor=load('carre1.coor');
@@ -167 +146 @@
 md.numberofgrids=size(Coor,1);
 md.numberofelements=size(Tria,1);
+
+%clean up:
+system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');