estimate.m

C++编译指导

function [ eta, res, osc]= estimate(mesh,eps,a1,a2,b,f,g_N,bta)
% ESTIMATE computes residual error estimator on each element
% 
% USAGE
%  [ eta, res, osc]= estimate(mesh,eps,a1,a2,b,f,bta)
%
% INPUT 
%    mesh:  current mesh
%    data:  eps,a1,a2,b,f,bta
% OUTPUT
%    eta:  error estimator on each element
%    res:  residual element error estimator on each element
%    osc:  oscillation error estimator on each element

N = size(mesh.node,1); NT = size(mesh.elem,1);
eta=sparse(NT,1);res=sparse(NT,1);osc=sparse(NT,1);center=sparse(NT,1);
tlength=sparse(NT,1);
dualEdge = sparse(mesh.elem(:,[1,2,3]),mesh.elem(:,[2,3,1]),[1:NT,1:NT,1:NT]);


% %
% %
if(size(mesh.Neumann)~=0)
     edge = [mesh.Neumann(:,[1,2]);mesh.Dirichlet(:,[1,2])];
     edge = unique(sort(edge,2),'rows');
     edgeN = mesh.Neumann(:,[1,2]);
     edgeN = unique(sort(edgeN,2),'rows');
     NN=max(max(edge(:,1),max(edge(:,2))));
     d2pN=sparse(NN,NN);
   
      NN=size(edgeN,1);
    for i=1:NN
       d2pN(edgeN(i,1),edgeN(i,2))=i; d2pN(edgeN(i,2),edgeN(i,1))=i;
    end 
end
% %
% %

center = ( mesh.node(mesh.elem(:,1),:) + mesh.node(mesh.elem(:,2),:)...
                                  + mesh.node(mesh.elem(:,3),:) )/3;
ve(:,:,1) = mesh.node(mesh.elem(:,3),:)-mesh.node(mesh.elem(:,2),:);
ve(:,:,2) = mesh.node(mesh.elem(:,1),:)-mesh.node(mesh.elem(:,3),:);
ve(:,:,3) = mesh.node(mesh.elem(:,2),:)-mesh.node(mesh.elem(:,1),:);
area = 0.5*abs(-ve(:,1,3).*ve(:,2,2)+ve(:,2,3).*ve(:,1,2));
for i=1:NT
tlength(i)=max([ve(i,1,1)^2 + ve(i,2,1)^2 , ve(i,1,2)^2 + ve(i,2,2)^2 ,ve(i,1,3)^2 + ve(i,2,3)^2]);
end 
tlength=4*(area.^2)./tlength;   
% for i=1:NT
% tlength(i)=min([ve(i,1,1)^2 + ve(i,2,1)^2 , ve(i,1,2)^2 + ve(i,2,2)^2 ,ve(i,1,3)^2 + ve(i,2,3)^2]);
% end 

%--------------------------------------------------------------------------
%residual element error estimator
%--------------------------------------------------------------------------
for i=1:3
    res = res +(mesh.solu(mesh.elem(:,i)).*(-feval(a1,center)'.*ve(:,2,i)...
                                +feval(a2,center)'.*ve(:,1,i)))./(2*area)...
                           -(feval(b,center)'.*mesh.solu(mesh.elem(:,i)))/3;
end
res =feval(f,center)'+res;
%--------------------------------------------------------------------------
%edge error estimator
%--------------------------------------------------------------------------
re=sparse(NT,1);
for k=1:NT
    for i=1:3
        if (mod(i+2,3)==0)
            e1=3;
        else
            e1=mod(i+2,3);
        end
        if (mod(i+1,3)==0)
            e2=3;
        else
            e2=mod(i+1,3);
        end    
        ct=dualEdge(mesh.elem(k,e1),mesh.elem(k,e2));
      if(ct~=0)
             e1=0;e2=0; 
          for j=1:3
               e1=e1+(mesh.solu(mesh.elem(k,j))*(-ve(k,2,j)*ve(k,2,i)-ve(k,1,j)...
                  *ve(k,1,i)))/(2*area(k)); 
               e2=e2+(mesh.solu(mesh.elem(ct,j))*(-ve(ct,2,j)*ve(k,2,i)-ve(ct,1,j)...
                  *ve(k,1,i)))/(2*area(ct)); 
          end
          re(k) =re(k)+(sqrt(tlength(k))/sqrt(ve(k,1,i)^2 + ve(k,2,i)^2 ))*(eps*(e1-e2))^2/sqrt(ve(k,1,i)^2 + ve(k,2,i)^2 );
%--------------------------------------------------------------------------
%neumann condition
       else
          if(size(mesh.Neumann)~=0)
             if( d2pN(mesh.elem(k,e1),mesh.elem(k,e2))~=0 )    
                midpoint=1/2*(mesh.node(mesh.elem(k,e1),:)+mesh.node(mesh.elem(k,e2),:));
                e1=0;e2=0; 
                for j=1:3
                    e1=e1+(mesh.solu(mesh.elem(k,j))*(-ve(k,2,j)*ve(k,2,i)-ve(k,1,j)...
                        *ve(k,1,i)))/(2*area(k)); 
                end
                re(k)=re(k)+2*(sqrt(tlength(k))/sqrt(ve(k,1,i)^2 + ve(k,2,i)^2 ))*(feval(g_N, midpoint)*sqrt(ve(k,1,i)^2 +ve(k,2,i)^2 )...
                       +eps*e1)^2/sqrt(ve(k,1,i)^2 + ve(k,2,i)^2 );
             end  
          end
%--------------------------------------------------------------------------

      end  
    end
end
%--------------------------------------------------------------------------
%oscallition error ,because I donnot know how to integrate the L2
%projection,so use the following simple method,but .......
%--------------------------------------------------------------------------
for i=1:3
    osc=osc+(-feval(a1,mesh.node(mesh.elem(:,1),:))'.*ve(:,2,i)...
              +feval(a2,mesh.node(mesh.elem(:,1),:))'.*ve(:,1,i))/3 ...
                            +feval(b,mesh.node(mesh.elem(:,1),:))'/3;
end
osc=area.*((res-osc).^2);
res =area.*(res.^2 );
if (bta==0)
 eta=(tlength(:)/eps).*res + eps*(sqrt(tlength(:)).*re)/2;  
else 
  eta=min( tlength(:)/eps,ones(NT,1)/bta ).*res +min( sqrt(tlength(:))/eps,(1/(sqrt(bta*eps)))*ones(NT,1)).*re/2;
end    
res=(tlength(:)/eps).*res;
osc=(tlength(:)/eps).*res;