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_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 8449 : QMonomial (unsigned int dim,
66 8449 : Order order=INVALID_ORDER) :
67 8449 : QBase(dim,order)
68 : {
69 8449 : if (dim == 1)
70 1207 : init(EDGE2);
71 8449 : }
72 :
73 : /**
74 : * Copy/move ctor, copy/move assignment operator, and destructor are
75 : * all explicitly defaulted for this simple class.
76 : */
77 0 : QMonomial (const QMonomial &) = default;
78 : QMonomial (QMonomial &&) = default;
79 : QMonomial & operator= (const QMonomial &) = default;
80 : QMonomial & operator= (QMonomial &&) = default;
81 8687 : virtual ~QMonomial() = default;
82 :
83 : /**
84 : * \returns \p QMONOMIAL.
85 : */
86 : virtual QuadratureType type() const override;
87 :
88 : virtual std::unique_ptr<QBase> clone() const override;
89 :
90 : private:
91 :
92 : /**
93 : * Uses a Gauss rule in 1D. More efficient rules for non tensor
94 : * product bases on quadrilaterals and hexahedra.
95 : */
96 : virtual void init_1D () override;
97 : virtual void init_2D () override;
98 : virtual void init_3D () override;
99 :
100 : /**
101 : * Wissmann published three interesting "partially symmetric" rules
102 : * for integrating degree 4, 6, and 8 polynomials exactly on QUADs.
103 : * These rules have all positive weights, all points inside the
104 : * reference element, and have fewer points than tensor-product
105 : * rules of equivalent order, making them superior to those rules
106 : * for monomial bases.
107 : *
108 : * J. W. Wissman and T. Becker, Partially symmetric cubature
109 : * formulas for even degrees of exactness, SIAM J. Numer. Anal. 23
110 : * (1986), 676--685.
111 : */
112 : void wissmann_rule(const Real rule_data[][3],
113 : const unsigned int n_pts);
114 :
115 : /**
116 : * Stroud's rules for quads and hexes can have one of several
117 : * different types of symmetry. The rule_symmetry array describes
118 : * how the different lines of the rule_data array are to be
119 : * applied. The different rule_symmetry possibilities are:
120 : * 0) Origin or single-point: (x,y)
121 : * Fully-symmetric, 3 cases:
122 : * 1) (x,y) -> (x,y), (-x,y), (x,-y), (-x,-y)
123 : * (y,x), (-y,x), (y,-x), (-y,-x)
124 : * 2) (x,x) -> (x,x), (-x,x), (x,-x), (-x,-x)
125 : * 3) (x,0) -> (x,0), (-x,0), (0, x), ( 0,-x)
126 : * 4) Rotational Invariant, (x,y) -> (x,y), (-x,-y), (-y, x), (y,-x)
127 : * 5) Partial Symmetry, (x,y) -> (x,y), (-x, y) [x!=0]
128 : * 6) Rectangular Symmetry, (x,y) -> (x,y), (-x, y), (-x,-y), (x,-y)
129 : * 7) Central Symmetry, (0,y) -> (0,y), ( 0,-y)
130 : *
131 : * Not all rules with these symmetries are due to Stroud, however,
132 : * his book is probably the most frequently-cited compendium of
133 : * quadrature rules and later authors certainly built upon his work.
134 : */
135 : void stroud_rule(const Real rule_data[][3],
136 : const unsigned int * rule_symmetry,
137 : const unsigned int n_pts);
138 :
139 : /**
140 : * Rules from Kim and Song, Comm. Korean Math. Soc vol. 13, no. 4, 1998, pp. 913-931.
141 : * The rules are obtained by considering the group G^{rot} of rotations of the reference
142 : * hex, and the invariant polynomials of this group.
143 : *
144 : * In Kim and Song's rules, quadrature points are described by the following points
145 : * and their unique permutations under the G^{rot} group:
146 : *
147 : * 0.) (0,0,0) ( 1 perm ) -> [0, 0, 0]
148 : * 1.) (x,0,0) ( 6 perms) -> [x, 0, 0], [0, -x, 0], [-x, 0, 0], [0, x, 0], [0, 0, -x], [0, 0, x]
149 : * 2.) (x,x,0) (12 perms) -> [x, x, 0], [x, -x, 0], [-x, -x, 0], [-x, x, 0], [x, 0, -x], [x, 0, x],
150 : * [0, x, -x], [0, x, x], [0, -x, -x], [-x, 0, -x], [0, -x, x], [-x, 0, x]
151 : * 3.) (x,y,0) (24 perms) -> [x, y, 0], [y, -x, 0], [-x, -y, 0], [-y, x, 0], [x, 0, -y], [x, -y, 0],
152 : * [x, 0, y], [0, y, -x], [-x, y, 0],
153 : * [0, y, x], [y, 0, -x], [0, -y, -x], [-y, 0, -x], [y, x, 0], [-y, -x, 0],
154 : * [y, 0, x], [0, -y, x], [-y, 0, x],
155 : * [-x, 0, y], [0, -x, -y], [0, -x, y], [-x, 0, -y], [0, x, y], [0, x, -y]
156 : * 4.) (x,x,x) ( 8 perms) -> [x, x, x], [x, -x, x], [-x, -x, x], [-x, x, x], [x, x, -x], [x, -x, -x],
157 : * [-x, x, -x], [-x, -x, -x]
158 : * 5.) (x,x,z) (24 perms) -> [x, x, z], [x, -x, z], [-x, -x, z], [-x, x, z], [x, z, -x], [x, -x, -z],
159 : * [x, -z, x], [z, x, -x], [-x, x, -z],
160 : * [-z, x, x], [x, -z, -x], [-z, -x, -x], [-x, z, -x], [x, x, -z],
161 : * [-x, -x, -z], [x, z, x], [z, -x, x],
162 : * [-x, -z, x], [-x, z, x], [z, -x, -x], [-z, -x, x], [-x, -z, -x],
163 : * [z, x, x], [-z, x, -x]
164 : * 6.) (x,y,z) (24 perms) -> [x, y, z], [y, -x, z], [-x, -y, z], [-y, x, z], [x, z, -y], [x, -y, -z],
165 : * [x, -z, y], [z, y, -x], [-x, y, -z],
166 : * [-z, y, x], [y, -z, -x], [-z, -y, -x], [-y, z, -x], [y, x, -z], [-y, -x, -z],
167 : * [y, z, x], [z, -y, x], [-y, -z, x], [-x, z, y], [z, -x, -y],
168 : * [-z, -x, y], [-x, -z, -y], [z, x, y], [-z, x, -y]
169 : *
170 : * Only two of Kim and Song's rules are particularly useful for FEM
171 : * calculations: the degree 7, 38-point rule and their degree 8,
172 : * 47-point rule. The others either contain negative weights or
173 : * points outside the reference interval. The points and weights,
174 : * to 32 digits, were obtained from: Ronald Cools' website
175 : * (http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html) and
176 : * the unique permutations of G^{rot} were computed by me [JWP]
177 : * using Maple.
178 : */
179 : void kim_rule(const Real rule_data[][4],
180 : const unsigned int * rule_id,
181 : const unsigned int n_pts);
182 : };
183 :
184 : } // namespace libMesh
185 :
186 : #endif // LIBMESH_QUADRATURE_MONOMIAL_H
|
