二次元P2和三次元P3Neumann边界条件的处理

二次元P2Neumann边界条件的处理

    % Neumann boundary condition
    if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 4;   % default order
        end
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        % quadratic bases on an edge (1---3---2)
        bdphi = zeros(nQuadgN,3);        
        bdphi(:,1) = (2*lambdagN(:,1)-1).*lambdagN(:,1);
        bdphi(:,2) = (2*lambdagN(:,2)-1).*lambdagN(:,2);
        bdphi(:,3) = 4*lambdagN(:,1).*lambdagN(:,2);
        % length of edge
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),3);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNp = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNp*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNp*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNp*bdphi(pp,3); % interior bubble
        end
        % update RHS
        ge = ge.*repmat(el,1,3);        
        b(1:N) = b(1:N) + accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b(N+Neumannidx) = b(N+Neumannidx) + ge(:,3);
    end

三次元P3Neumann边界条件的处理

 1     % Neumann boundary condition
 2     if ~isempty(pde.g_N) && ~isempty(Neumann) && ~(isnumeric(pde.g_N) && (pde.g_N == 0))
 3         if ~isfield(option,'gNquadorder')
 4             option.gNquadorder = 6;  
 5         end
 6         [lambdagN,weightgN] = quadpts1(option.gNquadorder);
 7         nQuadgN = size(lambdagN,1);
 8         % quadratic bases (1---3---4--2)
 9         bdphi = zeros(nQuadgN,4);        
10         bdphi(:,1) = 0.5*(3*lambdagN(:,1)-1).*(3*lambdagN(:,1)-2).*lambdagN(:,1); 
11         bdphi(:,2) = 0.5*(3*lambdagN(:,2)-1).*(3*lambdagN(:,2)-2).*lambdagN(:,2);
12         bdphi(:,3) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,1)-1);
13         bdphi(:,4) = 9/2*lambdagN(:,1).*lambdagN(:,2).*(3*lambdagN(:,2)-1);
14         % length of edge
15         el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
16         ge = zeros(size(Neumann,1),4);
17         for pp = 1:nQuadgN
18             ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
19                  + lambdagN(pp,2)*node(Neumann(:,2),:);
20             gNp = pde.g_N(ppxy);
21             ge(:,1) = ge(:,1) + weightgN(pp)*gNp*bdphi(pp,1);
22             ge(:,2) = ge(:,2) + weightgN(pp)*gNp*bdphi(pp,2);
23             ge(:,3) = ge(:,3) + weightgN(pp)*gNp*bdphi(pp,3);    
24             ge(:,4) = ge(:,4) + weightgN(pp)*gNp*bdphi(pp,4);
25         end
26         ge = ge.*repmat(el,1,4);
27         b(1:N) = b(1:N) + accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
28         b(N+2*Neumannidx-1) = b(N+2*Neumannidx-1) + ge(:,3);
29         b(N+2*Neumannidx)   = b(N+2*Neumannidx) + ge(:,4);
30     end