Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00: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\}}\)
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 
18 #include <deal.II/base/utilities.h>
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
31 template <>
32 Quadrature<0>::Quadrature(const unsigned int n_q)
33  : quadrature_points(n_q)
34  , weights(n_q, 0)
35  , is_tensor_product_flag(false)
36 {}
37 #endif
38 
39 
40 
41 template <int dim>
42 Quadrature<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 
50 template <int dim>
51 void
53  const std::vector<double> &w)
54 {
55  AssertDimension(w.size(), p.size());
57  weights = w;
58  is_tensor_product_flag = dim == 1;
59 }
60 
61 
62 
63 template <int dim>
64 Quadrature<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 
76 template <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 
89 template <int dim>
90 Quadrature<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 
101 template <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
118 template <>
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 
127 template <>
129  : is_tensor_product_flag(false)
130 {
131  Assert(false, ExcImpossibleInDim(0));
132 }
133 
134 
135 
136 template <>
137 Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
138 {
139  Assert(false, ExcImpossibleInDim(0));
140 }
141 #endif // DOXYGEN
142 
143 
144 
145 template <int dim>
147  : quadrature_points(q1.size() * q2.size())
148  , weights(q1.size() * q2.size())
149  , is_tensor_product_flag(q1.is_tensor_product())
150 {
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
189 template <>
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)
197  {
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  }
206 
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 
222 template <>
225  , quadrature_points(1)
226  , weights(1, 1.)
227  , is_tensor_product_flag(false)
228 {}
229 
230 
231 template <>
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 
243 template <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 
281 template <int dim>
283  : Subscriptor()
285  , weights(q.weights)
286  , is_tensor_product_flag(q.is_tensor_product_flag)
287 {
288  if (dim > 1 && is_tensor_product_flag)
289  tensor_basis =
290  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
291 }
292 
293 
294 
295 template <int dim>
298 {
299  weights = q.weights;
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 
315 template <int dim>
316 bool
318 {
319  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
320 }
321 
322 
323 
324 template <int dim>
325 std::size_t
327 {
330 }
331 
332 
333 
334 template <int dim>
335 typename 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
350 template <>
351 std::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 //---------------------------------------------------------------------------
365 template <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 
382 template <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 
410 template <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 
445 namespace 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 
484 template <>
485 QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
486  : Quadrature<0>()
487 {
488  Assert(false, ExcNotImplemented());
489 }
490 
491 
492 
493 template <>
494 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
495  : Quadrature<0>()
496 {
497  Assert(false, ExcNotImplemented());
498 }
499 
500 
501 
502 template <>
503 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
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*/,
570  ExcInternalError());
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 
613 template <>
614 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
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
629 template <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 
638 template <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
648 template class Quadrature<0>;
649 template class Quadrature<1>;
650 template class Quadrature<2>;
651 template class Quadrature<3>;
652 template class QAnisotropic<1>;
653 template class QAnisotropic<2>;
654 template class QAnisotropic<3>;
655 template class QIterated<1>;
656 template class QIterated<2>;
657 template class QIterated<3>;
658 
Definition: point.h:112
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:366
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:639
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:339
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:360
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:297
const std::vector< Point< dim > > & get_points() const
std::size_t memory_consumption() const
Definition: quadrature.cc:326
void initialize(const std::vector< Point< dim >> &points, const std::vector< double > &weights)
Definition: quadrature.cc:52
bool is_tensor_product_flag
Definition: quadrature.h:354
Quadrature(const unsigned int n_quadrature_points=0)
Definition: quadrature.cc:42
double weight(const unsigned int i) const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:338
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:317
const Point< dim > & point(const unsigned int i) const
std::vector< double > weights
Definition: quadrature.h:345
unsigned int size() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression fabs(const Expression &x)
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:192
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={})
Definition: generators.cc:490
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
T sum(const T &t, const MPI_Comm mpi_communicator)
T fixed_power(const T t)
Definition: utilities.h:975
void copy(const T *begin, const T *end, U *dest)