build_operators.m

ZEUS is a family of Eulerian (grid based) Magneto-Hydrodynamic codes (MHD) for use in astrophysics,

Gx = sparse(i,j,sx,numn,numn);

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function L = del(data,bc)

% Build the Laplacian operator using a hybrid FE/FV method. The gradients
% are obtained at the cell centroids and then used to construct the face
% fluxes on the median cell edges.
%
% Only uses nearest neighbours and seems to be fully 2nd order accuracte.

p    = data.p;          % Nodes
e    = data.e;          % Edges
t    = data.t;          % Triangles
n2n  = data.n2n;        % Node to node connectivity
n2e  = data.n2e;        % Node to edge connectivity
e2t  = data.e2t;        % Edge to triangle connectivity
bnd  = data.bnd;        % True for boundary nodes
hnx  = data.hnx;        % Median cell
hny  = data.hny;        % normals
numn = size(p,1);       % Number of nodes


% FORM FEM GRADIENTS
% Evaluate centroidal gradients (piecewise-linear interpolants)
x23 = p(t(:,2),1)-p(t(:,3),1);  y23 = p(t(:,2),2)-p(t(:,3),2);
x21 = p(t(:,2),1)-p(t(:,1),1);  y21 = p(t(:,2),2)-p(t(:,1),2);

% Denominators
detx = x23.*y21-x21.*y23;
dety = y23.*x21-y21.*x23;

% xy gradient coefficients at vertices [n1,n2,n3]
sx = [y23./detx, (y21-y23)./detx, -y21./detx];
sy = [x23./dety, (x21-x23)./dety, -x21./dety];


% FORM LAPLACIAN
% Max number of entries in [i,j,s]
num = 6*sum(sum(n2n>0,2));

% Pre-alloc
i = repmat(uint32(0),num,1); j = i; s = double(i); 

next = 1;
for k = 1:numn
    if ~bnd(k)||(bc(k)==0)      % Node k not Dirichlet
        
        % Loop around neighbours of node k
        m = 1;
        while n2e(k,m)>0
            
            cn = n2n(k,m);      % Neighbour
            ce = n2e(k,m);      % Current edge
            t1 = e2t(ce,1);
            t2 = e2t(ce,2);
            
            % ce has 2 median edges associated, deal with separately
            if k==e(ce,1)
                hnx1 = hnx(ce,1); hny1 = hny(ce,1);
                hnx2 = hnx(ce,2); hny2 = hny(ce,2);
            else
                hnx1 = -hnx(ce,1); hny1 = -hny(ce,1);
                hnx2 = -hnx(ce,2); hny2 = -hny(ce,2);
            end

            % Loop over associated median edges and assemble flux
            for z = 1:3
                % Edge 1
                i(next) = k;
                j(next) = t(t1,z);
                s(next) = hnx1*sx(t1,z) + hny1*sy(t1,z);
                next    = next+1;
                % Edge 2
                if t2>0
                    i(next) = k;
                    j(next) = t(t2,z);
                    s(next) = hnx2*sx(t2,z) + hny2*sy(t2,z);
                    next    = next+1;  
                end
            end
            
            % Counter
            m = m+1;
            
        end
        
    else    % Node k is Dirichlet
        i(next) = k;
        j(next) = k;
        s(next) = 1;
        next    = next+1;
    end
end

% Sparse operator
L = sparse(double(i(1:next-1)),double(j(1:next-1)),s(1:next-1),numn,numn);

return


% Old FV Laplacian
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function L = del(data,bc,A,B,C)
% 
% % Build Laplacian operators via a 2nd order FV discretisation on the median
% % mesh. Incorporate BC's directly into the coefficient matrix.
% 
% p    = data.p;          % Nodes
% e    = data.e;          % Edges
% n2n  = data.n2n;        % Node to node connectivity
% n2e  = data.n2e;        % Node to edge connectivity
% e2t  = data.e2t;        % Edge to triangle connectivity
% bnd  = data.bnd;        % True for boundary nodes
% d    = data.d;          % Edge length
% hnx  = data.hnx;        % Median cell
% hny  = data.hny;        % normals
% numn = size(p,1);       % Number of nodes
% 
% % Max number of entries in [i,j,s]
% numN = sum(n2n>0,2);
% num  = sum(numN.*(2*numN+2));
% 
% % Pre-alloc
% i = repmat(uint32(0),num,1); j = i; s = double(i); 
% 
% next = 1;
% for k = 1:numn
%     if ~bnd(k)||(bc(k)==0)      % Node k not Dirichlet
%         
%         % Loop around neighbours of node k
%         m = 1;
%         while n2e(k,m)>0
%             
%             cn = n2n(k,m);      % Neighbour
%             ce = n2e(k,m);      % Current edge
%             t2 = e2t(ce,2);
%             
%             % xy nodes
%             xk = p(k,1);  yk = p(k,2);
%             xc = p(cn,1); yc = p(cn,2);
%             
%             % Normal vector for ce (different from median normals)
%             nxw = (xc-xk)/d(ce);
%             nyw = (yc-yk)/d(ce);
%             
%             % ce has 2 median edges associated, deal with separately
%             if k==e(ce,1)
%                 hnx1 = hnx(ce,1); hny1 = hny(ce,1);
%                 hnx2 = hnx(ce,2); hny2 = hny(ce,2);
%             else
%                 hnx1 = -hnx(ce,1); hny1 = -hny(ce,1);
%                 hnx2 = -hnx(ce,2); hny2 = -hny(ce,2);
%             end
%             % Coefficients
%             ax1 = hnx1*(1-nxw^2); bx1 = hnx1*nxw*nyw;
%             ax2 = hnx2*(1-nxw^2); bx2 = hnx2*nxw*nyw;
%             ay1 = hny1*(1-nyw^2); by1 = hny1*nxw*nyw;
%             ay2 = hny2*(1-nyw^2); by2 = hny2*nxw*nyw;
%             
%             % Centre coefficients
%             sum1 = -(hnx1*nxw+hny1*nyw)/d(ce);
%             if t2>0
%                 sum1 = sum1 - (hnx2*nxw+hny2*nyw)/d(ce);    
%             end
%             sum2 = -sum1;
%             
%             % Gradient around k
%             Ak = 0.5*A(k); Bk = 0.5*B(k); Ck = 0.5*C(k);
%             z = 1;
%             while n2n(k,z)>0 
%                 % Neighbour
%                 nn = n2n(k,z);
%                 % xy increments
%                 x = xk-p(nn,1); 
%                 y = yk-p(nn,2);
%                 w = 1/(x^2+y^2);
%                 % Gradient coefficients
%                 gy = Bk*w*x-Ak*w*y;
%                 gx = Ck*w*x-Bk*w*y;
%                 % Centre term
%                 temp = (ax1*-gx)-(bx1*gy)+(ay1*gy)-(by1*-gx);
%                 sum1 = sum1 - temp;
%                 if t2>0
%                     cont = (ax2*gx)-(bx2*-gy)+(ay2*-gy)-(by2*gx);
%                     sum1 = sum1 + cont; 
%                     temp = temp - cont;
%                 end
%                 % Neighbour term
%                 i(next) = k;
%                 j(next) = nn;
%                 s(next) = temp;
%                 next    = next+1;
%                 % Counter
%                 z = z+1;
%             end
%             
%             % Gradient around cn
%             Ac = 0.5*A(cn); Bc = 0.5*B(cn); Cc = 0.5*C(cn);
%             z = 1;
%             while n2n(cn,z)>0
%                 % Neighbour
%                 nn = n2n(cn,z);
%                 % xy increment
%                 x = xc-p(nn,1); 
%                 y = yc-p(nn,2);
%                 w = 1/(x^2+y^2);
%                 % Gradient coefficients
%                 gy = Bc*w*x-Ac*w*y;
%                 gx = Cc*w*x-Bc*w*y;
%                 % Centre term
%                 temp = (ax1*-gx)-(bx1*gy)+(ay1*gy)-(by1*-gx);
%                 sum2 = sum2 - temp;
%                 if t2>0
%                     cont = (ax2*gx)-(bx2*-gy)+(ay2*-gy)-(by2*gx);
%                     sum2 = sum2 + cont; 
%                     temp = temp - cont;
%                 end
%                 % Neighbour term
%                 i(next) = k;
%                 j(next) = nn;
%                 s(next) = temp;
%                 next    = next+1;
%                 % Counter
%                 z = z+1;
%             end
%             
%             % Centre term
%             i(next) = k;
%             j(next) = k;
%             s(next) = sum1;
%             next    = next+1;
%             % Neighbour term
%             i(next) = k;
%             j(next) = cn;
%             s(next) = sum2;
%             next    = next+1;
%             
%             % Counter
%             m = m+1;
%             
%         end
%         
%     else    % Node k is Dirichlet
%         i(next) = k;
%         j(next) = k;
%         s(next) = 1;
%         next    = next+1;
%     end
% end
% 
% % Sparse operator
% L = sparse(double(i(1:next-1)),double(j(1:next-1)),s(1:next-1),numn,numn);
% 
% return
% 
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [A,B,C] = least_squares(data)
% 
% % Form least squares gradient coefficients.
% 
% p    = data.p;          % Nodes
% n2n  = data.n2n;        % Node to node connectivity
% numn = size(p,1);       % Number of nodes
% 
% % LEAST SQUARES GRADIENT COEFFICIENTS (inverse distance weighted)
% A = repmat(0,numn,1); B = A; C = A;
% for k = 1:numn  
%     % Loop around neighbours of node k
%     a = 0; b = 0; c = 0; m = 1;
%     while n2n(k,m)>0
%         % xy increments
%         x = p(k,1)-p(n2n(k,m),1); 
%         y = p(k,2)-p(n2n(k,m),2);
%         w = 1/(x^2+y^2);
%         % Coefficients
%         a = a+w*x^2; b = b+w*x*y; c = c+w*y^2;
%         % Counter
%         m = m+1;
%     end
%     % Save
%     A(k) = a; B(k) = b; C(k) = c;
% end
% 
% % Divide by det
% detM = A.*C-B.^2; 
% A    = A./detM; 
% B    = B./detM; 
% C    = C./detM;
% 
% return