Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2988-ge5f3bfba9d 2025-04-01 14:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
quadrature.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
18
19#include <algorithm>
20#include <array>
21#include <cmath>
22#include <limits>
23#include <memory>
24#include <vector>
25
27
28
29#ifndef DOXYGEN
30template <>
32 : is_tensor_product_flag(false)
33{}
34
35template <>
36Quadrature<0>::Quadrature(const unsigned int n_q)
38 , weights(n_q, 0)
39 , is_tensor_product_flag(false)
40{}
41#endif
42
43
44
45template <int dim>
47 : is_tensor_product_flag(dim == 1)
48{}
49
50
51
52template <int dim>
53Quadrature<dim>::Quadrature(const unsigned int n_q)
54 : quadrature_points(n_q, Point<dim>())
55 , weights(n_q, 0)
56 , is_tensor_product_flag(dim == 1)
57{}
58
59
60
61template <int dim>
62void
64 const ArrayView<const double> &weights)
65{
66 this->weights.clear();
67 if (weights.size() > 0)
68 {
69 AssertDimension(weights.size(), points.size());
70 this->weights.insert(this->weights.end(), weights.begin(), weights.end());
71 }
72 else
73 this->weights.resize(points.size(),
74 std::numeric_limits<double>::infinity());
75
76 quadrature_points.clear();
77 quadrature_points.insert(quadrature_points.end(),
78 points.begin(),
79 points.end());
80
81 is_tensor_product_flag = dim == 1;
82}
83
84
85
86template <int dim>
87Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points,
88 const std::vector<double> &weights)
89 : quadrature_points(points)
90 , weights(weights)
91 , is_tensor_product_flag(dim == 1)
92{
93 Assert(weights.size() == points.size(),
94 ExcDimensionMismatch(weights.size(), points.size()));
95}
96
97
98
99template <int dim>
101 std::vector<double> &&weights)
102 : quadrature_points(std::move(points))
103 , weights(std::move(weights))
104 , is_tensor_product_flag(dim == 1)
105{
106 Assert(weights.size() == points.size(),
107 ExcDimensionMismatch(weights.size(), points.size()));
108}
109
110
111
112template <int dim>
113Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points)
114 : quadrature_points(points)
115 , weights(points.size(), std::numeric_limits<double>::infinity())
116 , is_tensor_product_flag(dim == 1)
117{
118 Assert(weights.size() == points.size(),
119 ExcDimensionMismatch(weights.size(), points.size()));
120}
121
122
123
124template <int dim>
126 : quadrature_points(std::vector<Point<dim>>(1, point))
127 , weights(std::vector<double>(1, 1.))
128 , is_tensor_product_flag(true)
129 , tensor_basis(new std::array<Quadrature<1>, dim>())
130{
131 for (unsigned int i = 0; i < dim; ++i)
132 {
133 const std::vector<Point<1>> quad_vec_1d(1, Point<1>(point[i]));
134 (*tensor_basis)[i] = Quadrature<1>(quad_vec_1d, weights);
135 }
137
138
139
140#ifndef DOXYGEN
141template <>
143 : quadrature_points(std::vector<Point<1>>(1, point))
144 , weights(std::vector<double>(1, 1.))
145 , is_tensor_product_flag(true)
147
148
149
150template <>
152 : is_tensor_product_flag(false)
153{
154 Assert(false, ExcImpossibleInDim(0));
155}
156
157
158
159template <>
160Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
161{
162 Assert(false, ExcImpossibleInDim(0));
163}
164#endif // DOXYGEN
165
166
168template <int dim>
170 : quadrature_points(q1.size() * q2.size())
171 , weights(q1.size() * q2.size())
172 , is_tensor_product_flag(q1.is_tensor_product())
173{
174 unsigned int present_index = 0;
175 for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
176 for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
177 {
178 // compose coordinates of new quadrature point by tensor product in the
179 // last component
180 for (unsigned int d = 0; d < dim - 1; ++d)
181 quadrature_points[present_index][d] = q1.point(i1)[d];
182 quadrature_points[present_index][dim - 1] = q2.point(i2)[0];
183
184 weights[present_index] = q1.weight(i1) * q2.weight(i2);
186 ++present_index;
187 }
188
189 if constexpr (running_in_debug_mode())
190 {
191 if (size() > 0)
192 {
193 double sum = 0;
194 for (unsigned int i = 0; i < size(); ++i)
195 sum += weights[i];
196 // we cannot guarantee the sum of weights to be exactly one, but it
197 // should be near that.
198 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
199 }
200 }
201
202 if (is_tensor_product_flag)
203 {
204 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
205 for (unsigned int i = 0; i < dim - 1; ++i)
206 (*tensor_basis)[i] = q1.get_tensor_basis()[i];
207 (*tensor_basis)[dim - 1] = q2;
209}
210
211
212#ifndef DOXYGEN
213template <>
214Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2)
215 : quadrature_points(q2.size())
216 , weights(q2.size())
217 , is_tensor_product_flag(true)
218{
219 unsigned int present_index = 0;
220 for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
221 {
222 // compose coordinates of new quadrature point by tensor product in the
223 // last component
224 quadrature_points[present_index][0] = q2.point(i2)[0];
225
226 weights[present_index] = q2.weight(i2);
227
228 ++present_index;
229 }
230
231 if constexpr (running_in_debug_mode())
232 {
233 if (size() > 0)
234 {
235 double sum = 0;
236 for (unsigned int i = 0; i < size(); ++i)
237 sum += weights[i];
238 // we cannot guarantee the sum of weights to be exactly one, but it
239 // should be near that.
240 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
241 }
242 }
243}
245
246
247template <>
250 , quadrature_points(1)
251 , weights(1, 1.)
252 , is_tensor_product_flag(false)
253{}
254
255
256template <>
259{
260 // this function should never be called -- this should be the copy constructor
261 // in 1d...
262 Assert(false, ExcImpossibleInDim(1));
263}
264#endif // DOXYGEN
265
266
267
268template <int dim>
271 , quadrature_points(Utilities::fixed_power<dim>(q.size()))
272 , weights(Utilities::fixed_power<dim>(q.size()))
273 , is_tensor_product_flag(true)
274{
275 Assert(dim <= 3, ExcNotImplemented());
276
277 const unsigned int n0 = q.size();
278 const unsigned int n1 = (dim > 1) ? n0 : 1;
279 const unsigned int n2 = (dim > 2) ? n0 : 1;
280
281 unsigned int k = 0;
282 for (unsigned int i2 = 0; i2 < n2; ++i2)
283 for (unsigned int i1 = 0; i1 < n1; ++i1)
284 for (unsigned int i0 = 0; i0 < n0; ++i0)
285 {
286 quadrature_points[k][0] = q.point(i0)[0];
287 if (dim > 1)
288 quadrature_points[k][1] = q.point(i1)[0];
289 if (dim > 2)
290 quadrature_points[k][2] = q.point(i2)[0];
291 weights[k] = q.weight(i0);
292 if (dim > 1)
293 weights[k] *= q.weight(i1);
294 if (dim > 2)
295 weights[k] *= q.weight(i2);
296 ++k;
297 }
298
299 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
300 for (unsigned int i = 0; i < dim; ++i)
301 (*tensor_basis)[i] = q;
302}
303
304
305
306template <int dim>
309 , quadrature_points(q.quadrature_points)
310 , weights(q.weights)
311 , is_tensor_product_flag(q.is_tensor_product_flag)
312{
313 if (dim > 1 && is_tensor_product_flag)
315 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
316}
317
318
319
320template <int dim>
323{
324 weights = q.weights;
325 quadrature_points = q.quadrature_points;
326 is_tensor_product_flag = q.is_tensor_product_flag;
327 if (dim > 1 && is_tensor_product_flag)
328 {
329 if (tensor_basis == nullptr)
330 tensor_basis =
331 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
332 else
333 *tensor_basis = *q.tensor_basis;
334 }
335 return *this;
336}
337
338
339
340template <int dim>
341bool
343{
344 return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
345}
346
347
349template <int dim>
350std::size_t
356
357
358
359template <int dim>
360typename std::conditional_t<dim == 1,
361 std::array<Quadrature<1>, dim>,
362 const std::array<Quadrature<1>, dim> &>
364{
365 Assert(this->is_tensor_product_flag == true,
366 ExcMessage("This function only makes sense if "
367 "this object represents a tensor product!"));
368 Assert(tensor_basis != nullptr, ExcInternalError());
369
370 return *tensor_basis;
371}
372
373
374#ifndef DOXYGEN
375template <>
376std::array<Quadrature<1>, 1>
378{
379 Assert(this->is_tensor_product_flag == true,
380 ExcMessage("This function only makes sense if "
381 "this object represents a tensor product!"));
382
383 return std::array<Quadrature<1>, 1>{{*this}};
384}
385#endif
386
387
388
389//---------------------------------------------------------------------------
390template <int dim>
392 : Quadrature<dim>(qx.size())
393{
394 Assert(dim == 1, ExcImpossibleInDim(dim));
395 unsigned int k = 0;
396 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
397 {
398 this->quadrature_points[k][0] = qx.point(k1)[0];
399 this->weights[k++] = qx.weight(k1);
400 }
401 Assert(k == this->size(), ExcInternalError());
402 this->is_tensor_product_flag = true;
403}
404
405
406
407template <int dim>
409 const Quadrature<1> &qy)
410 : Quadrature<dim>(qx.size() * qy.size())
411{
412 Assert(dim == 2, ExcImpossibleInDim(dim));
413
414 // placate compiler in the dim == 1 case
415 constexpr int dim_1 = dim == 2 ? 1 : 0;
416
417 unsigned int k = 0;
418 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
419 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
420 {
421 this->quadrature_points[k][0] = qx.point(k1)[0];
422 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
423 this->weights[k++] = qx.weight(k1) * qy.weight(k2);
424 }
425 Assert(k == this->size(), ExcInternalError());
426
427 this->is_tensor_product_flag = true;
428 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
429 (*this->tensor_basis)[0] = qx;
430 (*this->tensor_basis)[dim_1] = qy;
431}
432
433
434
435template <int dim>
437 const Quadrature<1> &qy,
438 const Quadrature<1> &qz)
439 : Quadrature<dim>(qx.size() * qy.size() * qz.size())
440{
441 Assert(dim == 3, ExcImpossibleInDim(dim));
442
443 // placate compiler in lower dimensions
444 constexpr int dim_1 = dim == 3 ? 1 : 0;
445 constexpr int dim_2 = dim == 3 ? 2 : 0;
446
447 unsigned int k = 0;
448 for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
449 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
450 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
451 {
452 this->quadrature_points[k][0] = qx.point(k1)[0];
453 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
454 this->quadrature_points[k][dim_2] = qz.point(k3)[0];
455 this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
456 }
457 Assert(k == this->size(), ExcInternalError());
458
459 this->is_tensor_product_flag = true;
460 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
461 (*this->tensor_basis)[0] = qx;
462 (*this->tensor_basis)[dim_1] = qy;
463 (*this->tensor_basis)[dim_2] = qz;
464}
465
466
467
468// ------------------------------------------------------------ //
469
470namespace internal
471{
472 namespace QIteratedImplementation
473 {
474 namespace
475 {
476 bool
477 uses_both_endpoints(const Quadrature<1> &base_quadrature)
478 {
479 const bool at_left =
480 std::any_of(base_quadrature.get_points().cbegin(),
481 base_quadrature.get_points().cend(),
482 [](const Point<1> &p) { return p == Point<1>{0.}; });
483 const bool at_right =
484 std::any_of(base_quadrature.get_points().cbegin(),
485 base_quadrature.get_points().cend(),
486 [](const Point<1> &p) { return p == Point<1>{1.}; });
487 return (at_left && at_right);
488 }
489
490 std::vector<Point<1>>
491 create_equidistant_interval_points(const unsigned int n_copies)
492 {
493 std::vector<Point<1>> support_points(n_copies + 1);
494
495 for (unsigned int copy = 0; copy < n_copies; ++copy)
496 support_points[copy][0] =
497 static_cast<double>(copy) / static_cast<double>(n_copies);
498
499 support_points[n_copies][0] = 1.0;
500
501 return support_points;
502 }
503 } // namespace
504 } // namespace QIteratedImplementation
505} // namespace internal
506
507
508
509template <>
510QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
511 : Quadrature<0>()
512{
514}
515
516
517
518template <>
519QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
520 : Quadrature<0>()
521{
523}
524
525
526
527template <>
529 const std::vector<Point<1>> &intervals)
530 : Quadrature<1>(
531 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
532 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
533 base_quadrature.size() * (intervals.size() - 1))
534{
535 Assert(base_quadrature.size() > 0, ExcNotInitialized());
536 Assert(intervals.size() > 1, ExcZero());
537
538 const unsigned int n_copies = intervals.size() - 1;
539
540 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
541 // we don't have to skip some points in order to get a reasonable quadrature
542 // formula
543 {
544 unsigned int next_point = 0;
545 for (unsigned int copy = 0; copy < n_copies; ++copy)
546 for (unsigned int q_point = 0; q_point < base_quadrature.size();
547 ++q_point)
548 {
549 this->quadrature_points[next_point] =
550 Point<1>(base_quadrature.point(q_point)[0] *
551 (intervals[copy + 1][0] - intervals[copy][0]) +
552 intervals[copy][0]);
553 this->weights[next_point] =
554 base_quadrature.weight(q_point) *
555 (intervals[copy + 1][0] - intervals[copy][0]);
556
557 ++next_point;
558 }
559 }
560 else
561 // skip doubly available points
562 {
563 const unsigned int left_index =
564 std::distance(base_quadrature.get_points().begin(),
565 std::find_if(base_quadrature.get_points().cbegin(),
566 base_quadrature.get_points().cend(),
567 [](const Point<1> &p) {
568 return p == Point<1>{0.};
569 }));
570
571 const unsigned int right_index =
572 std::distance(base_quadrature.get_points().begin(),
573 std::find_if(base_quadrature.get_points().cbegin(),
574 base_quadrature.get_points().cend(),
575 [](const Point<1> &p) {
576 return p == Point<1>{1.};
577 }));
578
579 const unsigned double_point_offset =
580 left_index + (base_quadrature.size() - right_index);
581
582 for (unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
583 for (unsigned int q_point = 0; q_point < base_quadrature.size();
584 ++q_point)
585 {
586 // skip the left point of this copy since we have already entered it
587 // the last time
588 if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
589 {
590 Assert(this->quadrature_points[next_point - double_point_offset]
591 .distance(Point<1>(
592 base_quadrature.point(q_point)[0] *
593 (intervals[copy + 1][0] - intervals[copy][0]) +
594 intervals[copy][0])) < 1e-10 /*tolerance*/,
596
597 this->weights[next_point - double_point_offset] +=
598 base_quadrature.weight(q_point) *
599 (intervals[copy + 1][0] - intervals[copy][0]);
600
601 continue;
602 }
603
604 this->quadrature_points[next_point] =
605 Point<1>(base_quadrature.point(q_point)[0] *
606 (intervals[copy + 1][0] - intervals[copy][0]) +
607 intervals[copy][0]);
608
609 // if this is the rightmost point of one of the non-last copies:
610 // give it the double weight
611 this->weights[next_point] =
612 base_quadrature.weight(q_point) *
613 (intervals[copy + 1][0] - intervals[copy][0]);
614
615 ++next_point;
616 }
617 }
618
619 // make sure that there is no rounding error for 0.0 and 1.0, since there
620 // are multiple asserts in the library checking for equality without
621 // tolerances
622 for (auto &i : this->quadrature_points)
623 if (std::abs(i[0] - 0.0) < 1e-12)
624 i[0] = 0.0;
625 else if (std::abs(i[0] - 1.0) < 1e-12)
626 i[0] = 1.0;
627
628 if constexpr (running_in_debug_mode())
629 {
630 double sum_of_weights = 0;
631 for (unsigned int i = 0; i < this->size(); ++i)
632 sum_of_weights += this->weight(i);
633 Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
634 }
635}
636
637
638
639template <>
641 const unsigned int n_copies)
642 : QIterated<1>(
643 base_quadrature,
644 internal::QIteratedImplementation::create_equidistant_interval_points(
645 n_copies))
646{
647 Assert(base_quadrature.size() > 0, ExcNotInitialized());
648 Assert(n_copies > 0, ExcZero());
649}
650
651
652
653// construct higher dimensional quadrature formula by tensor product
654// of lower dimensional iterated quadrature formulae
655template <int dim>
657 const std::vector<Point<1>> &intervals)
658 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, intervals),
659 QIterated<1>(base_quadrature, intervals))
660{}
661
662
663
664template <int dim>
666 const unsigned int n_copies)
667 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, n_copies),
668 QIterated<1>(base_quadrature, n_copies))
669{}
670
671
672
673// explicit instantiations; note: we need them all for all dimensions
674template class Quadrature<0>;
675template class Quadrature<1>;
676template class Quadrature<2>;
677template class Quadrature<3>;
678template class QAnisotropic<1>;
679template class QAnisotropic<2>;
680template class QAnisotropic<3>;
681template class QIterated<1>;
682template class QIterated<2>;
683template class QIterated<3>;
684
iterator begin() const
Definition array_view.h:707
iterator end() const
Definition array_view.h:716
std::size_t size() const
Definition array_view.h:689
Definition point.h:113
QAnisotropic(const Quadrature< 1 > &qx)
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
std::vector< Point< dim > > quadrature_points
Definition quadrature.h:354
void initialize(const ArrayView< const Point< dim > > &points, const ArrayView< const double > &weights={})
Definition quadrature.cc:63
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition quadrature.h:375
Quadrature & operator=(const Quadrature< dim > &)
std::size_t memory_consumption() const
const Point< dim > & point(const unsigned int i) const
bool is_tensor_product_flag
Definition quadrature.h:369
double weight(const unsigned int i) const
bool operator==(const Quadrature< dim > &p) const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
std::vector< double > weights
Definition quadrature.h:360
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::size_t size
Definition mpi.cc:745
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void copy(const T *begin, const T *end, U *dest)
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)