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