Reference documentation for deal.II version GIT relicensing-1308-g50a4ffc717 2024-07-25 10:40:01+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\}}\)
Loading...
Searching...
No Matches
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);
136}
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)
146{}
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
167
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);
185
186 ++present_index;
187 }
188
189#ifdef DEBUG
190 if (size() > 0)
191 {
192 double sum = 0;
193 for (unsigned int i = 0; i < size(); ++i)
194 sum += weights[i];
195 // we cannot guarantee the sum of weights to be exactly one, but it should
196 // be near that.
197 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
198 }
199#endif
200
201 if (is_tensor_product_flag)
202 {
203 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
204 for (unsigned int i = 0; i < dim - 1; ++i)
205 (*tensor_basis)[i] = q1.get_tensor_basis()[i];
206 (*tensor_basis)[dim - 1] = q2;
208}
209
210
211#ifndef DOXYGEN
212template <>
213Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2)
214 : quadrature_points(q2.size())
215 , weights(q2.size())
216 , is_tensor_product_flag(true)
217{
218 unsigned int present_index = 0;
219 for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
220 {
221 // compose coordinates of new quadrature point by tensor product in the
222 // last component
223 quadrature_points[present_index][0] = q2.point(i2)[0];
224
225 weights[present_index] = q2.weight(i2);
226
227 ++present_index;
228 }
229
230# ifdef DEBUG
231 if (size() > 0)
233 double sum = 0;
234 for (unsigned int i = 0; i < size(); ++i)
235 sum += weights[i];
236 // we cannot guarantee the sum of weights to be exactly one, but it should
237 // be near that.
238 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
239 }
240# endif
241}
242
244
245template <>
247 : Subscriptor()
248 , quadrature_points(1)
249 , weights(1, 1.)
250 , is_tensor_product_flag(false)
251{}
252
253
254template <>
256 : Subscriptor()
257{
258 // this function should never be called -- this should be the copy constructor
259 // in 1d...
260 Assert(false, ExcImpossibleInDim(1));
261}
262#endif // DOXYGEN
263
264
265
266template <int dim>
268 : Subscriptor()
269 , quadrature_points(Utilities::fixed_power<dim>(q.size()))
270 , weights(Utilities::fixed_power<dim>(q.size()))
271 , is_tensor_product_flag(true)
272{
273 Assert(dim <= 3, ExcNotImplemented());
274
275 const unsigned int n0 = q.size();
276 const unsigned int n1 = (dim > 1) ? n0 : 1;
277 const unsigned int n2 = (dim > 2) ? n0 : 1;
278
279 unsigned int k = 0;
280 for (unsigned int i2 = 0; i2 < n2; ++i2)
281 for (unsigned int i1 = 0; i1 < n1; ++i1)
282 for (unsigned int i0 = 0; i0 < n0; ++i0)
283 {
284 quadrature_points[k][0] = q.point(i0)[0];
285 if (dim > 1)
286 quadrature_points[k][1] = q.point(i1)[0];
287 if (dim > 2)
288 quadrature_points[k][2] = q.point(i2)[0];
289 weights[k] = q.weight(i0);
290 if (dim > 1)
291 weights[k] *= q.weight(i1);
292 if (dim > 2)
293 weights[k] *= q.weight(i2);
294 ++k;
295 }
296
297 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
298 for (unsigned int i = 0; i < dim; ++i)
299 (*tensor_basis)[i] = q;
300}
301
302
303
304template <int dim>
306 : Subscriptor()
307 , quadrature_points(q.quadrature_points)
308 , weights(q.weights)
309 , is_tensor_product_flag(q.is_tensor_product_flag)
310{
311 if (dim > 1 && is_tensor_product_flag)
313 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
314}
315
316
317
318template <int dim>
321{
322 weights = q.weights;
323 quadrature_points = q.quadrature_points;
324 is_tensor_product_flag = q.is_tensor_product_flag;
325 if (dim > 1 && is_tensor_product_flag)
326 {
327 if (tensor_basis == nullptr)
328 tensor_basis =
329 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
330 else
331 *tensor_basis = *q.tensor_basis;
332 }
333 return *this;
334}
335
336
337
338template <int dim>
339bool
341{
342 return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
343}
344
345
346
347template <int dim>
348std::size_t
354
355
356
357template <int dim>
358typename std::conditional_t<dim == 1,
359 std::array<Quadrature<1>, dim>,
360 const std::array<Quadrature<1>, dim> &>
362{
363 Assert(this->is_tensor_product_flag == true,
364 ExcMessage("This function only makes sense if "
365 "this object represents a tensor product!"));
366 Assert(tensor_basis != nullptr, ExcInternalError());
367
368 return *tensor_basis;
369}
370
371
372#ifndef DOXYGEN
373template <>
374std::array<Quadrature<1>, 1>
376{
377 Assert(this->is_tensor_product_flag == true,
378 ExcMessage("This function only makes sense if "
379 "this object represents a tensor product!"));
380
381 return std::array<Quadrature<1>, 1>{{*this}};
382}
383#endif
384
385
386
387//---------------------------------------------------------------------------
388template <int dim>
390 : Quadrature<dim>(qx.size())
391{
392 Assert(dim == 1, ExcImpossibleInDim(dim));
393 unsigned int k = 0;
394 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
395 {
396 this->quadrature_points[k][0] = qx.point(k1)[0];
397 this->weights[k++] = qx.weight(k1);
398 }
399 Assert(k == this->size(), ExcInternalError());
400 this->is_tensor_product_flag = true;
401}
402
403
404
405template <int dim>
407 const Quadrature<1> &qy)
408 : Quadrature<dim>(qx.size() * qy.size())
409{
410 Assert(dim == 2, ExcImpossibleInDim(dim));
411
412 // placate compiler in the dim == 1 case
413 constexpr int dim_1 = dim == 2 ? 1 : 0;
414
415 unsigned int k = 0;
416 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
417 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
418 {
419 this->quadrature_points[k][0] = qx.point(k1)[0];
420 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
421 this->weights[k++] = qx.weight(k1) * qy.weight(k2);
422 }
423 Assert(k == this->size(), ExcInternalError());
424
425 this->is_tensor_product_flag = true;
426 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
427 (*this->tensor_basis)[0] = qx;
428 (*this->tensor_basis)[dim_1] = qy;
429}
430
431
432
433template <int dim>
435 const Quadrature<1> &qy,
436 const Quadrature<1> &qz)
437 : Quadrature<dim>(qx.size() * qy.size() * qz.size())
438{
439 Assert(dim == 3, ExcImpossibleInDim(dim));
440
441 // placate compiler in lower dimensions
442 constexpr int dim_1 = dim == 3 ? 1 : 0;
443 constexpr int dim_2 = dim == 3 ? 2 : 0;
444
445 unsigned int k = 0;
446 for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
447 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
448 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
449 {
450 this->quadrature_points[k][0] = qx.point(k1)[0];
451 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
452 this->quadrature_points[k][dim_2] = qz.point(k3)[0];
453 this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
454 }
455 Assert(k == this->size(), ExcInternalError());
456
457 this->is_tensor_product_flag = true;
458 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
459 (*this->tensor_basis)[0] = qx;
460 (*this->tensor_basis)[dim_1] = qy;
461 (*this->tensor_basis)[dim_2] = qz;
462}
463
464
465
466// ------------------------------------------------------------ //
467
468namespace internal
469{
470 namespace QIteratedImplementation
471 {
472 namespace
473 {
474 bool
475 uses_both_endpoints(const Quadrature<1> &base_quadrature)
476 {
477 const bool at_left =
478 std::any_of(base_quadrature.get_points().cbegin(),
479 base_quadrature.get_points().cend(),
480 [](const Point<1> &p) { return p == Point<1>{0.}; });
481 const bool at_right =
482 std::any_of(base_quadrature.get_points().cbegin(),
483 base_quadrature.get_points().cend(),
484 [](const Point<1> &p) { return p == Point<1>{1.}; });
485 return (at_left && at_right);
486 }
487
488 std::vector<Point<1>>
489 create_equidistant_interval_points(const unsigned int n_copies)
490 {
491 std::vector<Point<1>> support_points(n_copies + 1);
492
493 for (unsigned int copy = 0; copy < n_copies; ++copy)
494 support_points[copy][0] =
495 static_cast<double>(copy) / static_cast<double>(n_copies);
496
497 support_points[n_copies][0] = 1.0;
498
499 return support_points;
500 }
501 } // namespace
502 } // namespace QIteratedImplementation
503} // namespace internal
504
505
506
507template <>
508QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
509 : Quadrature<0>()
510{
512}
513
514
515
516template <>
517QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
518 : Quadrature<0>()
519{
521}
522
523
524
525template <>
527 const std::vector<Point<1>> &intervals)
528 : Quadrature<1>(
529 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
530 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
531 base_quadrature.size() * (intervals.size() - 1))
532{
533 Assert(base_quadrature.size() > 0, ExcNotInitialized());
534 Assert(intervals.size() > 1, ExcZero());
535
536 const unsigned int n_copies = intervals.size() - 1;
537
538 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
539 // we don't have to skip some points in order to get a reasonable quadrature
540 // formula
541 {
542 unsigned int next_point = 0;
543 for (unsigned int copy = 0; copy < n_copies; ++copy)
544 for (unsigned int q_point = 0; q_point < base_quadrature.size();
545 ++q_point)
546 {
547 this->quadrature_points[next_point] =
548 Point<1>(base_quadrature.point(q_point)[0] *
549 (intervals[copy + 1][0] - intervals[copy][0]) +
550 intervals[copy][0]);
551 this->weights[next_point] =
552 base_quadrature.weight(q_point) *
553 (intervals[copy + 1][0] - intervals[copy][0]);
554
555 ++next_point;
556 }
557 }
558 else
559 // skip doubly available points
560 {
561 const unsigned int left_index =
562 std::distance(base_quadrature.get_points().begin(),
563 std::find_if(base_quadrature.get_points().cbegin(),
564 base_quadrature.get_points().cend(),
565 [](const Point<1> &p) {
566 return p == Point<1>{0.};
567 }));
568
569 const unsigned int right_index =
570 std::distance(base_quadrature.get_points().begin(),
571 std::find_if(base_quadrature.get_points().cbegin(),
572 base_quadrature.get_points().cend(),
573 [](const Point<1> &p) {
574 return p == Point<1>{1.};
575 }));
576
577 const unsigned double_point_offset =
578 left_index + (base_quadrature.size() - right_index);
579
580 for (unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
581 for (unsigned int q_point = 0; q_point < base_quadrature.size();
582 ++q_point)
583 {
584 // skip the left point of this copy since we have already entered it
585 // the last time
586 if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
587 {
588 Assert(this->quadrature_points[next_point - double_point_offset]
589 .distance(Point<1>(
590 base_quadrature.point(q_point)[0] *
591 (intervals[copy + 1][0] - intervals[copy][0]) +
592 intervals[copy][0])) < 1e-10 /*tolerance*/,
594
595 this->weights[next_point - double_point_offset] +=
596 base_quadrature.weight(q_point) *
597 (intervals[copy + 1][0] - intervals[copy][0]);
598
599 continue;
600 }
601
602 this->quadrature_points[next_point] =
603 Point<1>(base_quadrature.point(q_point)[0] *
604 (intervals[copy + 1][0] - intervals[copy][0]) +
605 intervals[copy][0]);
606
607 // if this is the rightmost point of one of the non-last copies:
608 // give it the double weight
609 this->weights[next_point] =
610 base_quadrature.weight(q_point) *
611 (intervals[copy + 1][0] - intervals[copy][0]);
612
613 ++next_point;
614 }
615 }
616
617 // make sure that there is no rounding error for 0.0 and 1.0, since there
618 // are multiple asserts in the library checking for equality without
619 // tolerances
620 for (auto &i : this->quadrature_points)
621 if (std::abs(i[0] - 0.0) < 1e-12)
622 i[0] = 0.0;
623 else if (std::abs(i[0] - 1.0) < 1e-12)
624 i[0] = 1.0;
625
626#ifdef DEBUG
627 double sum_of_weights = 0;
628 for (unsigned int i = 0; i < this->size(); ++i)
629 sum_of_weights += this->weight(i);
630 Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
631#endif
632}
633
634
635
636template <>
638 const unsigned int n_copies)
639 : QIterated<1>(
640 base_quadrature,
641 internal::QIteratedImplementation::create_equidistant_interval_points(
642 n_copies))
643{
644 Assert(base_quadrature.size() > 0, ExcNotInitialized());
645 Assert(n_copies > 0, ExcZero());
646}
647
648
649
650// construct higher dimensional quadrature formula by tensor product
651// of lower dimensional iterated quadrature formulae
652template <int dim>
654 const std::vector<Point<1>> &intervals)
655 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, intervals),
656 QIterated<1>(base_quadrature, intervals))
657{}
658
659
660
661template <int dim>
663 const unsigned int n_copies)
664 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, n_copies),
665 QIterated<1>(base_quadrature, n_copies))
666{}
667
668
669
670// explicit instantiations; note: we need them all for all dimensions
671template class Quadrature<0>;
672template class Quadrature<1>;
673template class Quadrature<2>;
674template class Quadrature<3>;
675template class QAnisotropic<1>;
676template class QAnisotropic<2>;
677template class QAnisotropic<3>;
678template class QIterated<1>;
679template class QIterated<2>;
680template class QIterated<3>;
681
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
Definition point.h:111
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:353
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:374
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:368
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:359
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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)
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
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 > &)