1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 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 QGrundmann_Moller (unsigned int dim, 104 Order order=INVALID_ORDER) : QBase(dim,order)105 QBase(dim,order) 106 {} 107 108 /** 109 * Copy/move ctor, copy/move assignment operator, and destructor are 110 * all explicitly defaulted for this simple class. 111 */ 112 QGrundmann_Moller (const QGrundmann_Moller &) = default; 113 QGrundmann_Moller (QGrundmann_Moller &&) = default; 114 QGrundmann_Moller & operator= (const QGrundmann_Moller &) = default; 115 QGrundmann_Moller & operator= (QGrundmann_Moller &&) = default; 116 virtual ~QGrundmann_Moller() = default; 117 118 /** 119 * \returns \p QGRUNDMANN_MOLLER. 120 */ 121 virtual QuadratureType type() const override; 122 123 124 private: 125 126 /** 127 * In 1D, use a Gauss rule. 128 * In 2D, the GM product rule is only defined for Tris. 129 * In 3D, the GM product rule is only defined for Tets. 130 */ 131 virtual void init_1D (const ElemType, unsigned int) override; 132 virtual void init_2D (const ElemType, unsigned int) override; 133 virtual void init_3D (const ElemType, unsigned int) override; 134 135 /** 136 * This routine is called from init_2D() and init_3D(). It actually 137 * fills the _points and _weights vectors for a given rule index, s 138 * and dimension, dim. 139 */ 140 void gm_rule(unsigned int s, unsigned int dim); 141 142 /** 143 * Routine which generates p-compositions of a given order, s, 144 * as well as permutations thereof. This routine is called internally by 145 * the gm_rule() routine, you should not call this yourself! 146 */ 147 void compose_all(unsigned int s, // number to be compositioned 148 unsigned int p, // # of partitions 149 std::vector<std::vector<unsigned int>> & result); 150 }; 151 152 } // namespace libMesh 153 154 #endif // LIBMESH_QUADRATURE_GM_H 155