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_H
21 : #define LIBMESH_QUADRATURE_H
22 :
23 : // Local includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/reference_counted_object.h"
26 : #include "libmesh/point.h"
27 : #include "libmesh/enum_elem_type.h" // INVALID_ELEM
28 : #include "libmesh/enum_order.h" // INVALID_ORDER
29 :
30 : // C++ includes
31 : #include <vector>
32 : #include <string>
33 : #include <utility>
34 : #include <memory>
35 :
36 : // We used to have additional arguments to the init_0D, init_1D, etc.
37 : // classes, but they became unused after QBase members made them
38 : // redundant, so they were deprecated and eventually removed. Users'
39 : // derived classes that need to compile against both old and new
40 : // libMesh version can test this macro to determine what signature
41 : // their overrides should have.
42 : #define LIBMESH_QBASE_INIT_ARGUMENTS_REMOVED
43 :
44 : namespace libMesh
45 : {
46 :
47 : // forward declarations
48 : class Elem;
49 : enum QuadratureType : int;
50 :
51 : /**
52 : * The \p QBase class provides the basic functionality from which
53 : * various quadrature rules can be derived. It computes and stores
54 : * the quadrature points (in reference element space) and associated
55 : * weights.
56 : *
57 : * \author Benjamin S. Kirk
58 : * \date 2002
59 : * \brief Base class for all quadrature families and orders.
60 : */
61 : class QBase : public ReferenceCountedObject<QBase>
62 : {
63 : protected:
64 :
65 : /**
66 : * Constructor. Protected to prevent instantiation of this base
67 : * class. Use the build() method instead.
68 : */
69 : QBase (unsigned int dim,
70 : Order order=INVALID_ORDER);
71 :
72 : public:
73 :
74 : /**
75 : * Copy/move ctor, copy/move assignment operator, and destructor are
76 : * all explicitly defaulted for this simple class.
77 : */
78 0 : QBase (const QBase &) = default;
79 : QBase (QBase &&) = default;
80 : QBase & operator= (const QBase &) = default;
81 : QBase & operator= (QBase &&) = default;
82 4002233 : virtual ~QBase() = default;
83 :
84 : /**
85 : * \returns A copy of this quadrature rule wrapped in a smart pointer.
86 : */
87 : virtual std::unique_ptr<QBase> clone() const;
88 :
89 : /**
90 : * \returns The quadrature type in derived classes.
91 : */
92 : virtual QuadratureType type() const = 0;
93 :
94 : /**
95 : * Builds a specific quadrature rule based on the \p name
96 : * string. This enables selection of the quadrature rule at
97 : * run-time. The input parameter \p name must be mappable through
98 : * the \p Utility::string_to_enum<>() function.
99 : *
100 : * This function allocates memory, therefore a \p std::unique_ptr<QBase>
101 : * is returned so that the user does not accidentally leak it.
102 : */
103 : static std::unique_ptr<QBase> build (std::string_view name,
104 : const unsigned int dim,
105 : const Order order=INVALID_ORDER);
106 :
107 : /**
108 : * Builds a specific quadrature rule based on the QuadratureType.
109 : * This enables selection of the quadrature rule at run-time.
110 : *
111 : * This function allocates memory, therefore a \p std::unique_ptr<QBase>
112 : * is returned so that the user does not accidentally leak it.
113 : */
114 : static std::unique_ptr<QBase> build (const QuadratureType qt,
115 : const unsigned int dim,
116 : const Order order=INVALID_ORDER);
117 :
118 : /**
119 : * \returns The element type we're currently using.
120 : */
121 : ElemType get_elem_type() const { return _type; }
122 :
123 : /**
124 : * \returns The p-refinement level we're currently using.
125 : */
126 : unsigned int get_p_level() const { return _p_level; }
127 :
128 : /**
129 : * \returns The number of points associated with the quadrature rule.
130 : */
131 6039 : unsigned int n_points() const
132 : {
133 6039 : libmesh_assert (!_points.empty());
134 3490469 : return cast_int<unsigned int>(_points.size());
135 : }
136 :
137 : /**
138 : * Alias for n_points() to enable use in index_range
139 : *
140 : * \returns The number of points associated with the quadrature rule.
141 : */
142 : unsigned int size() const
143 : {
144 : return n_points();
145 : }
146 :
147 : /**
148 : * \returns The spatial dimension of the quadrature rule.
149 : */
150 30 : unsigned int get_dim() const { return _dim; }
151 :
152 : /**
153 : * \returns A \p std::vector containing the quadrature point locations
154 : * in reference element space.
155 : */
156 758 : const std::vector<Point> & get_points() const { return _points; }
157 :
158 : /**
159 : * \returns A \p std::vector containing the quadrature point locations
160 : * in reference element space as a writable reference.
161 : */
162 121414216 : std::vector<Point> & get_points() { return _points; }
163 :
164 : /**
165 : * \returns A constant reference to a \p std::vector containing the
166 : * quadrature weights.
167 : */
168 164918 : const std::vector<Real> & get_weights() const { return _weights; }
169 :
170 : /**
171 : * \returns A writable references to a \p std::vector containing the
172 : * quadrature weights.
173 : */
174 145683261 : std::vector<Real> & get_weights() { return _weights; }
175 :
176 : /**
177 : * \returns The \f$ i^{th} \f$ quadrature point in reference element space.
178 : */
179 2210070 : Point qp(const unsigned int i) const
180 : {
181 2210070 : libmesh_assert_less (i, _points.size());
182 28610674 : return _points[i];
183 : }
184 :
185 : /**
186 : * \returns The \f$ i^{th} \f$ quadrature weight.
187 : */
188 2207257 : Real w(const unsigned int i) const
189 : {
190 2207257 : libmesh_assert_less (i, _weights.size());
191 26833724 : return _weights[i];
192 : }
193 :
194 : /**
195 : * Initializes the data structures for a quadrature rule for the
196 : * element \p e. If \p p_level is specified it overrides the element
197 : * p_level() elevation to use.
198 : */
199 : virtual void init (const Elem & e,
200 : unsigned int p_level=invalid_uint);
201 :
202 : /**
203 : * Initializes the data structures for a quadrature rule for an
204 : * element of type \p type. Some types, such as Polygon subclasses,
205 : * might require more detailed element information and so might not
206 : * be compatible with this API.
207 : *
208 : * New code should use the Elem-based API for most use cases, but
209 : * some code may initialize quadrature rules with simple ElemType
210 : * values like triangles and edges for use in tensor product or
211 : * conical product constructions; this code can set the
212 : * \p simple_type_only flag to avoid being identified as deprecated.
213 : */
214 : virtual void init (const ElemType type=INVALID_ELEM,
215 : unsigned int p_level=0,
216 : bool simple_type_only=false);
217 :
218 : /**
219 : * Initializes the data structures for a quadrature rule based on
220 : * the element, element type, and p_level settings of \p other_rule
221 : */
222 : virtual void init (const QBase & other_rule);
223 :
224 : /**
225 : * Initializes the data structures for an element potentially "cut"
226 : * by a signed distance function. The array \p vertex_distance_func
227 : * contains vertex values of the signed distance function. If the
228 : * signed distance function changes sign on the vertices, then the
229 : * element is considered to be cut.) This interface can be extended
230 : * by derived classes in order to subdivide the element and construct
231 : * a composite quadrature rule.
232 : */
233 : virtual void init (const Elem & elem,
234 : const std::vector<Real> & vertex_distance_func,
235 : unsigned int p_level=0);
236 :
237 : /**
238 : * \returns The current "total" order of the quadrature rule which
239 : * can vary element by element, depending on the Elem::p_level(),
240 : * which gets passed to us during init().
241 : *
242 : * Each additional power of p increases the quadrature order
243 : * required to integrate the mass matrix by 2, hence the formula
244 : * below.
245 : *
246 : * \todo This function should also be used in all of the Order
247 : * switch statements in the rules themselves.
248 : */
249 6300292 : Order get_order() const { return static_cast<Order>(_order + 2 * _p_level); }
250 :
251 : /**
252 : * \returns The "base" order of the quadrature rule, independent of
253 : * element.
254 : *
255 : * This function should be used when comparing quadrature objects
256 : * independently of their last initialization.
257 : */
258 : Order get_base_order() const { return static_cast<Order>(_order); }
259 :
260 : /**
261 : * Prints information relevant to the quadrature rule, by default to
262 : * libMesh::out.
263 : */
264 : void print_info(std::ostream & os=libMesh::out) const;
265 :
266 : /**
267 : * Maps the points of a 1D quadrature rule defined by "old_range" to
268 : * another 1D interval defined by "new_range" and scales the weights
269 : * accordingly.
270 : */
271 : void scale(std::pair<Real, Real> old_range,
272 : std::pair<Real, Real> new_range);
273 :
274 : /**
275 : * Same as above, but allows you to use the stream syntax.
276 : */
277 : friend std::ostream & operator << (std::ostream & os, const QBase & q);
278 :
279 : /**
280 : * \returns \p true if the shape functions need to be recalculated,
281 : * \p false otherwise.
282 : *
283 : * This may be required if the number of quadrature points or their
284 : * position changes.
285 : */
286 143912918 : virtual bool shapes_need_reinit() { return false; }
287 :
288 : /**
289 : * Flag (default true) controlling the use of quadrature rules with
290 : * negative weights. Set this to false to require rules with all
291 : * positive weights.
292 : *
293 : * Rules with negative weights can be unsuitable for some problems.
294 : * For example, it is possible for a rule with negative weights to
295 : * obtain a negative result when integrating a positive function.
296 : *
297 : * A particular example: if rules with negative weights are not allowed,
298 : * a request for TET,THIRD (5 points) will return the TET,FIFTH (14 points)
299 : * rule instead, nearly tripling the computational effort required!
300 : */
301 : bool allow_rules_with_negative_weights;
302 :
303 : /**
304 : * The flag's value defaults to false so that one does not accidentally
305 : * use a nodal quadrature rule on Pyramid elements, since evaluating the
306 : * inverse element Jacobian (e.g. dphi) is not well-defined at the
307 : * Pyramid apex because the element Jacobian is zero there.
308 : *
309 : * We do not want to completely prevent someone from using a nodal
310 : * quadrature rule on Pyramids, however, since there are legitimate
311 : * use cases (lumped mass matrix) so the flag can be set to true to
312 : * override this behavior.
313 : */
314 : bool allow_nodal_pyramid_quadrature;
315 :
316 : protected:
317 :
318 :
319 : /**
320 : * Initializes the 0D quadrature rule by filling the points and
321 : * weights vectors with the appropriate values. Generally this
322 : * is just one point with weight 1.
323 : */
324 : virtual void init_0D ();
325 :
326 : /**
327 : * Initializes the 1D quadrature rule by filling the points and
328 : * weights vectors with the appropriate values. The order of
329 : * the rule will be defined by the implementing class.
330 : * It is assumed that derived quadrature rules will at least
331 : * define the init_1D function, therefore it is pure virtual.
332 : */
333 : virtual void init_1D () = 0;
334 :
335 : /**
336 : * Initializes the 2D quadrature rule by filling the points and
337 : * weights vectors with the appropriate values. The order of
338 : * the rule will be defined by the implementing class.
339 : * Should not be pure virtual since a derived quadrature rule
340 : * may only be defined in 1D. If not overridden, throws an
341 : * error.
342 : */
343 : virtual void init_2D ();
344 :
345 : /**
346 : * Initializes the 3D quadrature rule by filling the points and
347 : * weights vectors with the appropriate values. The order of
348 : * the rule will be defined by the implementing class.
349 : * Should not be pure virtual since a derived quadrature rule
350 : * may only be defined in 1D. If not overridden, throws an
351 : * error.
352 : */
353 : virtual void init_3D ();
354 :
355 : /**
356 : * Constructs a 2D rule from the tensor product of \p q1D with
357 : * itself. Used in the \p init_2D() routines for quadrilateral
358 : * element types.
359 : */
360 : void tensor_product_quad (const QBase & q1D);
361 :
362 : /**
363 : * Computes the tensor product quadrature rule [q1D x q1D x q1D]
364 : * from the 1D rule q1D. Used in the init_3D routines for
365 : * hexahedral element types.
366 : */
367 : void tensor_product_hex (const QBase & q1D);
368 :
369 : /**
370 : * Computes the tensor product of a 1D quadrature rule and a 2D
371 : * quadrature rule. Used in the init_3D routines for prismatic
372 : * element types.
373 : */
374 : void tensor_product_prism (const QBase & q1D, const QBase & q2D);
375 :
376 : /**
377 : * The spatial dimension of the quadrature rule.
378 : */
379 : unsigned int _dim;
380 :
381 : /**
382 : * The polynomial order which the quadrature rule is capable of
383 : * integrating exactly.
384 : */
385 : Order _order;
386 :
387 : /**
388 : * The type of element for which the current values have been
389 : * computed.
390 : */
391 : ElemType _type;
392 :
393 : /**
394 : * The element for which the current values were computed, or
395 : * nullptr if values were computed without a specific element.
396 : */
397 : const Elem * _elem;
398 :
399 : /**
400 : * The p-level of the element for which the current values have
401 : * been computed.
402 : */
403 : unsigned int _p_level;
404 :
405 : /**
406 : * The locations of the quadrature points in reference element
407 : * space.
408 : */
409 : std::vector<Point> _points;
410 :
411 : /**
412 : * The quadrature weights. The order of the weights matches the
413 : * ordering of the _points vector.
414 : */
415 : std::vector<Real> _weights;
416 : };
417 :
418 : } // namespace libMesh
419 :
420 : #endif // LIBMESH_QUADRATURE_H
|
