Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
\(\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  ,
210  // quadrature_points(1),
211  weights(1, 1.)
212  , is_tensor_product_flag(false)
213 {}
214 
215 
216 template <>
218  : Subscriptor()
219 {
220  // this function should never be called -- this should be the copy constructor
221  // in 1d...
222  Assert(false, ExcImpossibleInDim(1));
223 }
224 #endif // DOXYGEN
225 
226 
227 
228 template <int dim>
230  : Subscriptor()
232  , weights(Utilities::fixed_power<dim>(q.size()))
233  , is_tensor_product_flag(true)
234 {
235  Assert(dim <= 3, ExcNotImplemented());
236 
237  const unsigned int n0 = q.size();
238  const unsigned int n1 = (dim > 1) ? n0 : 1;
239  const unsigned int n2 = (dim > 2) ? n0 : 1;
240 
241  unsigned int k = 0;
242  for (unsigned int i2 = 0; i2 < n2; ++i2)
243  for (unsigned int i1 = 0; i1 < n1; ++i1)
244  for (unsigned int i0 = 0; i0 < n0; ++i0)
245  {
246  quadrature_points[k](0) = q.point(i0)(0);
247  if (dim > 1)
248  quadrature_points[k](1) = q.point(i1)(0);
249  if (dim > 2)
250  quadrature_points[k](2) = q.point(i2)(0);
251  weights[k] = q.weight(i0);
252  if (dim > 1)
253  weights[k] *= q.weight(i1);
254  if (dim > 2)
255  weights[k] *= q.weight(i2);
256  ++k;
257  }
258 
259  tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
260  for (unsigned int i = 0; i < dim; ++i)
261  (*tensor_basis)[i] = q;
262 }
263 
264 
265 
266 template <int dim>
268  : Subscriptor()
270  , weights(q.weights)
272 {
273  if (dim > 1 && is_tensor_product_flag)
274  tensor_basis =
275  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
276 }
277 
278 
279 
280 template <int dim>
283 {
284  weights = q.weights;
287  if (dim > 1 && is_tensor_product_flag)
288  {
289  if (tensor_basis == nullptr)
290  tensor_basis =
291  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
292  else
294  }
295  return *this;
296 }
297 
298 
299 
300 template <int dim>
301 bool
303 {
304  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
305 }
306 
307 
308 
309 template <int dim>
310 std::size_t
312 {
315 }
316 
317 
318 
319 template <int dim>
320 typename std::conditional<dim == 1,
321  std::array<Quadrature<1>, dim>,
322  const std::array<Quadrature<1>, dim> &>::type
324 {
325  Assert(this->is_tensor_product_flag == true,
326  ExcMessage("This function only makes sense if "
327  "this object represents a tensor product!"));
328  Assert(tensor_basis != nullptr, ExcInternalError());
329 
330  return *tensor_basis;
331 }
332 
333 
334 #ifndef DOXYGEN
335 template <>
336 std::array<Quadrature<1>, 1>
338 {
339  Assert(this->is_tensor_product_flag == true,
340  ExcMessage("This function only makes sense if "
341  "this object represents a tensor product!"));
342 
343  return std::array<Quadrature<1>, 1>{{*this}};
344 }
345 #endif
346 
347 
348 
349 //---------------------------------------------------------------------------
350 template <int dim>
352  : Quadrature<dim>(qx.size())
353 {
354  Assert(dim == 1, ExcImpossibleInDim(dim));
355  unsigned int k = 0;
356  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
357  {
358  this->quadrature_points[k](0) = qx.point(k1)(0);
359  this->weights[k++] = qx.weight(k1);
360  }
361  Assert(k == this->size(), ExcInternalError());
362  this->is_tensor_product_flag = true;
363 }
364 
365 
366 
367 template <int dim>
369  const Quadrature<1> &qy)
370  : Quadrature<dim>(qx.size() * qy.size())
371 {
372  Assert(dim == 2, ExcImpossibleInDim(dim));
373 }
374 
375 
376 
377 template <>
379  : Quadrature<2>(qx.size() * qy.size())
380 {
381  unsigned int k = 0;
382  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
383  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
384  {
385  this->quadrature_points[k](0) = qx.point(k1)(0);
386  this->quadrature_points[k](1) = qy.point(k2)(0);
387  this->weights[k++] = qx.weight(k1) * qy.weight(k2);
388  }
389  Assert(k == this->size(), ExcInternalError());
390  this->is_tensor_product_flag = true;
391  const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
392  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 2>>(q_array);
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 
406 
407 
408 template <>
410  const Quadrature<1> &qy,
411  const Quadrature<1> &qz)
412  : Quadrature<3>(qx.size() * qy.size() * qz.size())
413 {
414  unsigned int k = 0;
415  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
416  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
417  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
418  {
419  this->quadrature_points[k](0) = qx.point(k1)(0);
420  this->quadrature_points[k](1) = qy.point(k2)(0);
421  this->quadrature_points[k](2) = qz.point(k3)(0);
422  this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
423  }
424  Assert(k == this->size(), ExcInternalError());
425  this->is_tensor_product_flag = true;
426  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
427  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 3>>(q_array);
428 }
429 
430 
431 
432 // ------------------------------------------------------------ //
433 
434 namespace internal
435 {
436  namespace QIteratedImplementation
437  {
438  namespace
439  {
440  bool
441  uses_both_endpoints(const Quadrature<1> &base_quadrature)
442  {
443  const bool at_left =
444  std::any_of(base_quadrature.get_points().cbegin(),
445  base_quadrature.get_points().cend(),
446  [](const Point<1> &p) { return p == Point<1>{0.}; });
447  const bool at_right =
448  std::any_of(base_quadrature.get_points().cbegin(),
449  base_quadrature.get_points().cend(),
450  [](const Point<1> &p) { return p == Point<1>{1.}; });
451  return (at_left && at_right);
452  }
453  } // namespace
454  } // namespace QIteratedImplementation
455 } // namespace internal
456 
457 
458 
459 template <>
460 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
461  : Quadrature<0>()
462 {
463  Assert(false, ExcNotImplemented());
464 }
465 
466 
467 
468 template <>
469 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
470  const unsigned int n_copies)
471  : Quadrature<1>(
472  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
473  (base_quadrature.size() - 1) * n_copies + 1 :
474  base_quadrature.size() * n_copies)
475 {
476  Assert(base_quadrature.size() > 0, ExcNotInitialized());
477  Assert(n_copies > 0, ExcZero());
478 
479  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
480  // we don't have to skip some points in order to get a reasonable quadrature
481  // formula
482  {
483  unsigned int next_point = 0;
484  for (unsigned int copy = 0; copy < n_copies; ++copy)
485  for (unsigned int q_point = 0; q_point < base_quadrature.size();
486  ++q_point)
487  {
488  this->quadrature_points[next_point] =
489  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
490  (1.0 * copy) / n_copies);
491  this->weights[next_point] =
492  base_quadrature.weight(q_point) / n_copies;
493 
494  ++next_point;
495  }
496  }
497  else
498  // skip doubly available points
499  {
500  unsigned int next_point = 0;
501 
502  // first find out the weights of the left and the right boundary points.
503  // note that these usually are but need not necessarily be the same
504  double double_point_weight = 0;
505  unsigned int n_end_points = 0;
506  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
507  // add up the weight if this is an endpoint
508  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
509  (base_quadrature.point(i) == Point<1>(1.0)))
510  {
511  double_point_weight += base_quadrature.weight(i);
512  ++n_end_points;
513  }
514  // scale the weight correctly
515  double_point_weight /= n_copies;
516 
517  // make sure the base quadrature formula has only one quadrature point per
518  // end point
519  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
520 
521 
522  for (unsigned int copy = 0; copy < n_copies; ++copy)
523  for (unsigned int q_point = 0; q_point < base_quadrature.size();
524  ++q_point)
525  {
526  // skip the left point of this copy since we have already entered it
527  // the last time
528  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
529  continue;
530 
531  this->quadrature_points[next_point] =
532  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
533  (1.0 * copy) / n_copies);
534 
535  // if this is the rightmost point of one of the non-last copies:
536  // give it the double weight
537  if ((copy != n_copies - 1) &&
538  (base_quadrature.point(q_point) == Point<1>(1.0)))
539  this->weights[next_point] = double_point_weight;
540  else
541  this->weights[next_point] =
542  base_quadrature.weight(q_point) / n_copies;
543 
544  ++next_point;
545  }
546  }
547 
548 #if DEBUG
549  double sum_of_weights = 0;
550  for (unsigned int i = 0; i < this->size(); ++i)
551  sum_of_weights += this->weight(i);
552  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
553 #endif
554 }
555 
556 
557 
558 // construct higher dimensional quadrature formula by tensor product
559 // of lower dimensional iterated quadrature formulae
560 template <int dim>
562  const unsigned int N)
563  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
564  QIterated<1>(base_quadrature, N))
565 {}
566 
567 
568 
569 // explicit instantiations; note: we need them all for all dimensions
570 template class Quadrature<0>;
571 template class Quadrature<1>;
572 template class Quadrature<2>;
573 template class Quadrature<3>;
574 template class QAnisotropic<1>;
575 template class QAnisotropic<2>;
576 template class QAnisotropic<3>;
577 template class QIterated<1>;
578 template class QIterated<2>;
579 template class QIterated<3>;
580 
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
std::vector< double > weights
Definition: quadrature.h:288
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=StaticMappingQ1< dim, spacedim >::mapping, const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:444
const std::vector< Point< dim > > & get_points() const
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:323
static ::ExceptionBase & ExcNotInitialized()
Definition: point.h:110
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:282
T fixed_power(const T t)
Definition: utilities.h:1045
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:351
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:1466
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
std::size_t memory_consumption() const
Definition: quadrature.cc:311
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:297
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:282
unsigned int size() const
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: cuda.h:32
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:302
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:561
static const char N
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:303
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()