0001 function md2=meshadaptation(md,field,scale,epsilon)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
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
0025 numberofelements=md.numberofelements;
0026 numberofgrids=md.numberofgrids;
0027 index=md.elements;
0028 x=md.x; y=md.y; z=md.z;
0029
0030
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
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
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
0050 grad_elx=sum(field(index).*alpha,2);
0051 grad_ely=sum(field(index).*beta,2);
0052
0053
0054 areas=area(md);
0055
0056
0057 weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofgrids,1);
0058
0059
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
0066 hessian=[sum(gradx(index).*alpha,2) sum(grady(index).*beta,2) sum(grady(index).*alpha,2)];
0067
0068
0069 hmin=min(sqrt(areas))/scale;
0070 hmax=max(sqrt(areas));
0071
0072
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
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
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092 metric=20*1./sqrt(lambda1.^2+lambda2.^2);
0093
0094
0095 metric(find(metric<10^6))=10^6;
0096
0097
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))