Index: sm/trunk/src/m/classes/public/mesh/meshquality.m =================================================================== --- /issm/trunk/src/m/classes/public/mesh/meshquality.m (revision 1341) +++ (revision ) @@ -1,164 +1,0 @@ -function meshquality(md); -%MESHQUALITY - compute mesh quality -% -% Usage: -% meshquality(md); - -t1=clock; fprintf('%s',' computing Hessian...'); -%initialization -field=md.vel_obs; -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; - -%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; - pos=find(md.gridonwater); - metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1); -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)']); - -%Get Areas -areas=GetAreas(index,x,y); - -%length edges vectors -e1x=[x(index(:,2))-x(index(:,1))]; -e1y=[y(index(:,2))-y(index(:,1))]; -e2x=[x(index(:,3))-x(index(:,2))]; -e2y=[y(index(:,3))-y(index(:,2))]; -e3x=[x(index(:,1))-x(index(:,3))]; -e3y=[y(index(:,1))-y(index(:,3))]; - -%metric of each the 3 grids for each element -M1=metric(index(:,1),:); -M2=metric(index(:,2),:); -M3=metric(index(:,3),:); - -%Get edge length in the metric -L1=1/2*(sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y))+sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y))); -L2=1/2*(sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y))+sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y))); -L3=1/2*(sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y))+sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y))); - -%area in the metric -V=1/3*areas.*(sqrt(M1(:,1).*M1(:,3)-M1(:,2).^2)+sqrt(M2(:,1).*M2(:,3)-M2(:,2).^2)+sqrt(M3(:,1).*M3(:,3)-M3(:,2).^2)); - -%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'); - -%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; -b=b(index)*[1;1;1]/3; -d=d(index)*[1;1;1]/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)); -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; - pos=find(md.gridonwater); - lambda1(pos)=0; - lambda2(pos)=0; -end -%} -if length(md.gridonwater)==numberofgrids; - lambda1(pos)=0; - lambda2(pos)=0; -end -lambda1=lambda1(index)*[1;1;1]/3; -lambda2=lambda2(index)*[1;1;1]/3; - -lambdamax=max(lambda1,lambda2); -hmax=max(max(sqrt(e1x.^2+e1y.^2),sqrt(e2x.^2+e2y.^2)),sqrt(e3x.^2+e3y.^2)); -epsilon=scale*hmax.^2.*lambdamax; - -%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']); Index: /issm/trunk/src/m/utils/Mesh/ComputeHessian.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/ComputeHessian.m (revision 1342) +++ /issm/trunk/src/m/utils/Mesh/ComputeHessian.m (revision 1342) @@ -0,0 +1,60 @@ +function hessian=ComputeHessian(index,x,y,field,type) +%COMPUTEHESSIAN - compute hessian matrix from a field +% +% Compute the hessian matrix of a given field +% return the three components Hxx Hxy Hyy +% for each element or each node +% +% Usage: +% hessian=ComputeHessian(index,x,y,field,type) +% +% Example: +% hessian=ComputeHessian(md.elements,md.x,md.y,md.vel_obs,'node') + +%some variables +numberofgrids=length(x); +numberofelements=size(index,1); + +%some checks +if length(field)~=numberofgrids & length(field)~=numberofelements, + error('ComputeHessian error message: the given field size not supported yet'); +end +if strcmpi(type,'node') & strcmpi(type,'element'), + error('ComputeHessian error message: only ''node'' or ''element'' type supported yet'); +end + +%initialization +line=index(:); +linesize=3*numberofelements; + +%get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma +[alpha beta]=GetNodalFunctionsCoeff(index,x,y); +areas=GetAreas(index,x,y); + +%comput 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 field on nodes if on elements +if length(field)==numberofelements, + field=sparse(line,ones(linesize,1),repmat(areas.*field,3,1),numberofgrids,1)./weights ; +end + +%Compute gradient for each element +grad_elx=sum(field(index).*alpha,2); +grad_ely=sum(field(index).*beta,2); + +%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)]; + +if strcmpi(type,'node') + %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 ]; +end Index: /issm/trunk/src/m/utils/Mesh/ComputeMetric.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/ComputeMetric.m (revision 1342) +++ /issm/trunk/src/m/utils/Mesh/ComputeMetric.m (revision 1342) @@ -0,0 +1,44 @@ +function metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos) +%COMPUTEMETRIC - compute metric from an Hessian +% +% Usage: +% metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos) +% pos is contains the positions where the metric is wished to be maximized (water?) +% +% Example: +% metric=ComputeMetric(hessian,2/9,10^-1,100,10^5,find(md.gridonwater) + +%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 +metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1); + +if any(isnan(metric)), + error('ComputeMetric error message: NaN in the metric...') +end Index: /issm/trunk/src/m/utils/Mesh/MeshQuality.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/MeshQuality.m (revision 1342) +++ /issm/trunk/src/m/utils/Mesh/MeshQuality.m (revision 1342) @@ -0,0 +1,82 @@ +function MeshQuality(md,epsilon,hmin,hmax); +%MESHQUALITY - compute mesh quality +% +% Usage: +% MeshQuality(md,epsilon,hmin,hmax); + +%Get some variables from the model +index=md.elements; +x=md.x; +y=md.y; + +%2d geometric parameter (do not change) +scale=2/9; + +%Compute Hessian +hessian=ComputeHessian(index,x,y,md.vel_obs,'node'); + +%Compute metric +if length(md.gridonwater)==md.numberofgrids, + pos=find(md.gridonwater); +else + pos=[]; +end +metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos); + +%Get Areas +areas=GetAreas(index,x,y); + +%length edges vectors +e1x=[x(index(:,2))-x(index(:,1))]; +e1y=[y(index(:,2))-y(index(:,1))]; +e2x=[x(index(:,3))-x(index(:,2))]; +e2y=[y(index(:,3))-y(index(:,2))]; +e3x=[x(index(:,1))-x(index(:,3))]; +e3y=[y(index(:,1))-y(index(:,3))]; + +%metric of each the 3 grids for each element +M1=metric(index(:,1),:); +M2=metric(index(:,2),:); +M3=metric(index(:,3),:); + +%Get edge length in the metric +L1=1/2*(sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y))+sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y))); +L2=1/2*(sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y))+sqrt(e2x.*(M2(:,1).*e2x+M2(:,2).*e2y)+e2y.*(M2(:,2).*e2x+M2(:,3).*e2y))); +L3=1/2*(sqrt(e1x.*(M1(:,1).*e1x+M1(:,2).*e1y)+e1y.*(M1(:,2).*e1x+M1(:,3).*e1y))+sqrt(e3x.*(M3(:,1).*e3x+M3(:,2).*e3y)+e3y.*(M3(:,2).*e3x+M3(:,3).*e3y))); + +%area in the metric +V=1/3*areas.*(sqrt(M1(:,1).*M1(:,3)-M1(:,2).^2)+sqrt(M2(:,1).*M2(:,3)-M2(:,2).^2)+sqrt(M3(:,1).*M3(:,3)-M3(:,2).^2)); + +%compute quality: +quality=4*sqrt(3)*V./(L1+L2+L3); + +%compute error +a=hessian(:,1); b=hessian(:,2); d=hessian(:,3); +a=a(index)*[1;1;1]/3; +b=b(index)*[1;1;1]/3; +d=d(index)*[1;1;1]/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)); +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.gridonwater)==md.numberofgrids; + pos=find(md.gridonwater); + lambda1(pos)=0; + lambda2(pos)=0; +end +lambda1=lambda1(index)*[1;1;1]/3; +lambda2=lambda2(index)*[1;1;1]/3; + +lambdamax=max(lambda1,lambda2); +hmax=max(max(sqrt(e1x.^2+e1y.^2),sqrt(e2x.^2+e2y.^2)),sqrt(e3x.^2+e3y.^2)); +epsilon=scale*hmax.^2.*lambdamax; + +%display +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')); Index: /issm/trunk/src/m/utils/Mesh/YamsCall.m =================================================================== --- /issm/trunk/src/m/utils/Mesh/YamsCall.m (revision 1341) +++ /issm/trunk/src/m/utils/Mesh/YamsCall.m (revision 1342) @@ -19,89 +19,23 @@ 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,9 +66,7 @@ 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