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 11823642 : QBase::QBase(unsigned int d,
28 11823642 : Order o) :
29 10695853 : allow_rules_with_negative_weights(true),
30 10695853 : allow_nodal_pyramid_quadrature(false),
31 10695853 : _dim(d),
32 10695853 : _order(o),
33 10695853 : _type(INVALID_ELEM),
34 10695853 : _elem(nullptr),
35 11823642 : _p_level(0)
36 11823642 : {}
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 145893607 : void QBase::init(const Elem & elem,
66 : unsigned int p)
67 : {
68 12209300 : libmesh_assert_equal_to(elem.dim(), _dim);
69 :
70 : // Default to the element p_level() value
71 145893607 : if (p == invalid_uint)
72 19334826 : p = elem.p_level();
73 :
74 145893607 : 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 145893607 : if (t == _type && p == _p_level && !elem.runtime_topology())
84 12129410 : return;
85 : else
86 : {
87 1149438 : _elem = &elem;
88 1149438 : _type = t;
89 1149438 : _p_level = p;
90 : }
91 :
92 1149438 : switch(_elem->dim())
93 : {
94 11909 : case 0:
95 11909 : this->init_0D();
96 :
97 11909 : return;
98 :
99 149339 : case 1:
100 149339 : this->init_1D();
101 :
102 149339 : return;
103 :
104 661082 : case 2:
105 661082 : this->init_2D();
106 :
107 661082 : return;
108 :
109 327108 : case 3:
110 327108 : this->init_3D();
111 :
112 327108 : return;
113 :
114 0 : default:
115 0 : libmesh_error_msg("Invalid dimension _dim = " << _dim);
116 : }
117 : }
118 :
119 :
120 :
121 5481711 : 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 5481711 : 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 278910 : 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 5481711 : if (t == _type && p == _p_level)
141 2306 : return;
142 : else
143 : {
144 5408966 : _elem = nullptr;
145 5408966 : _type = t;
146 5408966 : _p_level = p;
147 : }
148 :
149 5408966 : switch(_dim)
150 : {
151 5 : case 0:
152 5 : this->init_0D();
153 :
154 5 : return;
155 :
156 5371962 : case 1:
157 5371962 : this->init_1D();
158 :
159 5371962 : return;
160 :
161 36412 : case 2:
162 36412 : this->init_2D();
163 :
164 36412 : 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 144749 : void QBase::init(const QBase & other_rule)
179 : {
180 144749 : if (other_rule._elem)
181 78195 : this->init(*other_rule._elem, other_rule._p_level);
182 : else
183 66554 : this->init(other_rule._type, other_rule._p_level, true);
184 144749 : }
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 11914 : void QBase::init_0D()
199 : {
200 11914 : _points.resize(1);
201 11914 : _weights.resize(1);
202 11914 : _points[0] = Point(0.);
203 11914 : _weights[0] = 1.0;
204 11914 : }
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 492480 : void QBase::tensor_product_quad(const QBase & q1D)
257 : {
258 :
259 35045 : const unsigned int np = q1D.n_points();
260 :
261 492480 : _points.resize(np * np);
262 :
263 492480 : _weights.resize(np * np);
264 :
265 35045 : unsigned int q=0;
266 :
267 2404651 : for (unsigned int j=0; j<np; j++)
268 14551130 : for (unsigned int i=0; i<np; i++)
269 : {
270 13779954 : _points[q](0) = q1D.qp(i)(0);
271 12638959 : _points[q](1) = q1D.qp(j)(0);
272 :
273 12638959 : _weights[q] = q1D.w(i)*q1D.w(j);
274 :
275 12638959 : q++;
276 : }
277 492480 : }
278 :
279 :
280 :
281 :
282 :
283 40249 : void QBase::tensor_product_hex(const QBase & q1D)
284 : {
285 2085 : const unsigned int np = q1D.n_points();
286 :
287 40249 : _points.resize(np * np * np);
288 :
289 40249 : _weights.resize(np * np * np);
290 :
291 2085 : unsigned int q=0;
292 :
293 175761 : for (unsigned int k=0; k<np; k++)
294 860814 : for (unsigned int j=0; j<np; j++)
295 8068404 : for (unsigned int i=0; i<np; i++)
296 : {
297 7584007 : _points[q](0) = q1D.qp(i)(0);
298 7343102 : _points[q](1) = q1D.qp(j)(0);
299 7343102 : _points[q](2) = q1D.qp(k)(0);
300 :
301 7584007 : _weights[q] = q1D.w(i) * q1D.w(j) * q1D.w(k);
302 :
303 7343102 : q++;
304 : }
305 40249 : }
306 :
307 :
308 :
309 :
310 28776 : void QBase::tensor_product_prism(const QBase & q1D, const QBase & q2D)
311 : {
312 2097 : const unsigned int n_points1D = q1D.n_points();
313 2097 : const unsigned int n_points2D = q2D.n_points();
314 :
315 28776 : _points.resize (n_points1D * n_points2D);
316 28776 : _weights.resize (n_points1D * n_points2D);
317 :
318 2097 : unsigned int q=0;
319 :
320 131136 : for (unsigned int j=0; j<n_points1D; j++)
321 1669451 : for (unsigned int i=0; i<n_points2D; i++)
322 : {
323 1668572 : _points[q](0) = q2D.qp(i)(0);
324 1567091 : _points[q](1) = q2D.qp(i)(1);
325 1567091 : _points[q](2) = q1D.qp(j)(0);
326 :
327 1567091 : _weights[q] = q2D.w(i) * q1D.w(j);
328 :
329 1567091 : q++;
330 : }
331 :
332 28776 : }
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
|
