quadpts

QUAPTS quadrature points
Order 1, nQuad 1
Order 2, nQuad 3
Order 3, nQuad 4
Order 4, nQuad 6
Order 5, nQuad 7
Verification

function [lambda,weight] = quadpts(order)

QUAPTS quadrature points

[lambda,weight] = quadpts(order) return quadrature points with given order (up to 5) in the barycentric coordinates.

The output lambda is a matrix of size nQ by 3, where nQ is the number of quadrature points. lambda(i,:) is three barycentric coordinate of the i-th quadrature point and lambda(:,j) is the j-th barycentric coordinate of all quadrature points. The x-y coordinate of the p-th quadrature point can be computed as

   pxy = lambda(p,1)*node(elem(:,1),:) ...
       + lambda(p,2)*node(elem(:,2),:) ...
       + lambda(p,3)*node(elem(:,3),:);

The weight of p-th quadrature point is given by weight(p). See verifyquadpts for the usage of qudrature rules to compute integrals over triangles.

References:

David Dunavant. High degree efficient symmetrical Gaussian quadrature rules for the triangle. International journal for numerical methods in engineering. 21(6):1129--1148, 1985.
John Burkardt. DUNAVANT Quadrature Rules for the Triangle. http://people.sc.fsu.edu/~burkardt/m_src/dunavant/dunavant.html

See also verifyquadpts, getmatrix, dunavant/dunavant_subrule

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

if order>5
    order = 5;
end
switch order

Order 1, nQuad 1

    case 1

        lambda = [1/3, 1/3, 1/3];
        weight = 1;

Order 2, nQuad 3

    case 2

        lambda = [2/3, 1/6, 1/6; ...
                  1/6, 2/3, 1/6; ...
                  1/6, 1/6, 2/3];
        weight = [1/3, 1/3, 1/3];

Order 3, nQuad 4

    case 3

        lambda = [1/3, 1/3, 1/3; ...
                  0.6, 0.2, 0.2; ...
                  0.2, 0.6, 0.2; ...
                  0.2, 0.2, 0.6];
        weight = [-27/48, 25/48, 25/48, 25/48];

Order 4, nQuad 6

    case 4

        lambda = [0.108103018168070, 0.445948490915965, 0.445948490915965; ...
                  0.445948490915965, 0.108103018168070, 0.445948490915965; ...
                  0.445948490915965, 0.445948490915965, 0.108103018168070; ...
                  0.816847572980459, 0.091576213509771, 0.091576213509771; ...
                  0.091576213509771, 0.816847572980459, 0.091576213509771; ...
                  0.091576213509771, 0.091576213509771, 0.816847572980459];
        weight = [0.223381589678011, 0.223381589678011, 0.223381589678011, ...
                  0.109951743655322, 0.109951743655322, 0.109951743655322];

Order 5, nQuad 7

    case 5

        alpha1 = 0.059715871789770;      beta1 = 0.470142064105115;
        alpha2 = 0.797426985353087;      beta2 = 0.101286507323456;
        lambda = [   1/3,    1/3,    1/3; ...
                  alpha1,  beta1,  beta1; ...
                   beta1, alpha1,  beta1; ...
                   beta1,  beta1, alpha1; ...
                  alpha2,  beta2,  beta2; ...
                   beta2, alpha2,  beta2; ...
                   beta2,  beta2, alpha2];
        weight = [0.225, 0.132394152788506, 0.132394152788506, 0.132394152788506, ...
                         0.125939180544827, 0.125939180544827, 0.125939180544827];
end

Verification

The order of the quadrature rule is verified by the function verifyquadpts. See verifyquadpts.