
<!DOCTYPE html
  PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   <!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
      --><title>getquadpts</title><meta name="generator" content="MATLAB 7.9"><meta name="date" content="2009-12-25"><meta name="m-file" content="getquadpts"><style type="text/css">

body {
  background-color: white;
  margin:10px;
}

h1 {
  color: #990000; 
  font-size: x-large;
}

h2 {
  color: #990000;
  font-size: medium;
}

/* Make the text shrink to fit narrow windows, but not stretch too far in 
wide windows. */ 
p,h1,h2,div.content div {
  max-width: 600px;
  /* Hack for IE6 */
  width: auto !important; width: 600px;
}

pre.codeinput {
  background: #EEEEEE;
  padding: 10px;
}
@media print {
  pre.codeinput {word-wrap:break-word; width:100%;}
} 

span.keyword {color: #0000FF}
span.comment {color: #228B22}
span.string {color: #A020F0}
span.untermstring {color: #B20000}
span.syscmd {color: #B28C00}

pre.codeoutput {
  color: #666666;
  padding: 10px;
}

pre.error {
  color: red;
}

p.footer {
  text-align: right;
  font-size: xx-small;
  font-weight: lighter;
  font-style: italic;
  color: gray;
}

  </style></head><body><div class="content"><h2>Contents</h2><div><ul><li><a href="#2">GETQUAPTS quadrature points</a></li><li><a href="#3">Order 1, nQuad 1</a></li><li><a href="#5">Order 2, nQuad 3</a></li><li><a href="#8">Order 3, nQuad 4</a></li><li><a href="#11">Order 4, nQuad 6</a></li><li><a href="#14">Order 5, nQuad 7</a></li><li><a href="#17">Verification</a></li></ul></div><pre class="codeinput"><span class="keyword">function</span> [lambda,weight] = getquadpts(order)
</pre><h2>GETQUAPTS quadrature points<a name="2"></a></h2><p>[lambda,weight] = getquadpts(order) return quadrature points with given order (up to 5) in the barycentric coordinates.</p><p>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</p><pre>   pxy = lambda(p,1)*node(elem(:,1),:) ...
       + lambda(p,2)*node(elem(:,2),:) ...
       + lambda(p,3)*node(elem(:,3),:);</pre><p>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.</p><p>References:</p><div><ul><li>David Dunavant. High degree efficient symmetrical Gaussian    quadrature rules for the triangle. International journal for numerical    methods in engineering. 21(6):1129--1148, 1985.</li><li>John Burkardt. DUNAVANT Quadrature Rules for the Triangle.    <a href="http://people.sc.fsu.edu/~burkardt/m_src/dunavant/dunavant.html">http://people.sc.fsu.edu/~burkardt/m_src/dunavant/dunavant.html</a></li></ul></div><p>See also verifyquadpts, getmatrix, dunavant/dunavant_subrule</p><pre class="codeinput"><span class="comment">% Copyright (C) Long Chen. See COPYRIGHT.txt for details.</span>

<span class="keyword">if</span> order&gt;5
    order = 5;
<span class="keyword">end</span>
<span class="keyword">switch</span> order
</pre><h2>Order 1, nQuad 1<a name="3"></a></h2><pre class="codeinput">    <span class="keyword">case</span> 1
</pre><pre class="codeinput">        lambda = [1/3, 1/3, 1/3];
        weight = 1;
</pre><h2>Order 2, nQuad 3<a name="5"></a></h2><pre class="codeinput">    <span class="keyword">case</span> 2
</pre><pre class="codeinput">        lambda = [2/3, 1/6, 1/6; <span class="keyword">...</span>
                  1/6, 2/3, 1/6; <span class="keyword">...</span>
                  1/6, 1/6, 2/3];
        weight = [1/3, 1/3, 1/3];
</pre><h2>Order 3, nQuad 4<a name="8"></a></h2><pre class="codeinput">    <span class="keyword">case</span> 3
</pre><pre class="codeinput">        lambda = [1/3, 1/3, 1/3; <span class="keyword">...</span>
                  0.6, 0.2, 0.2; <span class="keyword">...</span>
                  0.2, 0.6, 0.2; <span class="keyword">...</span>
                  0.2, 0.2, 0.6];
        weight = [-27/48, 25/48, 25/48, 25/48];
</pre><h2>Order 4, nQuad 6<a name="11"></a></h2><pre class="codeinput">    <span class="keyword">case</span> 4
</pre><pre class="codeinput">        lambda = [0.108103018168070, 0.445948490915965, 0.445948490915965; <span class="keyword">...</span>
                  0.445948490915965, 0.108103018168070, 0.445948490915965; <span class="keyword">...</span>
                  0.445948490915965, 0.445948490915965, 0.108103018168070; <span class="keyword">...</span>
                  0.816847572980459, 0.091576213509771, 0.091576213509771; <span class="keyword">...</span>
                  0.091576213509771, 0.816847572980459, 0.091576213509771; <span class="keyword">...</span>
                  0.091576213509771, 0.091576213509771, 0.816847572980459];
        weight = [0.223381589678011, 0.223381589678011, 0.223381589678011, <span class="keyword">...</span>
                  0.109951743655322, 0.109951743655322, 0.109951743655322];
</pre><h2>Order 5, nQuad 7<a name="14"></a></h2><pre class="codeinput">    <span class="keyword">case</span> 5
</pre><pre class="codeinput">        alpha1 = 0.059715871789770;      beta1 = 0.470142064105115;
        alpha2 = 0.797426985353087;      beta2 = 0.101286507323456;
        lambda = [   1/3,    1/3,    1/3; <span class="keyword">...</span>
                  alpha1,  beta1,  beta1; <span class="keyword">...</span>
                   beta1, alpha1,  beta1; <span class="keyword">...</span>
                   beta1,  beta1, alpha1; <span class="keyword">...</span>
                  alpha2,  beta2,  beta2; <span class="keyword">...</span>
                   beta2, alpha2,  beta2; <span class="keyword">...</span>
                   beta2,  beta2, alpha2];
        weight = [0.225, 0.132394152788506, 0.132394152788506, 0.132394152788506, <span class="keyword">...</span>
                         0.125939180544827, 0.125939180544827, 0.125939180544827];
<span class="keyword">end</span>
</pre><h2>Verification<a name="17"></a></h2><p>The order of the quadrature rule is verified by the function verifyquadpts. See <a href="matlab:ifemdoc('verifyquadpts')">verifyquadpts</a>.</p><p class="footer"><br>
      Published with MATLAB&reg; 7.9<br></p></div><!--
##### SOURCE BEGIN #####
function [lambda,weight] = getquadpts(order)
%% GETQUAPTS quadrature points
%
% [lambda,weight] = getquadpts(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):1129REPLACE_WITH_DASH_DASH1148, 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 <matlab:ifemdoc('verifyquadpts') verifyquadpts>.
##### SOURCE END #####
--></body></html>