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_MONOMIAL_H 21 #define LIBMESH_QUADRATURE_MONOMIAL_H 22 23 // Local includes 24 #include "libmesh/quadrature.h" 25 26 namespace libMesh 27 { 28 29 /** 30 * This class defines alternate quadrature rules on "tensor-product" 31 * elements (quadrilaterals and hexahedra) which can be useful when 32 * integrating monomial finite element bases. 33 * 34 * While tensor product rules are optimal for integrating bi/tri-linear, 35 * bi/tri-quadratic, etc. (i.e. tensor product) bases (which consist 36 * of incomplete polynomials up to degree=dim*p) they are not optimal 37 * for the MONOMIAL or FEXYZ bases, which consist of complete 38 * polynomials of degree=p. 39 * 40 * This class provides quadrature rules which are more efficient than 41 * tensor product rules when they are available, and falls back on 42 * Gaussian quadrature rules otherwise. 43 * 44 * A number of these rules have been helpfully collected in electronic form by: 45 * Prof. Ronald Cools 46 * Katholieke Universiteit Leuven, 47 * Dept. Computerwetenschappen 48 * http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 49 * A username and password to access the tables is available by request. 50 * 51 * We also provide the original reference for each rule when it is 52 * available. 53 * 54 * \author John W. Peterson 55 * \date 2008 56 * \brief Implements quadrature rules for non-tensor polynomials. 57 */ 58 class QMonomial final : public QBase 59 { 60 public: 61 62 /** 63 * Constructor. Declares the order of the quadrature rule. 64 */ 65 QMonomial (unsigned int dim, 66 Order order=INVALID_ORDER) : QBase(dim,order)67 QBase(dim,order) 68 {} 69 70 /** 71 * Copy/move ctor, copy/move assignment operator, and destructor are 72 * all explicitly defaulted for this simple class. 73 */ 74 QMonomial (const QMonomial &) = default; 75 QMonomial (QMonomial &&) = default; 76 QMonomial & operator= (const QMonomial &) = default; 77 QMonomial & operator= (QMonomial &&) = default; 78 virtual ~QMonomial() = default; 79 80 /** 81 * \returns \p QMONOMIAL. 82 */ 83 virtual QuadratureType type() const override; 84 85 86 private: 87 88 /** 89 * Uses a Gauss rule in 1D. More efficient rules for non tensor 90 * product bases on quadrilaterals and hexahedra. 91 */ 92 virtual void init_1D (const ElemType, unsigned int) override; 93 virtual void init_2D (const ElemType, unsigned int) override; 94 virtual void init_3D (const ElemType, unsigned int) override; 95 96 /** 97 * Wissmann published three interesting "partially symmetric" rules 98 * for integrating degree 4, 6, and 8 polynomials exactly on QUADs. 99 * These rules have all positive weights, all points inside the 100 * reference element, and have fewer points than tensor-product 101 * rules of equivalent order, making them superior to those rules 102 * for monomial bases. 103 * 104 * J. W. Wissman and T. Becker, Partially symmetric cubature 105 * formulas for even degrees of exactness, SIAM J. Numer. Anal. 23 106 * (1986), 676--685. 107 */ 108 void wissmann_rule(const Real rule_data[][3], 109 const unsigned int n_pts); 110 111 /** 112 * Stroud's rules for quads and hexes can have one of several 113 * different types of symmetry. The rule_symmetry array describes 114 * how the different lines of the rule_data array are to be 115 * applied. The different rule_symmetry possibilities are: 116 * 0) Origin or single-point: (x,y) 117 * Fully-symmetric, 3 cases: 118 * 1) (x,y) -> (x,y), (-x,y), (x,-y), (-x,-y) 119 * (y,x), (-y,x), (y,-x), (-y,-x) 120 * 2) (x,x) -> (x,x), (-x,x), (x,-x), (-x,-x) 121 * 3) (x,0) -> (x,0), (-x,0), (0, x), ( 0,-x) 122 * 4) Rotational Invariant, (x,y) -> (x,y), (-x,-y), (-y, x), (y,-x) 123 * 5) Partial Symmetry, (x,y) -> (x,y), (-x, y) [x!=0] 124 * 6) Rectangular Symmetry, (x,y) -> (x,y), (-x, y), (-x,-y), (x,-y) 125 * 7) Central Symmetry, (0,y) -> (0,y), ( 0,-y) 126 * 127 * Not all rules with these symmetries are due to Stroud, however, 128 * his book is probably the most frequently-cited compendium of 129 * quadrature rules and later authors certainly built upon his work. 130 */ 131 void stroud_rule(const Real rule_data[][3], 132 const unsigned int * rule_symmetry, 133 const unsigned int n_pts); 134 135 /** 136 * Rules from Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931. 137 * The rules are obtained by considering the group G^{rot} of rotations of the reference 138 * hex, and the invariant polynomials of this group. 139 * 140 * In Kim and Song's rules, quadrature points are described by the following points 141 * and their unique permutations under the G^{rot} group: 142 * 143 * 0.) (0,0,0) ( 1 perm ) -> [0, 0, 0] 144 * 1.) (x,0,0) ( 6 perms) -> [x, 0, 0], [0, -x, 0], [-x, 0, 0], [0, x, 0], [0, 0, -x], [0, 0, x] 145 * 2.) (x,x,0) (12 perms) -> [x, x, 0], [x, -x, 0], [-x, -x, 0], [-x, x, 0], [x, 0, -x], [x, 0, x], 146 * [0, x, -x], [0, x, x], [0, -x, -x], [-x, 0, -x], [0, -x, x], [-x, 0, x] 147 * 3.) (x,y,0) (24 perms) -> [x, y, 0], [y, -x, 0], [-x, -y, 0], [-y, x, 0], [x, 0, -y], [x, -y, 0], 148 * [x, 0, y], [0, y, -x], [-x, y, 0], 149 * [0, y, x], [y, 0, -x], [0, -y, -x], [-y, 0, -x], [y, x, 0], [-y, -x, 0], 150 * [y, 0, x], [0, -y, x], [-y, 0, x], 151 * [-x, 0, y], [0, -x, -y], [0, -x, y], [-x, 0, -y], [0, x, y], [0, x, -y] 152 * 4.) (x,x,x) ( 8 perms) -> [x, x, x], [x, -x, x], [-x, -x, x], [-x, x, x], [x, x, -x], [x, -x, -x], 153 * [-x, x, -x], [-x, -x, -x] 154 * 5.) (x,x,z) (24 perms) -> [x, x, z], [x, -x, z], [-x, -x, z], [-x, x, z], [x, z, -x], [x, -x, -z], 155 * [x, -z, x], [z, x, -x], [-x, x, -z], 156 * [-z, x, x], [x, -z, -x], [-z, -x, -x], [-x, z, -x], [x, x, -z], 157 * [-x, -x, -z], [x, z, x], [z, -x, x], 158 * [-x, -z, x], [-x, z, x], [z, -x, -x], [-z, -x, x], [-x, -z, -x], 159 * [z, x, x], [-z, x, -x] 160 * 6.) (x,y,z) (24 perms) -> [x, y, z], [y, -x, z], [-x, -y, z], [-y, x, z], [x, z, -y], [x, -y, -z], 161 * [x, -z, y], [z, y, -x], [-x, y, -z], 162 * [-z, y, x], [y, -z, -x], [-z, -y, -x], [-y, z, -x], [y, x, -z], [-y, -x, -z], 163 * [y, z, x], [z, -y, x], [-y, -z, x], [-x, z, y], [z, -x, -y], 164 * [-z, -x, y], [-x, -z, -y], [z, x, y], [-z, x, -y] 165 * 166 * Only two of Kim and Song's rules are particularly useful for FEM 167 * calculations: the degree 7, 38-point rule and their degree 8, 168 * 47-point rule. The others either contain negative weights or 169 * points outside the reference interval. The points and weights, 170 * to 32 digits, were obtained from: Ronald Cools' website 171 * (http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html) and 172 * the unique permutations of G^{rot} were computed by me [JWP] 173 * using Maple. 174 */ 175 void kim_rule(const Real rule_data[][4], 176 const unsigned int * rule_id, 177 const unsigned int n_pts); 178 }; 179 180 } // namespace libMesh 181 182 #endif // LIBMESH_QUADRATURE_MONOMIAL_H 183