Contents

function [A,F] = getmatrix(node,elem,d,b,c,f,qOrder)

GETMATRIX matrix discretization a PDE problem.

[A,F] = getmatrix(node,elem,d,b,c,f,[qd qb qc qf]) get the matrix equation Au=F of the linear finite element discretization of the following elliptic equation -div(d*grad(u)) + b*grad(u) + c*u = f, on a mesh described by node and elem. Quadrature rules with order qd,qb,qc,qf are used for different parts of the PDE. The qudarature order can be as high as 5. When qc==1, mass lumping is applied to assemble the mass matrix for the reaction term.

[A,F] = getmatrix(node,elem,d,b,c,f) computes A and F using quadrature rule qd=1,qb=1,qc=2,and qf=2.

The diffusion coefficient d could be a scalar or a symmetric and uniformly positive definite tensor (2 by 2 matrix) function. When d is a tensor function, the four entries are represented by d_11 = d(:,1), d_12 = d(:,2), d_21 = d(:,3), d_22 = d(:,4). The convection velocity b is a vector function. The reaction coefficient c is a scalar function and so is f.

The coefficients d, b, c, and f of the PDE problem can be given in a wide variety of ways:

See also ellipticpde, getbdcondition, getquadpts

% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = size(node,1); NT = size(elem,1); A = sparse(N,N);

Set up quadrature rules for different parts

if nargin<=3
    b = []; c = []; f = [];
end
if (b==0), b=[]; end
if (c==0), c=[]; end
if (f==0), f=[]; end
qd = 1;   qb = 1;   qc = 2;   qf = 2; % default choice
if naragin==7
   qd = qOrder(1);
   qb = qOrder(2);
   qc = qOrder(3);
   qf = qOrder(4);
end
if qc == 1
    masslumping = true;
end

Compute geometric quantities and gradient of local basis

ve1 = node(elem(:,3),:)-node(elem(:,2),:);
ve2 = node(elem(:,1),:)-node(elem(:,3),:);
ve3 = node(elem(:,2),:)-node(elem(:,1),:);
area = 0.5*abs(-ve3(:,1).*ve2(:,2) + ve3(:,2).*ve2(:,1));
Dphi(1:NT,:,1) = [-ve1(:,2)./(2*area), ve1(:,1)./(2*area)];
Dphi(1:NT,:,2) = [-ve2(:,2)./(2*area), ve2(:,1)./(2*area)];
Dphi(1:NT,:,3) = [-ve3(:,2)./(2*area), ve3(:,1)./(2*area)];

$$\nabla \phi_i = l_i^{\bot}/(2|\tau|)$$

clear ve1 ve2 ve3

Assemble matrices

for i = 1:3
    for j = 1:3

diffusion $D_{ij} = \int _{\tau} (d\nabla \phi _j)\cdot \nabla \phi _i\, dxdy$

        if isreal(d)  % d is an array
            if (size(d,2) == 1) && (size(d,1)==NT) % piecewise-constant scalar
                Dij = d.*(Dphi(:,1,i).*Dphi(:,1,j) + Dphi(:,2,i).*Dphi(:,2,j));
            elseif (size(d,2) == 4) && (size(d,1)==NT) % piecewise-constant tensor
                Dij = d(:,1).*Dphi(:,1,i).*Dphi(:,1,j) ...
                    + d(:,2).*Dphi(:,1,i).*Dphi(:,2,j) ...
                    + d(:,3).*Dphi(:,2,i).*Dphi(:,1,j) ...
                    + d(:,4).*Dphi(:,2,i).*Dphi(:,2,j);
            elseif (size(d,2) == 1) && (size(d,1)==N) % piecewise-linear scalar
                Dij = (d(elem(:,1))+ d(elem(:,2))+ d(elem(:,3)))/3.*...
                    (Dphi(:,1,i).*Dphi(:,1,j) + Dphi(:,2,i).*Dphi(:,2,j));
            elseif (size(d,2) == 4) && (size(d,1)==N) % piecewise-linear tensor
                dt = (d(elem(:,1),:)+ d(elem(:,2),:)+ d(elem(:,3),:))/3;
                Dij = dt(:,1).*Dphi(:,1,i).*Dphi(:,1,j) ...
                    + dt(:,2).*Dphi(:,1,i).*Dphi(:,2,j) ...
                    + dt(:,3).*Dphi(:,2,i).*Dphi(:,1,j) ...
                    + dt(:,4).*Dphi(:,2,i).*Dphi(:,2,j);
            else
                error('The diffusion coefficient is either a scalar or a 2 by 2 matrix');
            end
        else          % d is a function
            % compute quadrature points
            [lambda,weight] = getquadpts(qd);
            nQuad = size(lambda,1);
            % compute element-wise coefficient d
            dt = zeros(NT,size(d([0 0]),2));
            for p = 1:nQuad
                % quadrature points in the x-y coordinate
                pxy = lambda(p,1)*node(elem(:,1),:) ...
                    + lambda(p,2)*node(elem(:,2),:) ...
                    + lambda(p,3)*node(elem(:,3),:);
                dt = dt + weight(p)*d(pxy);
            end
            if size(dt,2) == 1     % scalar function
                Dij = dt.*(Dphi(:,1,i).*Dphi(:,1,j) + Dphi(:,2,i).*Dphi(:,2,j));
            elseif size(d,2) == 4  % tensor function
                Dij = dt(:,1).*Dphi(:,1,i).*Dphi(:,1,j) ...
                    + dt(:,2).*Dphi(:,1,i).*Dphi(:,2,j) ...
                    + dt(:,3).*Dphi(:,2,i).*Dphi(:,1,j) ...
                    + dt(:,4).*Dphi(:,2,i).*Dphi(:,2,j);
            end
        end
        Aij = Dij;

convection $B_{ij} = \int _{\tau} \phi_i b \cdot \nabla \phi_j \, dxdy$

        if ~isempty(b)
            if isreal(b)  % b is an array
                if size(b,1) == NT  % piecewise-constant vector
                    % 1-pt quadrature. phi_i(center) = 1/3;
                    Bij = (b(:,1).*Dphi(:,1,j)+ b(:,2).*Dphi(:,2,j))/3;
                elseif size(b,1) == N
                    % 3 middle points quadrature. phi_i = [1/2 1/2 0].
                    k = setdiff([1 2 3], i);
                    bt = 0.5*0.5*(b(elem(:,i),:) + b(elem(:,k(1)),:)) ...
                       + 0.5*0.5*(b(elem(:,i),:) + b(elem(:,k(2)),:));
                    Bij = (bt(:,1).*Dphi(:,1,j)+ bt(:,2).*Dphi(:,2,j))/3;
                else
                    error('The convection coefficient is a vector function');
                end
            else         % b is a function
               % compute quadrature points
                [lambda,weight] = getquadpts(qb);
                phi = lambda;
                nQuad = size(lambda,1);
                % compute element-wise coefficient b
                bt = zeros(NT,2);
                for p = 1:nQuad
                    % quadrature points in the x-y coordinate
                    pxy = lambda(p,1)*node(elem(:,1),:) ...
                        + lambda(p,2)*node(elem(:,2),:) ...
                        + lambda(p,3)*node(elem(:,3),:);
                    bt = bt + weight(p)*b(pxy)*phi(p,i);
                end
                Bij = bt(:,1).*Dphi(:,1,j) + bt(:,2).*Dphi(:,2,j);
            end
            Aij = Aij + Bij;
        end

reaction $C_{ij} = \int _{\tau} c \phi_j\phi_i\, dxdy$

        if ~isempty(c) && ~masslumping
            if isreal(c)  % c is an array
                if size(c,1) == NT % piecewise-constant
                    Cij = c.*area*((i==j)+1)/12;
                elseif size(c,1) == N % piecewise-constant
                    % 1-pt quadrature. Not exact.
                    Cij = (c(elem(:,1)) + c(elem(:,2)) + c(elem(:,3)))/3 ...
                          .*area*((i==j)+1)/12;
                end
            else          % c is a function
                % compute quadrature points
                [lambda,weight] = getquadpts(qc);
                phi = lambda;
                nQuad = size(lambda,1);
                % compute element-wise coefficient C
                Cij = zeros(NT,1);
                for p = 1:nQuad
                    % quadrature points in the x-y coordinate
                    pxy = lambda(p,1)*node(elem(:,1),:) ...
                        + lambda(p,2)*node(elem(:,2),:) ...
                        + lambda(p,3)*node(elem(:,3),:);
                    Cij = Cij + weight(p)*phi(p,i)*phi(p,j)*c(pxy);
                end
            end
            Aij = Aij + Cij;
        end

assemble element-wise quantity to node-wise one

        Aij = Aij.*area;
        A = A + sparse(elem(:,i),elem(:,j),Aij,N,N);

    end
end

assembel lumped mass matrix

if masslumping
   if isreal(c) && (length(c) == N)
      M = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,[N,1]);
      A = A + spdiags(c.*M,0,N,N);
   elseif isreal(c) && (length(c) == NT)
      ca = c.*area;
      M = accumarray([elem(:,1);elem(:,2);elem(:,3)],[ca;ca;ca]/3,[N,1]);
      A = A + spdiags(c.*M,0,N,N);
   else
      M = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,[N,1]);
      A = A + spdiags(c(node).*M,0,N,N);
   end
end

Assemble the right hand side

$F_i = \int _{\tau} f\phi_i\, dxdy$

if isreal(f) % f is an array
    if size(f,2)~=1
        f = f';   % transfer f to a column vector
    end
    if length(f)== N        % piecewise linear
       M = accumarray([elem(:,1);elem(:,2);elem(:,3)],[area;area;area]/3,[N,1]);
           F = f.*M;
    elseif length(f) == NT  % piecewise constant
           sf = f.*area;
           F = accumarray([elem(:,1);elem(:,2);elem(:,3)],[sf;sf;sf]/3,[N,1]);
    else
        error('The length of f should be either number of elements or number of nodes.')
    end
else         % f is a function
    % compute quadrature points
    [lambda,weight] = getquadpts(qf);
    phi = lambda;
    nQuad = size(lambda,1);
    % compute element-wise integral
    ft = zeros(NT,3);
    for j = 1:3
        for p = 1:nQuad
            % quadrature points in the x-y coordinate
            pxy = lambda(p,1)*node(elem(:,1),:) ...
                + lambda(p,2)*node(elem(:,2),:) ...
                + lambda(p,3)*node(elem(:,3),:);
            ft(:,j) = ft(:,j) + weight(p)*f(pxy)*phi(p,j);
        end
        ft(:,j) = ft(:,j).*area;
    end
    % assemble to node-wise quantity
    F = accumarray(elem(:),ft(:),[N 1]);
end


Published with MATLAB® 7.9