Reference documentation for deal.II version Git 8ba89c874e 2021-08-04 22:06:33 -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  } // namespace
453  } // namespace QIteratedImplementation
454 } // namespace internal
455 
456 
457 
458 template <>
459 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
460  : Quadrature<0>()
461 {
462  Assert(false, ExcNotImplemented());
463 }
464 
465 
466 
467 template <>
468 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
469  const unsigned int n_copies)
470  : Quadrature<1>(
471  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
472  (base_quadrature.size() - 1) * n_copies + 1 :
473  base_quadrature.size() * n_copies)
474 {
475  Assert(base_quadrature.size() > 0, ExcNotInitialized());
476  Assert(n_copies > 0, ExcZero());
477 
478  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
479  // we don't have to skip some points in order to get a reasonable quadrature
480  // formula
481  {
482  unsigned int next_point = 0;
483  for (unsigned int copy = 0; copy < n_copies; ++copy)
484  for (unsigned int q_point = 0; q_point < base_quadrature.size();
485  ++q_point)
486  {
487  this->quadrature_points[next_point] =
488  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
489  (1.0 * copy) / n_copies);
490  this->weights[next_point] =
491  base_quadrature.weight(q_point) / n_copies;
492 
493  ++next_point;
494  }
495  }
496  else
497  // skip doubly available points
498  {
499  unsigned int next_point = 0;
500 
501  // first find out the weights of the left and the right boundary points.
502  // note that these usually are but need not necessarily be the same
503  double double_point_weight = 0;
504  unsigned int n_end_points = 0;
505  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
506  // add up the weight if this is an endpoint
507  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
508  (base_quadrature.point(i) == Point<1>(1.0)))
509  {
510  double_point_weight += base_quadrature.weight(i);
511  ++n_end_points;
512  }
513  // scale the weight correctly
514  double_point_weight /= n_copies;
515 
516  // make sure the base quadrature formula has only one quadrature point per
517  // end point
518  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
519 
520 
521  for (unsigned int copy = 0; copy < n_copies; ++copy)
522  for (unsigned int q_point = 0; q_point < base_quadrature.size();
523  ++q_point)
524  {
525  // skip the left point of this copy since we have already entered it
526  // the last time
527  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
528  continue;
529 
530  this->quadrature_points[next_point] =
531  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
532  (1.0 * copy) / n_copies);
533 
534  // if this is the rightmost point of one of the non-last copies:
535  // give it the double weight
536  if ((copy != n_copies - 1) &&
537  (base_quadrature.point(q_point) == Point<1>(1.0)))
538  this->weights[next_point] = double_point_weight;
539  else
540  this->weights[next_point] =
541  base_quadrature.weight(q_point) / n_copies;
542 
543  ++next_point;
544  }
545  }
546 
547 #if DEBUG
548  double sum_of_weights = 0;
549  for (unsigned int i = 0; i < this->size(); ++i)
550  sum_of_weights += this->weight(i);
551  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
552 #endif
553 }
554 
555 
556 
557 // construct higher dimensional quadrature formula by tensor product
558 // of lower dimensional iterated quadrature formulae
559 template <int dim>
561  const unsigned int N)
562  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
563  QIterated<1>(base_quadrature, N))
564 {}
565 
566 
567 
568 // explicit instantiations; note: we need them all for all dimensions
569 template class Quadrature<0>;
570 template class Quadrature<1>;
571 template class Quadrature<2>;
572 template class Quadrature<3>;
573 template class QAnisotropic<1>;
574 template class QAnisotropic<2>;
575 template class QAnisotropic<3>;
576 template class QIterated<1>;
577 template class QIterated<2>;
578 template class QIterated<3>;
579 
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1659
std::vector< double > weights
Definition: quadrature.h:292
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:1081
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:1467
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
static ::ExceptionBase & ExcInvalidQuadratureFormula()
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:301
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:286
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:560
static const char N
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:307
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()