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