Index: /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m (revision 1328) +++ /issm/trunk/src/m/utils/Mesh/BuildAdaptedMeshYams.m (revision 1329) @@ -19,9 +19,9 @@ %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'; %%} Index: /issm/trunk/src/m/utils/Mesh/YamsCall.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/YamsCall.m (revision 1328) +++ /issm/trunk/src/m/utils/Mesh/YamsCall.m (revision 1329) @@ -19,4 +19,5 @@ global ISSM_DIR +t1=clock; fprintf('%s',' computing Hessian...'); %initialization index=md.elements; @@ -59,63 +60,43 @@ 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,33 +110,31 @@ %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,2 +146,5 @@ md.numberofgrids=size(Coor,1); md.numberofelements=size(Tria,1); + +%clean up: +system('rm carre0.mesh carre0.met carre1.tria carre1.coor carre1.meshb');