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 : // Local includes
20 : #include "libmesh/elem.h"
21 : #include "libmesh/quadrature.h"
22 : #include "libmesh/int_range.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 14184741 : QBase::QBase(unsigned int d,
28 14184741 : Order o) :
29 12877758 : allow_rules_with_negative_weights(true),
30 12877758 : allow_nodal_pyramid_quadrature(false),
31 12877758 : _dim(d),
32 12877758 : _order(o),
33 12877758 : _type(INVALID_ELEM),
34 12877758 : _elem(nullptr),
35 14184741 : _p_level(0)
36 14184741 : {}
37 :
38 0 : std::unique_ptr<QBase> QBase::clone() const
39 : {
40 0 : return QBase::build(this->type(), this->get_dim(), this->get_order());
41 : }
42 :
43 0 : void QBase::print_info(std::ostream & os) const
44 : {
45 0 : libmesh_assert(!_points.empty());
46 0 : libmesh_assert(!_weights.empty());
47 :
48 0 : Real summed_weights=0;
49 0 : os << "N_Q_Points=" << this->n_points() << std::endl << std::endl;
50 0 : for (auto qpoint: index_range(_points))
51 : {
52 0 : os << " Point " << qpoint << ":\n"
53 0 : << " "
54 0 : << _points[qpoint]
55 : << "\n Weight:\n "
56 0 : << " w=" << _weights[qpoint] << "\n" << std::endl;
57 :
58 0 : summed_weights += _weights[qpoint];
59 : }
60 0 : os << "Summed Weights: " << summed_weights << std::endl;
61 0 : }
62 :
63 :
64 :
65 148685485 : void QBase::init(const Elem & elem,
66 : unsigned int p)
67 : {
68 12434902 : libmesh_assert_equal_to(elem.dim(), _dim);
69 :
70 : // Default to the element p_level() value
71 148685485 : if (p == invalid_uint)
72 19788530 : p = elem.p_level();
73 :
74 148685485 : ElemType t = elem.type();
75 :
76 : // check to see if we have already
77 : // done the work for this quadrature rule
78 : //
79 : // If we have something like a Polygon subclass then we're going to
80 : // need to recompute to be safe; even if we're using the same
81 : // element, it might have been distorted enough that its subtriangle
82 : // triangulation has been changed.
83 148685485 : if (t == _type && p == _p_level && !elem.runtime_topology())
84 12341737 : return;
85 : else
86 : {
87 1437381 : _elem = &elem;
88 1437381 : _type = t;
89 1437381 : _p_level = p;
90 : }
91 :
92 1437381 : switch(_elem->dim())
93 : {
94 11906 : case 0:
95 11906 : this->init_0D();
96 :
97 11906 : return;
98 :
99 149278 : case 1:
100 149278 : this->init_1D();
101 :
102 149278 : return;
103 :
104 690353 : case 2:
105 690353 : this->init_2D();
106 :
107 690353 : return;
108 :
109 585844 : case 3:
110 585844 : this->init_3D();
111 :
112 585844 : return;
113 :
114 0 : default:
115 0 : libmesh_error_msg("Invalid dimension _dim = " << _dim);
116 : }
117 : }
118 :
119 :
120 :
121 6588450 : void QBase::init(const ElemType t,
122 : unsigned int p,
123 : bool simple_type_only)
124 : {
125 : // Some element types require data from a specific element, so can
126 : // only be used with newer APIs.
127 6588450 : if (t == C0POLYGON || t == C0POLYHEDRON)
128 0 : libmesh_error_msg("Code (see stack trace) used an outdated quadrature function overload.\n"
129 : "Quadrature rules on a C0Polygon are not defined by its ElemType alone.");
130 :
131 : // This API is dangerous to use on general meshes, which may include
132 : // element types where the desired quadrature depends on the
133 : // physical element, but we still want to be able to initialize
134 : // based on only a type for certain simple cases
135 323259 : if (t != EDGE2 && !simple_type_only)
136 : libmesh_deprecated();
137 :
138 : // check to see if we have already
139 : // done the work for this quadrature rule
140 6588450 : if (t == _type && p == _p_level)
141 2306 : return;
142 : else
143 : {
144 6515705 : _elem = nullptr;
145 6515705 : _type = t;
146 6515705 : _p_level = p;
147 : }
148 :
149 6515705 : switch(_dim)
150 : {
151 5 : case 0:
152 5 : this->init_0D();
153 :
154 5 : return;
155 :
156 6295767 : case 1:
157 6295767 : this->init_1D();
158 :
159 6295767 : return;
160 :
161 219346 : case 2:
162 219346 : this->init_2D();
163 :
164 219346 : return;
165 :
166 587 : case 3:
167 587 : this->init_3D();
168 :
169 587 : return;
170 :
171 0 : default:
172 0 : libmesh_error_msg("Invalid dimension _dim = " << _dim);
173 : }
174 : }
175 :
176 :
177 :
178 170192 : void QBase::init(const QBase & other_rule)
179 : {
180 170192 : if (other_rule._elem)
181 103638 : this->init(*other_rule._elem, other_rule._p_level);
182 : else
183 66554 : this->init(other_rule._type, other_rule._p_level, true);
184 170192 : }
185 :
186 :
187 :
188 0 : void QBase::init (const Elem & elem,
189 : const std::vector<Real> & /* vertex_distance_func */,
190 : unsigned int p_level)
191 : {
192 : // dispatch generic implementation
193 0 : this->init(elem.type(), p_level);
194 0 : }
195 :
196 :
197 :
198 11911 : void QBase::init_0D()
199 : {
200 11911 : _points.resize(1);
201 11911 : _weights.resize(1);
202 11911 : _points[0] = Point(0.);
203 11911 : _weights[0] = 1.0;
204 11911 : }
205 :
206 :
207 :
208 0 : void QBase::init_2D ()
209 : {
210 0 : libmesh_not_implemented();
211 : }
212 :
213 :
214 :
215 0 : void QBase::init_3D ()
216 : {
217 0 : libmesh_not_implemented();
218 : }
219 :
220 :
221 :
222 1349 : void QBase::scale(std::pair<Real, Real> old_range,
223 : std::pair<Real, Real> new_range)
224 : {
225 : // Make sure we are in 1D
226 38 : libmesh_assert_equal_to (_dim, 1);
227 :
228 : Real
229 1349 : h_new = new_range.second - new_range.first,
230 1349 : h_old = old_range.second - old_range.first;
231 :
232 : // Make sure that we have sane ranges
233 38 : libmesh_assert_greater (h_new, 0.);
234 38 : libmesh_assert_greater (h_old, 0.);
235 :
236 : // Make sure there are some points
237 38 : libmesh_assert_greater (_points.size(), 0);
238 :
239 : // Compute the scale factor
240 1349 : Real scfact = h_new/h_old;
241 :
242 : // We're mapping from old_range -> new_range
243 4899 : for (auto i : index_range(_points))
244 : {
245 3550 : _points[i](0) = new_range.first +
246 3550 : (_points[i](0) - old_range.first) * scfact;
247 :
248 : // Scale the weights
249 3650 : _weights[i] *= scfact;
250 : }
251 1349 : }
252 :
253 :
254 :
255 :
256 507140 : void QBase::tensor_product_quad(const QBase & q1D)
257 : {
258 :
259 35809 : const unsigned int np = q1D.n_points();
260 :
261 507140 : _points.resize(np * np);
262 :
263 507140 : _weights.resize(np * np);
264 :
265 35809 : unsigned int q=0;
266 :
267 2451236 : for (unsigned int j=0; j<np; j++)
268 14654952 : for (unsigned int i=0; i<np; i++)
269 : {
270 13855603 : _points[q](0) = q1D.qp(i)(0);
271 12710856 : _points[q](1) = q1D.qp(j)(0);
272 :
273 12710856 : _weights[q] = q1D.w(i)*q1D.w(j);
274 :
275 12710856 : q++;
276 : }
277 507140 : }
278 :
279 :
280 :
281 :
282 :
283 65165 : void QBase::tensor_product_hex(const QBase & q1D)
284 : {
285 3044 : const unsigned int np = q1D.n_points();
286 :
287 65165 : _points.resize(np * np * np);
288 :
289 65165 : _weights.resize(np * np * np);
290 :
291 3044 : unsigned int q=0;
292 :
293 256129 : for (unsigned int k=0; k<np; k++)
294 1044030 : for (unsigned int j=0; j<np; j++)
295 8502276 : for (unsigned int i=0; i<np; i++)
296 : {
297 7903231 : _points[q](0) = q1D.qp(i)(0);
298 7649210 : _points[q](1) = q1D.qp(j)(0);
299 7649210 : _points[q](2) = q1D.qp(k)(0);
300 :
301 7903231 : _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
302 :
303 7649210 : q++;
304 : }
305 65165 : }
306 :
307 :
308 :
309 :
310 211710 : void QBase::tensor_product_prism(const QBase & q1D, const QBase & q2D)
311 : {
312 10741 : const unsigned int n_points1D = q1D.n_points();
313 10741 : const unsigned int n_points2D = q2D.n_points();
314 :
315 211710 : _points.resize (n_points1D * n_points2D);
316 211710 : _weights.resize (n_points1D * n_points2D);
317 :
318 10741 : unsigned int q=0;
319 :
320 823453 : for (unsigned int j=0; j<n_points1D; j++)
321 5912377 : for (unsigned int i=0; i<n_points2D; i++)
322 : {
323 5619695 : _points[q](0) = q2D.qp(i)(0);
324 5300634 : _points[q](1) = q2D.qp(i)(1);
325 5300634 : _points[q](2) = q1D.qp(j)(0);
326 :
327 5300634 : _weights[q] = q2D.w(i) * q1D.w(j);
328 :
329 5300634 : q++;
330 : }
331 :
332 211710 : }
333 :
334 :
335 :
336 :
337 0 : std::ostream & operator << (std::ostream & os, const QBase & q)
338 : {
339 0 : q.print_info(os);
340 0 : return os;
341 : }
342 :
343 : } // namespace libMesh
|
