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