Reference documentation for deal.II version GIT relicensing-437-g81ec864850 2024-04-19 07:30:02+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 <>
31Quadrature<0>::Quadrature(const unsigned int n_q)
33 , weights(n_q, 0)
34 , is_tensor_product_flag(false)
35{}
36#endif
37
38
39
40template <int dim>
41Quadrature<dim>::Quadrature(const unsigned int n_q)
42 : quadrature_points(n_q, Point<dim>())
43 , weights(n_q, 0)
44 , is_tensor_product_flag(dim == 1)
45{}
46
47
48
49template <int dim>
50void
52 const ArrayView<const double> &weights)
53{
54 this->weights.clear();
55 if (weights.size() > 0)
56 {
57 AssertDimension(weights.size(), points.size());
58 this->weights.insert(this->weights.end(), weights.begin(), weights.end());
59 }
60 else
61 this->weights.resize(points.size(),
62 std::numeric_limits<double>::infinity());
63
64 quadrature_points.clear();
65 quadrature_points.insert(quadrature_points.end(),
66 points.begin(),
67 points.end());
68
69 is_tensor_product_flag = dim == 1;
70}
71
72
73
74template <int dim>
75Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points,
76 const std::vector<double> &weights)
77 : quadrature_points(points)
78 , weights(weights)
79 , is_tensor_product_flag(dim == 1)
80{
81 Assert(weights.size() == points.size(),
82 ExcDimensionMismatch(weights.size(), points.size()));
83}
84
85
86
87template <int dim>
89 std::vector<double> &&weights)
90 : quadrature_points(std::move(points))
91 , weights(std::move(weights))
92 , is_tensor_product_flag(dim == 1)
93{
94 Assert(weights.size() == points.size(),
95 ExcDimensionMismatch(weights.size(), points.size()));
96}
97
98
99
100template <int dim>
101Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points)
102 : quadrature_points(points)
103 , weights(points.size(), std::numeric_limits<double>::infinity())
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>
114 : quadrature_points(std::vector<Point<dim>>(1, point))
115 , weights(std::vector<double>(1, 1.))
116 , is_tensor_product_flag(true)
117 , tensor_basis(new std::array<Quadrature<1>, dim>())
118{
119 for (unsigned int i = 0; i < dim; ++i)
120 {
121 const std::vector<Point<1>> quad_vec_1d(1, Point<1>(point[i]));
122 (*tensor_basis)[i] = Quadrature<1>(quad_vec_1d, weights);
123 }
124}
125
126
127
128#ifndef DOXYGEN
129template <>
131 : quadrature_points(std::vector<Point<1>>(1, point))
132 , weights(std::vector<double>(1, 1.))
133 , is_tensor_product_flag(true)
134{}
135
136
137
138template <>
140 : is_tensor_product_flag(false)
141{
142 Assert(false, ExcImpossibleInDim(0));
143}
144
145
146
147template <>
148Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
149{
151}
152#endif // DOXYGEN
153
154
155
156template <int dim>
158 : quadrature_points(q1.size() * q2.size())
159 , weights(q1.size() * q2.size())
160 , is_tensor_product_flag(q1.is_tensor_product())
161{
162 unsigned int present_index = 0;
163 for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
164 for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
165 {
166 // compose coordinates of new quadrature point by tensor product in the
167 // last component
168 for (unsigned int d = 0; d < dim - 1; ++d)
169 quadrature_points[present_index][d] = q1.point(i1)[d];
170 quadrature_points[present_index][dim - 1] = q2.point(i2)[0];
172 weights[present_index] = q1.weight(i1) * q2.weight(i2);
173
174 ++present_index;
175 }
177#ifdef DEBUG
178 if (size() > 0)
179 {
180 double sum = 0;
181 for (unsigned int i = 0; i < size(); ++i)
182 sum += weights[i];
183 // we cannot guarantee the sum of weights to be exactly one, but it should
184 // be near that.
185 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
186 }
187#endif
188
190 {
191 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
192 for (unsigned int i = 0; i < dim - 1; ++i)
193 (*tensor_basis)[i] = q1.get_tensor_basis()[i];
194 (*tensor_basis)[dim - 1] = q2;
195 }
196}
198
199#ifndef DOXYGEN
200template <>
202 : quadrature_points(q2.size())
203 , weights(q2.size())
204 , is_tensor_product_flag(true)
205{
206 unsigned int present_index = 0;
207 for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
208 {
209 // compose coordinates of new quadrature point by tensor product in the
210 // last component
211 quadrature_points[present_index][0] = q2.point(i2)[0];
213 weights[present_index] = q2.weight(i2);
214
215 ++present_index;
216 }
217
218# ifdef DEBUG
219 if (size() > 0)
220 {
221 double sum = 0;
222 for (unsigned int i = 0; i < size(); ++i)
223 sum += weights[i];
224 // we cannot guarantee the sum of weights to be exactly one, but it should
225 // be near that.
226 Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
227 }
228# endif
229}
230
231
232
233template <>
235 : Subscriptor()
237 , weights(1, 1.)
238 , is_tensor_product_flag(false)
239{}
240
241
242template <>
244 : Subscriptor()
245{
246 // this function should never be called -- this should be the copy constructor
247 // in 1d...
249}
250#endif // DOXYGEN
251
252
253
254template <int dim>
256 : Subscriptor()
257 , quadrature_points(Utilities::fixed_power<dim>(q.size()))
258 , weights(Utilities::fixed_power<dim>(q.size()))
259 , is_tensor_product_flag(true)
260{
261 Assert(dim <= 3, ExcNotImplemented());
262
263 const unsigned int n0 = q.size();
264 const unsigned int n1 = (dim > 1) ? n0 : 1;
265 const unsigned int n2 = (dim > 2) ? n0 : 1;
266
267 unsigned int k = 0;
268 for (unsigned int i2 = 0; i2 < n2; ++i2)
269 for (unsigned int i1 = 0; i1 < n1; ++i1)
270 for (unsigned int i0 = 0; i0 < n0; ++i0)
271 {
272 quadrature_points[k][0] = q.point(i0)[0];
273 if (dim > 1)
274 quadrature_points[k][1] = q.point(i1)[0];
275 if (dim > 2)
276 quadrature_points[k][2] = q.point(i2)[0];
277 weights[k] = q.weight(i0);
278 if (dim > 1)
279 weights[k] *= q.weight(i1);
280 if (dim > 2)
281 weights[k] *= q.weight(i2);
282 ++k;
283 }
284
285 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
286 for (unsigned int i = 0; i < dim; ++i)
287 (*tensor_basis)[i] = q;
288}
289
290
291
292template <int dim>
294 : Subscriptor()
295 , quadrature_points(q.quadrature_points)
296 , weights(q.weights)
297 , is_tensor_product_flag(q.is_tensor_product_flag)
298{
299 if (dim > 1 && is_tensor_product_flag)
301 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
302}
303
304
305
306template <int dim>
309{
310 weights = q.weights;
311 quadrature_points = q.quadrature_points;
312 is_tensor_product_flag = q.is_tensor_product_flag;
313 if (dim > 1 && is_tensor_product_flag)
314 {
315 if (tensor_basis == nullptr)
316 tensor_basis =
317 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
318 else
319 *tensor_basis = *q.tensor_basis;
320 }
321 return *this;
322}
323
324
325
326template <int dim>
327bool
329{
330 return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
331}
332
333
334
335template <int dim>
336std::size_t
342
343
344
345template <int dim>
346typename std::conditional_t<dim == 1,
347 std::array<Quadrature<1>, dim>,
348 const std::array<Quadrature<1>, dim> &>
350{
351 Assert(this->is_tensor_product_flag == true,
352 ExcMessage("This function only makes sense if "
353 "this object represents a tensor product!"));
354 Assert(tensor_basis != nullptr, ExcInternalError());
355
356 return *tensor_basis;
357}
358
359
360#ifndef DOXYGEN
361template <>
362std::array<Quadrature<1>, 1>
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
369 return std::array<Quadrature<1>, 1>{{*this}};
370}
371#endif
372
373
374
375//---------------------------------------------------------------------------
376template <int dim>
378 : Quadrature<dim>(qx.size())
379{
380 Assert(dim == 1, ExcImpossibleInDim(dim));
381 unsigned int k = 0;
382 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
383 {
384 this->quadrature_points[k][0] = qx.point(k1)[0];
385 this->weights[k++] = qx.weight(k1);
386 }
387 Assert(k == this->size(), ExcInternalError());
388 this->is_tensor_product_flag = true;
389}
390
391
392
393template <int dim>
395 const Quadrature<1> &qy)
396 : Quadrature<dim>(qx.size() * qy.size())
397{
398 Assert(dim == 2, ExcImpossibleInDim(dim));
399
400 // placate compiler in the dim == 1 case
401 constexpr int dim_1 = dim == 2 ? 1 : 0;
402
403 unsigned int k = 0;
404 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
405 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
406 {
407 this->quadrature_points[k][0] = qx.point(k1)[0];
408 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
409 this->weights[k++] = qx.weight(k1) * qy.weight(k2);
410 }
411 Assert(k == this->size(), ExcInternalError());
412
413 this->is_tensor_product_flag = true;
414 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
415 (*this->tensor_basis)[0] = qx;
416 (*this->tensor_basis)[dim_1] = qy;
417}
418
419
420
421template <int dim>
423 const Quadrature<1> &qy,
424 const Quadrature<1> &qz)
425 : Quadrature<dim>(qx.size() * qy.size() * qz.size())
426{
427 Assert(dim == 3, ExcImpossibleInDim(dim));
428
429 // placate compiler in lower dimensions
430 constexpr int dim_1 = dim == 3 ? 1 : 0;
431 constexpr int dim_2 = dim == 3 ? 2 : 0;
432
433 unsigned int k = 0;
434 for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
435 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
436 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
437 {
438 this->quadrature_points[k][0] = qx.point(k1)[0];
439 this->quadrature_points[k][dim_1] = qy.point(k2)[0];
440 this->quadrature_points[k][dim_2] = qz.point(k3)[0];
441 this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
442 }
443 Assert(k == this->size(), ExcInternalError());
444
445 this->is_tensor_product_flag = true;
446 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
447 (*this->tensor_basis)[0] = qx;
448 (*this->tensor_basis)[dim_1] = qy;
449 (*this->tensor_basis)[dim_2] = qz;
450}
451
452
453
454// ------------------------------------------------------------ //
455
456namespace internal
457{
458 namespace QIteratedImplementation
459 {
460 namespace
461 {
462 bool
463 uses_both_endpoints(const Quadrature<1> &base_quadrature)
464 {
465 const bool at_left =
466 std::any_of(base_quadrature.get_points().cbegin(),
467 base_quadrature.get_points().cend(),
468 [](const Point<1> &p) { return p == Point<1>{0.}; });
469 const bool at_right =
470 std::any_of(base_quadrature.get_points().cbegin(),
471 base_quadrature.get_points().cend(),
472 [](const Point<1> &p) { return p == Point<1>{1.}; });
473 return (at_left && at_right);
474 }
475
476 std::vector<Point<1>>
477 create_equidistant_interval_points(const unsigned int n_copies)
478 {
479 std::vector<Point<1>> support_points(n_copies + 1);
480
481 for (unsigned int copy = 0; copy < n_copies; ++copy)
482 support_points[copy][0] =
483 static_cast<double>(copy) / static_cast<double>(n_copies);
484
485 support_points[n_copies][0] = 1.0;
486
487 return support_points;
488 }
489 } // namespace
490 } // namespace QIteratedImplementation
491} // namespace internal
492
493
494
495template <>
496QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
497 : Quadrature<0>()
498{
500}
501
502
503
504template <>
505QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
506 : Quadrature<0>()
507{
509}
510
511
512
513template <>
515 const std::vector<Point<1>> &intervals)
516 : Quadrature<1>(
517 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
518 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
519 base_quadrature.size() * (intervals.size() - 1))
520{
521 Assert(base_quadrature.size() > 0, ExcNotInitialized());
522 Assert(intervals.size() > 1, ExcZero());
523
524 const unsigned int n_copies = intervals.size() - 1;
525
526 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
527 // we don't have to skip some points in order to get a reasonable quadrature
528 // formula
529 {
530 unsigned int next_point = 0;
531 for (unsigned int copy = 0; copy < n_copies; ++copy)
532 for (unsigned int q_point = 0; q_point < base_quadrature.size();
533 ++q_point)
534 {
535 this->quadrature_points[next_point] =
536 Point<1>(base_quadrature.point(q_point)[0] *
537 (intervals[copy + 1][0] - intervals[copy][0]) +
538 intervals[copy][0]);
539 this->weights[next_point] =
540 base_quadrature.weight(q_point) *
541 (intervals[copy + 1][0] - intervals[copy][0]);
542
543 ++next_point;
544 }
545 }
546 else
547 // skip doubly available points
548 {
549 const unsigned int left_index =
550 std::distance(base_quadrature.get_points().begin(),
551 std::find_if(base_quadrature.get_points().cbegin(),
552 base_quadrature.get_points().cend(),
553 [](const Point<1> &p) {
554 return p == Point<1>{0.};
555 }));
556
557 const unsigned int right_index =
558 std::distance(base_quadrature.get_points().begin(),
559 std::find_if(base_quadrature.get_points().cbegin(),
560 base_quadrature.get_points().cend(),
561 [](const Point<1> &p) {
562 return p == Point<1>{1.};
563 }));
564
565 const unsigned double_point_offset =
566 left_index + (base_quadrature.size() - right_index);
567
568 for (unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
569 for (unsigned int q_point = 0; q_point < base_quadrature.size();
570 ++q_point)
571 {
572 // skip the left point of this copy since we have already entered it
573 // the last time
574 if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
575 {
576 Assert(this->quadrature_points[next_point - double_point_offset]
577 .distance(Point<1>(
578 base_quadrature.point(q_point)[0] *
579 (intervals[copy + 1][0] - intervals[copy][0]) +
580 intervals[copy][0])) < 1e-10 /*tolerance*/,
582
583 this->weights[next_point - double_point_offset] +=
584 base_quadrature.weight(q_point) *
585 (intervals[copy + 1][0] - intervals[copy][0]);
586
587 continue;
588 }
589
590 this->quadrature_points[next_point] =
591 Point<1>(base_quadrature.point(q_point)[0] *
592 (intervals[copy + 1][0] - intervals[copy][0]) +
593 intervals[copy][0]);
594
595 // if this is the rightmost point of one of the non-last copies:
596 // give it the double weight
597 this->weights[next_point] =
598 base_quadrature.weight(q_point) *
599 (intervals[copy + 1][0] - intervals[copy][0]);
600
601 ++next_point;
602 }
603 }
604
605 // make sure that there is no rounding error for 0.0 and 1.0, since there
606 // are multiple asserts in the library checking for equality without
607 // tolerances
608 for (auto &i : this->quadrature_points)
609 if (std::abs(i[0] - 0.0) < 1e-12)
610 i[0] = 0.0;
611 else if (std::abs(i[0] - 1.0) < 1e-12)
612 i[0] = 1.0;
613
614#ifdef DEBUG
615 double sum_of_weights = 0;
616 for (unsigned int i = 0; i < this->size(); ++i)
617 sum_of_weights += this->weight(i);
618 Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
619#endif
620}
621
622
623
624template <>
626 const unsigned int n_copies)
627 : QIterated<1>(
628 base_quadrature,
629 internal::QIteratedImplementation::create_equidistant_interval_points(
630 n_copies))
631{
632 Assert(base_quadrature.size() > 0, ExcNotInitialized());
633 Assert(n_copies > 0, ExcZero());
634}
635
636
637
638// construct higher dimensional quadrature formula by tensor product
639// of lower dimensional iterated quadrature formulae
640template <int dim>
642 const std::vector<Point<1>> &intervals)
643 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, intervals),
644 QIterated<1>(base_quadrature, intervals))
645{}
646
647
648
649template <int dim>
651 const unsigned int n_copies)
652 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, n_copies),
653 QIterated<1>(base_quadrature, n_copies))
654{}
655
656
657
658// explicit instantiations; note: we need them all for all dimensions
659template class Quadrature<0>;
660template class Quadrature<1>;
661template class Quadrature<2>;
662template class Quadrature<3>;
663template class QAnisotropic<1>;
664template class QAnisotropic<2>;
665template class QAnisotropic<3>;
666template class QIterated<1>;
667template class QIterated<2>;
668template class QIterated<3>;
669
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:343
void initialize(const ArrayView< const Point< dim > > &points, const ArrayView< const double > &weights={})
Definition quadrature.cc:51
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition quadrature.h:364
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:358
Quadrature(const unsigned int n_quadrature_points=0)
Definition quadrature.cc:41
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:349
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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()
if(marked_vertices.size() !=0) for(auto it
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 > &)