Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : #ifndef LIBMESH_QUADRATURE_GM_H
      21             : #define LIBMESH_QUADRATURE_GM_H
      22             : 
      23             : // Local includes
      24             : #include "libmesh/quadrature.h"
      25             : 
      26             : namespace libMesh
      27             : {
      28             : 
      29             : /**
      30             :  * This class implements the Grundmann-Moller quadrature rules for
      31             :  * tetrahedra.  The GM rules are well-defined for simplices of
      32             :  * arbitrary dimension and to any order, but the rules by Dunavant for
      33             :  * two-dimensional simplices are in general superior.  This is primarily
      34             :  * due to the fact that the GM rules contain a significant proportion
      35             :  * of negative weights, making them susceptible to round-off error
      36             :  * at high-order.
      37             :  *
      38             :  * The GM rules are interesting in 3D because they overlap with the
      39             :  * conical product rules at higher order while having significantly
      40             :  * fewer evaluation points, making them potentially much more
      41             :  * efficient.  The table below gives a comparison between the number
      42             :  * of points in a conical product (CP) rule and the GM rule of
      43             :  * equivalent order.  The GM rules are defined to be exact for
      44             :  * polynomials of degree d=2*s+1, s=0,1,2,3,... The table also gives
      45             :  * the percentage of each GM rule's weights which are negative.
      46             :  * The percentage of negative weights appears to approach 50, and the
      47             :  * amplification factor (a measure of the effect of round-off) defined
      48             :  * as
      49             :  *
      50             :  * amp. factor = \f$ \frac{1}{V} \sum \|w_i\|, \f$
      51             :  *
      52             :  * where V is the volume of the reference element, grows like exp(C*s).
      53             :  * (A rule with all positive weights has an amplification factor of
      54             :  * 1.0 by definition.)
      55             :  *
      56             :  * \verbatim
      57             :  * s   degree  n_pts(conical)  n_pts(GM)   % neg wts  amp. factor
      58             :  * ------------------------------------------------------------------------
      59             :  * 0   1       1               1            0.00      1.00e+00
      60             :  * 1   3       8               5           20.00      2.60e+00
      61             :  * 2   5       27              15          26.67      5.63e+00
      62             :  * 3   7       64              35          31.43      1.19e+01
      63             :  * 4   9       125             70          34.29      2.54e+01
      64             :  * 5   11      216             126         36.51      5.41e+01
      65             :  * 6   13      343             210         38.10      1.16e+02
      66             :  * 7   15      512             330         39.39      2.51e+02
      67             :  * 8   17      729             495         40.40      5.45e+02
      68             :  * 9   19      1000            715         41.26      1.19e+03
      69             :  * 10  21      1331            1001        41.96      2.59e+03
      70             :  * 11  23      1728            1365        42.56      5.68e+03
      71             :  * 12  25      2197            1820        43.08      1.25e+04
      72             :  * 13  27      2744            2380        43.53      2.75e+04
      73             :  * 14  29      3375            3060        43.92      6.07e+04
      74             :  * 15  31      4096            3876        44.27      1.34e+05
      75             :  * 16  33      4913            4845        44.58      2.97e+05
      76             :  * 17  35      5832            5985        44.86      6.59e+05 <= Conical rule has fewer points for degree >= 34
      77             :  * 18  37      6859            7315        45.11      1.46e+06
      78             :  * 19  39      8000            8855        45.34      3.25e+06
      79             :  * 20  41      9261            10626       45.55      7.23e+06
      80             :  * 21  43      10648           12650       45.74      1.61e+07
      81             :  * \endverbatim
      82             :  *
      83             :  * Reference:
      84             :  *    Axel Grundmann and Michael M\"{o}ller, "Invariant Integration
      85             :  *    Formulas for the N-Simplex by Combinatorial Methods," SIAM
      86             :  *    Journal on Numerical Analysis, Volume 15, Number 2, April 1978,
      87             :  *    pages 282-290.
      88             :  *
      89             :  * Reference LGPL Fortran90 code by John Burkardt can be found here:
      90             :  * http://people.scs.fsu.edu/~burkardt/f_src/gm_rules/gm_rules.html
      91             :  *
      92             :  * \author John W. Peterson
      93             :  * \date 2008
      94             :  * \brief Implements the quadrature rules of Grundmann and Moller in 2D and 3D.
      95             :  */
      96             : class QGrundmann_Moller final : public QBase
      97             : {
      98             : public:
      99             : 
     100             :   /**
     101             :    * Constructor.  Declares the order of the quadrature rule.
     102             :    */
     103        1467 :   QGrundmann_Moller (unsigned int dim,
     104        7317 :                      Order order=INVALID_ORDER) :
     105        7317 :     QBase(dim,order)
     106             :   {
     107        1467 :     if (dim == 1)
     108           0 :       init(EDGE2);
     109        1467 :   }
     110             : 
     111             :   /**
     112             :    * Copy/move ctor, copy/move assignment operator, and destructor are
     113             :    * all explicitly defaulted for this simple class.
     114             :    */
     115           0 :   QGrundmann_Moller (const QGrundmann_Moller &) = default;
     116             :   QGrundmann_Moller (QGrundmann_Moller &&) = default;
     117             :   QGrundmann_Moller & operator= (const QGrundmann_Moller &) = default;
     118             :   QGrundmann_Moller & operator= (QGrundmann_Moller &&) = default;
     119        7351 :   virtual ~QGrundmann_Moller() = default;
     120             : 
     121             :   /**
     122             :    * \returns \p QGRUNDMANN_MOLLER.
     123             :    */
     124             :   virtual QuadratureType type() const override;
     125             : 
     126             :   virtual std::unique_ptr<QBase> clone() const override;
     127             : 
     128             : private:
     129             : 
     130             :   /**
     131             :    * In 1D, use a Gauss rule.
     132             :    * In 2D, the GM product rule is only defined for Tris.
     133             :    * In 3D, the GM product rule is only defined for Tets.
     134             :    */
     135             :   virtual void init_1D () override;
     136             :   virtual void init_2D () override;
     137             :   virtual void init_3D () override;
     138             : 
     139             :   /**
     140             :    * This routine is called from init_2D() and init_3D().  It actually
     141             :    * fills the _points and _weights vectors for a given rule index, s
     142             :    * and dimension, dim.
     143             :    */
     144             :   void gm_rule(unsigned int s, unsigned int dim);
     145             : 
     146             :   /**
     147             :    * Routine which generates p-compositions of a given order, s,
     148             :    * as well as permutations thereof.  This routine is called internally by
     149             :    * the gm_rule() routine, you should not call this yourself!
     150             :    */
     151             :   void compose_all(unsigned int s, // number to be compositioned
     152             :                    unsigned int p, // # of partitions
     153             :                    std::vector<std::vector<unsigned int>> & result);
     154             : };
     155             : 
     156             : } // namespace libMesh
     157             : 
     158             : #endif // LIBMESH_QUADRATURE_GM_H