estimate.m

C++编译指导

function eta = estimate(node,elem,u)
[dudx,dudy] = ZZ(node,elem,u);
[dudxx,dudxy] = sgrad(node,elem,dudx);
[dudyx,dudyy] = sgrad(node,elem,dudy);
eta = abs(dudxx) + abs(dudxy) + abs(dudyx) + abs(dudyy);

% --------------------------------------------------------
function [dudx,dudy] = ZZ(node,elem,u)
N = size(node,1); NT = size(elem,1);
[sdudx,sdudy,area] = sgrad(node,elem,u);
patcharea = sparse(elem,ones(NT,3),area*[1,1,1],N,1);
dudx = sparse(elem,ones(NT,3),sdudx*[1,1,1],N,1);
dudy = sparse(elem,ones(NT,3),sdudy*[1,1,1],N,1);
dudx = dudx./patcharea;
dudy = dudy./patcharea;

function [sdudx,sdudy,area] = sgrad(node,elem,u)
% this function compute the scaled grad of u on the mesh
N = size(node,1); NT = size(elem,1);
loce(:,1,:) = node(elem(:,3),:)-node(elem(:,2),:);
loce(:,2,:) = node(elem(:,1),:)-node(elem(:,3),:);
loce(:,3,:) = node(elem(:,2),:)-node(elem(:,1),:);
area = 0.5*abs(-loce(:,3,1).*loce(:,2,2)+loce(:,3,2).*loce(:,2,1));
sdudx = -0.5*( u(elem(:,1)).*loce(:,1,2) + u(elem(:,2)).*loce(:,2,2) ...
           + u(elem(:,3)).*loce(:,3,2) );
sdudy = 0.5*( u(elem(:,1)).*loce(:,1,1) + u(elem(:,2)).*loce(:,2,1) ...
           + u(elem(:,3)).*loce(:,3,1) );