二次元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