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