Reference documentation for deal.II version GIT ffb4c3937f 2023-03-31 14: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 
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 <>
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  // placate compiler in the dim == 1 case
377  constexpr int dim_1 = dim == 2 ? 1 : 0;
378 
379  unsigned int k = 0;
380  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
381  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
382  {
383  this->quadrature_points[k](0) = qx.point(k1)(0);
384  this->quadrature_points[k](dim_1) = qy.point(k2)(0);
385  this->weights[k++] = qx.weight(k1) * qy.weight(k2);
386  }
387  Assert(k == this->size(), ExcInternalError());
388 
389  this->is_tensor_product_flag = true;
390  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
391  (*this->tensor_basis)[0] = qx;
392  (*this->tensor_basis)[dim_1] = qy;
393 }
394 
395 
396 
397 template <int dim>
399  const Quadrature<1> &qy,
400  const Quadrature<1> &qz)
401  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
402 {
403  Assert(dim == 3, ExcImpossibleInDim(dim));
404 
405  // placate compiler in lower dimensions
406  constexpr int dim_1 = dim == 3 ? 1 : 0;
407  constexpr int dim_2 = dim == 3 ? 2 : 0;
408 
409  unsigned int k = 0;
410  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
411  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
412  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
413  {
414  this->quadrature_points[k](0) = qx.point(k1)(0);
415  this->quadrature_points[k](dim_1) = qy.point(k2)(0);
416  this->quadrature_points[k](dim_2) = qz.point(k3)(0);
417  this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
418  }
419  Assert(k == this->size(), ExcInternalError());
420 
421  this->is_tensor_product_flag = true;
422  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
423  (*this->tensor_basis)[0] = qx;
424  (*this->tensor_basis)[dim_1] = qy;
425  (*this->tensor_basis)[dim_2] = qz;
426 }
427 
428 
429 
430 // ------------------------------------------------------------ //
431 
432 namespace internal
433 {
434  namespace QIteratedImplementation
435  {
436  namespace
437  {
438  bool
439  uses_both_endpoints(const Quadrature<1> &base_quadrature)
440  {
441  const bool at_left =
442  std::any_of(base_quadrature.get_points().cbegin(),
443  base_quadrature.get_points().cend(),
444  [](const Point<1> &p) { return p == Point<1>{0.}; });
445  const bool at_right =
446  std::any_of(base_quadrature.get_points().cbegin(),
447  base_quadrature.get_points().cend(),
448  [](const Point<1> &p) { return p == Point<1>{1.}; });
449  return (at_left && at_right);
450  }
451 
452  std::vector<Point<1>>
453  create_equidistant_interval_points(const unsigned int n_copies)
454  {
455  std::vector<Point<1>> support_points(n_copies + 1);
456 
457  for (unsigned int copy = 0; copy < n_copies; ++copy)
458  support_points[copy][0] =
459  static_cast<double>(copy) / static_cast<double>(n_copies);
460 
461  support_points[n_copies][0] = 1.0;
462 
463  return support_points;
464  }
465  } // namespace
466  } // namespace QIteratedImplementation
467 } // namespace internal
468 
469 
470 
471 template <>
472 QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
473  : Quadrature<0>()
474 {
475  Assert(false, ExcNotImplemented());
476 }
477 
478 
479 
480 template <>
481 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
482  : Quadrature<0>()
483 {
484  Assert(false, ExcNotImplemented());
485 }
486 
487 
488 
489 template <>
490 QIterated<1>::QIterated(const Quadrature<1> & base_quadrature,
491  const std::vector<Point<1>> &intervals)
492  : Quadrature<1>(
493  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
494  (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
495  base_quadrature.size() * (intervals.size() - 1))
496 {
497  Assert(base_quadrature.size() > 0, ExcNotInitialized());
498  Assert(intervals.size() > 1, ExcZero());
499 
500  const unsigned int n_copies = intervals.size() - 1;
501 
502  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
503  // we don't have to skip some points in order to get a reasonable quadrature
504  // formula
505  {
506  unsigned int next_point = 0;
507  for (unsigned int copy = 0; copy < n_copies; ++copy)
508  for (unsigned int q_point = 0; q_point < base_quadrature.size();
509  ++q_point)
510  {
511  this->quadrature_points[next_point] =
512  Point<1>(base_quadrature.point(q_point)(0) *
513  (intervals[copy + 1][0] - intervals[copy][0]) +
514  intervals[copy][0]);
515  this->weights[next_point] =
516  base_quadrature.weight(q_point) *
517  (intervals[copy + 1][0] - intervals[copy][0]);
518 
519  ++next_point;
520  }
521  }
522  else
523  // skip doubly available points
524  {
525  const unsigned int left_index =
526  std::distance(base_quadrature.get_points().begin(),
527  std::find_if(base_quadrature.get_points().cbegin(),
528  base_quadrature.get_points().cend(),
529  [](const Point<1> &p) {
530  return p == Point<1>{0.};
531  }));
532 
533  const unsigned int right_index =
534  std::distance(base_quadrature.get_points().begin(),
535  std::find_if(base_quadrature.get_points().cbegin(),
536  base_quadrature.get_points().cend(),
537  [](const Point<1> &p) {
538  return p == Point<1>{1.};
539  }));
540 
541  const unsigned double_point_offset =
542  left_index + (base_quadrature.size() - right_index);
543 
544  for (unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
545  for (unsigned int q_point = 0; q_point < base_quadrature.size();
546  ++q_point)
547  {
548  // skip the left point of this copy since we have already entered it
549  // the last time
550  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
551  {
552  Assert(this->quadrature_points[next_point - double_point_offset]
553  .distance(Point<1>(
554  base_quadrature.point(q_point)(0) *
555  (intervals[copy + 1][0] - intervals[copy][0]) +
556  intervals[copy][0])) < 1e-10 /*tolerance*/,
557  ExcInternalError());
558 
559  this->weights[next_point - double_point_offset] +=
560  base_quadrature.weight(q_point) *
561  (intervals[copy + 1][0] - intervals[copy][0]);
562 
563  continue;
564  }
565 
566  this->quadrature_points[next_point] =
567  Point<1>(base_quadrature.point(q_point)(0) *
568  (intervals[copy + 1][0] - intervals[copy][0]) +
569  intervals[copy][0]);
570 
571  // if this is the rightmost point of one of the non-last copies:
572  // give it the double weight
573  this->weights[next_point] =
574  base_quadrature.weight(q_point) *
575  (intervals[copy + 1][0] - intervals[copy][0]);
576 
577  ++next_point;
578  }
579  }
580 
581  // make sure that there is no rounding error for 0.0 and 1.0, since there
582  // are multiple asserts in the library checking for equality without
583  // tolerances
584  for (auto &i : this->quadrature_points)
585  if (std::abs(i[0] - 0.0) < 1e-12)
586  i[0] = 0.0;
587  else if (std::abs(i[0] - 1.0) < 1e-12)
588  i[0] = 1.0;
589 
590 #ifdef DEBUG
591  double sum_of_weights = 0;
592  for (unsigned int i = 0; i < this->size(); ++i)
593  sum_of_weights += this->weight(i);
594  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
595 #endif
596 }
597 
598 
599 
600 template <>
601 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
602  const unsigned int n_copies)
603  : QIterated<1>(
604  base_quadrature,
605  internal::QIteratedImplementation::create_equidistant_interval_points(
606  n_copies))
607 {
608  Assert(base_quadrature.size() > 0, ExcNotInitialized());
609  Assert(n_copies > 0, ExcZero());
610 }
611 
612 
613 
614 // construct higher dimensional quadrature formula by tensor product
615 // of lower dimensional iterated quadrature formulae
616 template <int dim>
617 QIterated<dim>::QIterated(const Quadrature<1> & base_quadrature,
618  const std::vector<Point<1>> &intervals)
619  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, intervals),
620  QIterated<1>(base_quadrature, intervals))
621 {}
622 
623 
624 
625 template <int dim>
627  const unsigned int n_copies)
628  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, n_copies),
629  QIterated<1>(base_quadrature, n_copies))
630 {}
631 
632 
633 
634 // explicit instantiations; note: we need them all for all dimensions
635 template class Quadrature<0>;
636 template class Quadrature<1>;
637 template class Quadrature<2>;
638 template class Quadrature<3>;
639 template class QAnisotropic<1>;
640 template class QAnisotropic<2>;
641 template class QAnisotropic<3>;
642 template class QIterated<1>;
643 template class QIterated<2>;
644 template class QIterated<3>;
645 
Definition: point.h:110
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:353
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:626
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:326
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:347
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:341
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:332
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
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:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
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:189
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:966
void copy(const T *begin, const T *end, U *dest)