meshadaptation

PURPOSE ^

MESHADAPTATION - remesh a model for a given field

SYNOPSIS ^

function md2=meshadaptation(md,field,scale,epsilon)

DESCRIPTION ^

MESHADAPTATION - remesh a model for a given field

   remesh a model for a given field, use the hessian computation to spread the error on the element
   scale: 1.2 large scale leads to higher refining.
   epsilon: .5  large epsilon will refine fewer elements

   Usage:
      md2=meshadaptation(md,field,scale,epsilon)

   Example:
      md2=meshadaptation(md,md.vel,1.2,0.5);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md2=meshadaptation(md,field,scale,epsilon)
0002 %MESHADAPTATION - remesh a model for a given field
0003 %
0004 %   remesh a model for a given field, use the hessian computation to spread the error on the element
0005 %   scale: 1.2 large scale leads to higher refining.
0006 %   epsilon: .5  large epsilon will refine fewer elements
0007 %
0008 %   Usage:
0009 %      md2=meshadaptation(md,field,scale,epsilon)
0010 %
0011 %   Example:
0012 %      md2=meshadaptation(md,md.vel,1.2,0.5);
0013 
0014 %some checks
0015 if ~strcmpi(md.type,'2d')
0016     error('meshadaptation error message: adaptation for 3d meshes not implemented yet')
0017 end
0018 if length(field)~=md.numberofgrids
0019     error('meshadaptation error message: input field length shoud be numberofgrids')
0020 end
0021 
0022 disp(sprintf('      metric computation') )
0023 
0024 %load some variables (it is much faster if the variab;es are loaded from md once for all)
0025 numberofelements=md.numberofelements;
0026 numberofgrids=md.numberofgrids;
0027 index=md.elements;
0028 x=md.x; y=md.y; z=md.z;
0029 
0030 %initialization
0031 alpha=zeros(md.numberofelements,3);
0032 beta=zeros(md.numberofelements,3);
0033 gradx=zeros(md.numberofgrids,1);
0034 grady=zeros(md.numberofgrids,1);
0035 metric=zeros(md.numberofelements,1);
0036 
0037 %build some usefull variables
0038 field=field(:);
0039 line=index(:);
0040 summation=1/3*ones(3,1);
0041 linesize=3*numberofelements;
0042 x1=x(index(:,1)); x2=x(index(:,2)); x3=x(index(:,3)); y1=y(index(:,1)); y2=y(index(:,2)); y3=y(index(:,3));
0043 
0044 %compute nodal functions coefficients N(x,y)=alpha x + beta y + gamma
0045 invdet=1./(x1.*(y2-y3)-x2.*(y1-y3)+x3.*(y1-y2));
0046 alpha=[invdet.*(y2-y3) invdet.*(y3-y1) invdet.*(y1-y2)];
0047 beta=[invdet.*(x3-x2) invdet.*(x1-x3) invdet.*(x2-x1)];
0048 
0049 %Compute gradient for each element
0050 grad_elx=sum(field(index).*alpha,2); 
0051 grad_ely=sum(field(index).*beta,2);
0052 
0053 %compute the volume of each element
0054 areas=area(md);
0055 
0056 %update weights that holds the volume of all the element holding the grid i
0057 weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
0058 
0059 %Compute gradient for each grid (average of the elements around)
0060 gradx=sparse(line,ones(linesize,1),repmat(areas.*grad_elx,3,1),numberofgrids,1);
0061 grady=sparse(line,ones(linesize,1),repmat(areas.*grad_ely,3,1),numberofgrids,1);
0062 gradx=gradx./weights;
0063 grady=grady./weights;
0064 
0065 %Compute hessian for each element
0066 hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*beta,2) sum(grady(index).*alpha,2)];
0067 
0068 %Compute the new metric
0069 hmin=min(sqrt(areas))/scale;
0070 hmax=max(sqrt(areas));
0071 
0072 %first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,3); hessian(i,3)  hessian(i,2)]
0073 a=hessian(:,1); b=hessian(:,3); d=hessian(:,2);
0074 lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2));
0075 lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2));
0076 
0077 %Modify the eigen values to control the shape of the elements
0078 lambda1=min(max(abs(lambda1)/epsilon,1/hmax^2),1/hmin^2);
0079 lambda2=min(max(abs(lambda2)/epsilon,1/hmax^2),1/hmin^2);
0080 
0081 %find the corresponding eigen vectors of the initial eigen values
0082 %NOT USED FOR NOW
0083 %div1=sqrt(1+(b./(lambda1-d)).^2);
0084 %div2=sqrt(1+(b./(lambda2-d)).^2);
0085 %vector1=[1./div1 (b./(lambda1-d))./div1];
0086 %vector2=[1./div2 (b./(lambda2-d))./div2];
0087 
0088 %compute the metric
0089 %R=[vector1(:,1) vector2(:,1) vector1(:,2) vector2(:,2)];
0090 %det=vector1(:,1).*vector2(:,2)-vector1(:,2).*vector1(:,2);
0091 %Rinv=[vector2(:,2)./det   -vector2(:,1)./det   -vector1(:,2)./det   vector1(:,1)./det];
0092 metric=20*1./sqrt(lambda1.^2+lambda2.^2);
0093 
0094 %ok, metric can't go under 1km resolution.
0095 metric(find(metric<10^6))=10^6;
0096 
0097 %Remesh with this new metric
0098 disp(sprintf('      initial number of element: %i', md.numberofelements))
0099 md2=meshrefine(md,full(metric));
0100 disp(sprintf('      new number of elements:    %i', md2.numberofelements))

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003