Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 1784128 : QBase::QBase(unsigned int d,
28 1784128 : Order o) :
29 731569 : allow_rules_with_negative_weights(true),
30 731569 : allow_nodal_pyramid_quadrature(false),
31 731569 : _dim(d),
32 731569 : _order(o),
33 731569 : _type(INVALID_ELEM),
34 731569 : _elem(nullptr),
35 1784128 : _p_level(0)
36 1784128 : {}
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 42726542 : void QBase::init(const Elem & elem,
66 : unsigned int p)
67 : {
68 12406182 : libmesh_assert_equal_to(elem.dim(), _dim);
69 :
70 : // Default to the element p_level() value
71 42726542 : if (p == invalid_uint)
72 19506214 : p = elem.p_level();
73 :
74 42726542 : 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 42726542 : if (t == _type && p == _p_level && !elem.runtime_topology())
84 12315640 : return;
85 : else
86 : {
87 314454 : _elem = &elem;
88 314454 : _type = t;
89 314454 : _p_level = p;
90 : }
91 :
92 314454 : switch(_elem->dim())
93 : {
94 2531 : case 0:
95 2531 : this->init_0D();
96 :
97 2531 : return;
98 :
99 59901 : case 1:
100 59901 : this->init_1D();
101 :
102 59901 : return;
103 :
104 161405 : case 2:
105 161405 : this->init_2D();
106 :
107 161405 : return;
108 :
109 90617 : case 3:
110 90617 : this->init_3D();
111 :
112 90617 : return;
113 :
114 0 : default:
115 0 : libmesh_error_msg("Invalid dimension _dim = " << _dim);
116 : }
117 : }
118 :
119 :
120 :
121 902088 : 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 902088 : 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 271641 : 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 902088 : if (t == _type && p == _p_level)
141 2086 : return;
142 : else
143 : {
144 894803 : _elem = nullptr;
145 894803 : _type = t;
146 894803 : _p_level = p;
147 : }
148 :
149 894803 : switch(_dim)
150 : {
151 5 : case 0:
152 5 : this->init_0D();
153 :
154 5 : return;
155 :
156 884896 : case 1:
157 884896 : this->init_1D();
158 :
159 884896 : return;
160 :
161 9560 : case 2:
162 9560 : this->init_2D();
163 :
164 9560 : return;
165 :
166 342 : case 3:
167 342 : this->init_3D();
168 :
169 342 : return;
170 :
171 0 : default:
172 0 : libmesh_error_msg("Invalid dimension _dim = " << _dim);
173 : }
174 : }
175 :
176 :
177 :
178 27287 : void QBase::init(const QBase & other_rule)
179 : {
180 27287 : if (other_rule._elem)
181 20689 : this->init(*other_rule._elem, other_rule._p_level);
182 : else
183 6598 : this->init(other_rule._type, other_rule._p_level, true);
184 27287 : }
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 2536 : void QBase::init_0D()
199 : {
200 2536 : _points.resize(1);
201 2536 : _weights.resize(1);
202 2536 : _points[0] = Point(0.);
203 2536 : _weights[0] = 1.0;
204 2536 : }
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 133 : 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 133 : h_new = new_range.second - new_range.first,
230 133 : 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 133 : Real scfact = h_new/h_old;
241 :
242 : // We're mapping from old_range -> new_range
243 483 : for (auto i : index_range(_points))
244 : {
245 350 : _points[i](0) = new_range.first +
246 350 : (_points[i](0) - old_range.first) * scfact;
247 :
248 : // Scale the weights
249 450 : _weights[i] *= scfact;
250 : }
251 133 : }
252 :
253 :
254 :
255 :
256 108460 : void QBase::tensor_product_quad(const QBase & q1D)
257 : {
258 :
259 33028 : const unsigned int np = q1D.n_points();
260 :
261 108460 : _points.resize(np * np);
262 :
263 108460 : _weights.resize(np * np);
264 :
265 33028 : unsigned int q=0;
266 :
267 502147 : for (unsigned int j=0; j<np; j++)
268 2647644 : for (unsigned int i=0; i<np; i++)
269 : {
270 2977798 : _points[q](0) = q1D.qp(i)(0);
271 2253957 : _points[q](1) = q1D.qp(j)(0);
272 :
273 2253957 : _weights[q] = q1D.w(i)*q1D.w(j);
274 :
275 2253957 : q++;
276 : }
277 108460 : }
278 :
279 :
280 :
281 :
282 :
283 8557 : void QBase::tensor_product_hex(const QBase & q1D)
284 : {
285 2243 : const unsigned int np = q1D.n_points();
286 :
287 8557 : _points.resize(np * np * np);
288 :
289 8557 : _weights.resize(np * np * np);
290 :
291 2243 : unsigned int q=0;
292 :
293 35772 : for (unsigned int k=0; k<np; k++)
294 143408 : for (unsigned int j=0; j<np; j++)
295 995522 : for (unsigned int i=0; i<np; i++)
296 : {
297 1123382 : _points[q](0) = q1D.qp(i)(0);
298 879329 : _points[q](1) = q1D.qp(j)(0);
299 879329 : _points[q](2) = q1D.qp(k)(0);
300 :
301 1123382 : _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
302 :
303 879329 : q++;
304 : }
305 8557 : }
306 :
307 :
308 :
309 :
310 8653 : void QBase::tensor_product_prism(const QBase & q1D, const QBase & q2D)
311 : {
312 2403 : const unsigned int n_points1D = q1D.n_points();
313 2403 : const unsigned int n_points2D = q2D.n_points();
314 :
315 8653 : _points.resize (n_points1D * n_points2D);
316 8653 : _weights.resize (n_points1D * n_points2D);
317 :
318 2403 : unsigned int q=0;
319 :
320 38928 : for (unsigned int j=0; j<n_points1D; j++)
321 424089 : for (unsigned int i=0; i<n_points2D; i++)
322 : {
323 503799 : _points[q](0) = q2D.qp(i)(0);
324 393814 : _points[q](1) = q2D.qp(i)(1);
325 393814 : _points[q](2) = q1D.qp(j)(0);
326 :
327 393814 : _weights[q] = q2D.w(i) * q1D.w(j);
328 :
329 393814 : q++;
330 : }
331 :
332 8653 : }
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
|
