Reference documentation for deal.II version Git 409ee4b167 2020-08-14 09:46:12 -0400
\(\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 
20 #include <deal.II/base/utilities.h>
21 
22 #include <cmath>
23 #include <cstdlib>
24 #include <iterator>
25 #include <memory>
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 #ifndef DOXYGEN
101 template <>
103  : quadrature_points(std::vector<Point<1>>(1, point))
104  , weights(std::vector<double>(1, 1.))
105  , is_tensor_product_flag(true)
106 {}
107 
108 
109 
110 template <>
112 {
113  Assert(false, ExcImpossibleInDim(0));
114 }
115 #endif // DOXYGEN
116 
117 
118 
119 template <int dim>
121  : quadrature_points(q1.size() * q2.size())
122  , weights(q1.size() * q2.size())
124 {
125  unsigned int present_index = 0;
126  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
127  for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
128  {
129  // compose coordinates of
130  // new quadrature point by tensor
131  // product in the last component
132  for (unsigned int d = 0; d < dim - 1; ++d)
133  quadrature_points[present_index](d) = q1.point(i1)(d);
134  quadrature_points[present_index](dim - 1) = q2.point(i2)(0);
135 
136  weights[present_index] = q1.weight(i1) * q2.weight(i2);
137 
138  ++present_index;
139  }
140 
141 #ifdef DEBUG
142  if (size() > 0)
143  {
144  double sum = 0;
145  for (unsigned int i = 0; i < size(); ++i)
146  sum += weights[i];
147  // we cannot guarantee the sum of weights
148  // to be exactly one, but it should be
149  // near that.
150  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
151  }
152 #endif
153 
155  {
156  tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
157  for (unsigned int i = 0; i < dim - 1; ++i)
158  (*tensor_basis)[i] = q1.get_tensor_basis()[i];
159  (*tensor_basis)[dim - 1] = q2;
160  }
161 }
162 
163 
164 #ifndef DOXYGEN
165 template <>
167  : quadrature_points(q2.size())
168  , weights(q2.size())
169  , is_tensor_product_flag(true)
170 {
171  unsigned int present_index = 0;
172  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
173  {
174  // compose coordinates of
175  // new quadrature point by tensor
176  // product in the last component
177  quadrature_points[present_index](0) = q2.point(i2)(0);
178 
179  weights[present_index] = q2.weight(i2);
180 
181  ++present_index;
182  }
183 
184 # ifdef DEBUG
185  if (size() > 0)
186  {
187  double sum = 0;
188  for (unsigned int i = 0; i < size(); ++i)
189  sum += weights[i];
190  // we cannot guarantee the sum of weights
191  // to be exactly one, but it should be
192  // near that.
193  Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
194  }
195 # endif
196 }
197 
198 
199 
200 template <>
202  : Subscriptor()
203  ,
204  // quadrature_points(1),
205  weights(1, 1.)
206  , is_tensor_product_flag(false)
207 {}
208 
209 
210 template <>
212  : Subscriptor()
213 {
214  // this function should never be
215  // called -- this should be the
216  // copy constructor in 1d...
217  Assert(false, ExcImpossibleInDim(1));
218 }
219 #endif // DOXYGEN
220 
221 
222 
223 template <int dim>
225  : Subscriptor()
227  , weights(Utilities::fixed_power<dim>(q.size()))
228  , is_tensor_product_flag(true)
229 {
230  Assert(dim <= 3, ExcNotImplemented());
231 
232  const unsigned int n0 = q.size();
233  const unsigned int n1 = (dim > 1) ? n0 : 1;
234  const unsigned int n2 = (dim > 2) ? n0 : 1;
235 
236  unsigned int k = 0;
237  for (unsigned int i2 = 0; i2 < n2; ++i2)
238  for (unsigned int i1 = 0; i1 < n1; ++i1)
239  for (unsigned int i0 = 0; i0 < n0; ++i0)
240  {
241  quadrature_points[k](0) = q.point(i0)(0);
242  if (dim > 1)
243  quadrature_points[k](1) = q.point(i1)(0);
244  if (dim > 2)
245  quadrature_points[k](2) = q.point(i2)(0);
246  weights[k] = q.weight(i0);
247  if (dim > 1)
248  weights[k] *= q.weight(i1);
249  if (dim > 2)
250  weights[k] *= q.weight(i2);
251  ++k;
252  }
253 
254  tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
255  for (unsigned int i = 0; i < dim; ++i)
256  (*tensor_basis)[i] = q;
257 }
258 
259 
260 
261 template <int dim>
263  : Subscriptor()
265  , weights(q.weights)
267 {
268  if (dim > 1 && is_tensor_product_flag)
269  tensor_basis =
270  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
271 }
272 
273 
274 
275 template <int dim>
278 {
279  weights = q.weights;
282  if (dim > 1 && is_tensor_product_flag)
283  {
284  if (tensor_basis == nullptr)
285  tensor_basis =
286  std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
287  else
289  }
290  return *this;
291 }
292 
293 
294 
295 template <int dim>
296 bool
298 {
299  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
300 }
301 
302 
303 
304 template <int dim>
305 std::size_t
307 {
310 }
311 
312 
313 
314 template <int dim>
315 typename std::conditional<dim == 1,
316  std::array<Quadrature<1>, dim>,
317  const std::array<Quadrature<1>, dim> &>::type
319 {
320  Assert(this->is_tensor_product_flag == true,
321  ExcMessage("This function only makes sense if "
322  "this object represents a tensor product!"));
323  Assert(tensor_basis != nullptr, ExcInternalError());
324 
325  return *tensor_basis;
326 }
327 
328 
329 #ifndef DOXYGEN
330 template <>
331 std::array<Quadrature<1>, 1>
333 {
334  Assert(this->is_tensor_product_flag == true,
335  ExcMessage("This function only makes sense if "
336  "this object represents a tensor product!"));
337 
338  return std::array<Quadrature<1>, 1>{{*this}};
339 }
340 #endif
341 
342 
343 
344 //---------------------------------------------------------------------------
345 template <int dim>
347  : Quadrature<dim>(qx.size())
348 {
349  Assert(dim == 1, ExcImpossibleInDim(dim));
350  unsigned int k = 0;
351  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
352  {
353  this->quadrature_points[k](0) = qx.point(k1)(0);
354  this->weights[k++] = qx.weight(k1);
355  }
356  Assert(k == this->size(), ExcInternalError());
357  this->is_tensor_product_flag = true;
358 }
359 
360 
361 
362 template <int dim>
364  const Quadrature<1> &qy)
365  : Quadrature<dim>(qx.size() * qy.size())
366 {
367  Assert(dim == 2, ExcImpossibleInDim(dim));
368 }
369 
370 
371 
372 template <>
374  : Quadrature<2>(qx.size() * qy.size())
375 {
376  unsigned int k = 0;
377  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
378  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
379  {
380  this->quadrature_points[k](0) = qx.point(k1)(0);
381  this->quadrature_points[k](1) = qy.point(k2)(0);
382  this->weights[k++] = qx.weight(k1) * qy.weight(k2);
383  }
384  Assert(k == this->size(), ExcInternalError());
385  this->is_tensor_product_flag = true;
386  const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
387  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 2>>(q_array);
388 }
389 
390 
391 
392 template <int dim>
394  const Quadrature<1> &qy,
395  const Quadrature<1> &qz)
396  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
397 {
398  Assert(dim == 3, ExcImpossibleInDim(dim));
399 }
400 
401 
402 
403 template <>
405  const Quadrature<1> &qy,
406  const Quadrature<1> &qz)
407  : Quadrature<3>(qx.size() * qy.size() * qz.size())
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](1) = qy.point(k2)(0);
416  this->quadrature_points[k](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  this->is_tensor_product_flag = true;
421  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
422  this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 3>>(q_array);
423 }
424 
425 
426 
427 // ------------------------------------------------------------ //
428 
429 namespace internal
430 {
431  namespace QIteratedImplementation
432  {
433  namespace
434  {
435  bool
436  uses_both_endpoints(const Quadrature<1> &base_quadrature)
437  {
438  const bool at_left =
439  std::any_of(base_quadrature.get_points().cbegin(),
440  base_quadrature.get_points().cend(),
441  [](const Point<1> &p) { return p == Point<1>{0.}; });
442  const bool at_right =
443  std::any_of(base_quadrature.get_points().cbegin(),
444  base_quadrature.get_points().cend(),
445  [](const Point<1> &p) { return p == Point<1>{1.}; });
446  return (at_left && at_right);
447  }
448  } // namespace
449  } // namespace QIteratedImplementation
450 } // namespace internal
451 
452 // template <>
453 // void
454 // QIterated<1>::fill(Quadrature<1>& dst,
455 // const Quadrature<1> &base_quadrature,
456 // const unsigned int n_copies)
457 // {
458 // Assert (n_copies > 0, ExcZero());
459 // Assert (base_quadrature.size() > 0, ExcZero());
460 
461 // const unsigned int np =
462 // uses_both_endpoints(base_quadrature)
463 // ? (base_quadrature.size()-1) * n_copies + 1
464 // : base_quadrature.size() * n_copies;
465 
466 // dst.quadrature_points.resize(np);
467 // dst.weights.resize(np);
468 
469 // if (!uses_both_endpoints(base_quadrature))
470 // // we don't have to skip some
471 // // points in order to get a
472 // // reasonable quadrature formula
473 // {
474 // unsigned int next_point = 0;
475 // for (unsigned int copy=0; copy<n_copies; ++copy)
476 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
477 // {
478 // dst.quadrature_points[next_point](0)
479 // = (copy + base_quadrature.point(q_point)(0)) / n_copies;
480 // dst.weights[next_point]
481 // = base_quadrature.weight(q_point) / n_copies;
482 // ++next_point;
483 // }
484 // }
485 // else
486 // // skip doubly available points
487 // {
488 // unsigned int next_point = 0;
489 
490 // // first find out the weights of
491 // // the left and the right boundary
492 // // points. note that these usually
493 // // are but need not necessarily be
494 // // the same
495 // double double_point_weight = 0;
496 // unsigned int n_end_points = 0;
497 // for (unsigned int i=0; i<base_quadrature.size(); ++i)
498 // // add up the weight if this
499 // // is an endpoint
500 // if ((base_quadrature.point(i)(0) == 0.) ||
501 // (base_quadrature.point(i)(0) == 1.))
502 // {
503 // double_point_weight += base_quadrature.weight(i);
504 // ++n_end_points;
505 // }
506 // // scale the weight correctly
507 // double_point_weight /= n_copies;
508 
509 // // make sure the base quadrature formula
510 // // has only one quadrature point
511 // // per end point
512 // Assert (n_end_points == 2, ExcInvalidQuadratureFormula());
513 
514 
515 // for (unsigned int copy=0; copy<n_copies; ++copy)
516 // for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
517 // {
518 // // skip the left point of
519 // // this copy since we
520 // // have already entered
521 // // it the last time
522 // if ((copy > 0) &&
523 // (base_quadrature.point(q_point)(0) == 0.))
524 // continue;
525 
526 // dst.quadrature_points[next_point](0)
527 // = (copy+base_quadrature.point(q_point)(0)) / n_copies;
528 
529 // // if this is the
530 // // rightmost point of one
531 // // of the non-last
532 // // copies: give it the
533 // // double weight
534 // if ((copy != n_copies-1) &&
535 // (base_quadrature.point(q_point)(0) == 1.))
536 // dst.weights[next_point] = double_point_weight;
537 // else
538 // dst.weights[next_point] = base_quadrature.weight(q_point) /
539 // n_copies;
540 
541 // ++next_point;
542 // }
543 // }
544 
545 // #if DEBUG
546 // double sum_of_weights = 0;
547 // for (unsigned int i=0; i<dst.size(); ++i)
548 // sum_of_weights += dst.weight(i);
549 // Assert (std::fabs(sum_of_weights-1) < 1e-15,
550 // ExcInternalError());
551 // #endif
552 
553 // }
554 
555 
556 template <>
557 QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
558  : Quadrature<0>()
559 {
560  Assert(false, ExcNotImplemented());
561 }
562 
563 
564 
565 template <>
566 QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
567  const unsigned int n_copies)
568  : Quadrature<1>(
569  internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
570  (base_quadrature.size() - 1) * n_copies + 1 :
571  base_quadrature.size() * n_copies)
572 {
573  // fill(*this, base_quadrature, n_copies);
574  Assert(base_quadrature.size() > 0, ExcNotInitialized());
575  Assert(n_copies > 0, ExcZero());
576 
577  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
578  // we don't have to skip some
579  // points in order to get a
580  // reasonable quadrature formula
581  {
582  unsigned int next_point = 0;
583  for (unsigned int copy = 0; copy < n_copies; ++copy)
584  for (unsigned int q_point = 0; q_point < base_quadrature.size();
585  ++q_point)
586  {
587  this->quadrature_points[next_point] =
588  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
589  (1.0 * copy) / n_copies);
590  this->weights[next_point] =
591  base_quadrature.weight(q_point) / n_copies;
592 
593  ++next_point;
594  }
595  }
596  else
597  // skip doubly available points
598  {
599  unsigned int next_point = 0;
600 
601  // first find out the weights of
602  // the left and the right boundary
603  // points. note that these usually
604  // are but need not necessarily be
605  // the same
606  double double_point_weight = 0;
607  unsigned int n_end_points = 0;
608  for (unsigned int i = 0; i < base_quadrature.size(); ++i)
609  // add up the weight if this
610  // is an endpoint
611  if ((base_quadrature.point(i) == Point<1>(0.0)) ||
612  (base_quadrature.point(i) == Point<1>(1.0)))
613  {
614  double_point_weight += base_quadrature.weight(i);
615  ++n_end_points;
616  }
617  // scale the weight correctly
618  double_point_weight /= n_copies;
619 
620  // make sure the base quadrature formula
621  // has only one quadrature point
622  // per end point
623  Assert(n_end_points == 2, ExcInvalidQuadratureFormula());
624 
625 
626  for (unsigned int copy = 0; copy < n_copies; ++copy)
627  for (unsigned int q_point = 0; q_point < base_quadrature.size();
628  ++q_point)
629  {
630  // skip the left point of
631  // this copy since we
632  // have already entered
633  // it the last time
634  if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
635  continue;
636 
637  this->quadrature_points[next_point] =
638  Point<1>(base_quadrature.point(q_point)(0) / n_copies +
639  (1.0 * copy) / n_copies);
640 
641  // if this is the
642  // rightmost point of one
643  // of the non-last
644  // copies: give it the
645  // double weight
646  if ((copy != n_copies - 1) &&
647  (base_quadrature.point(q_point) == Point<1>(1.0)))
648  this->weights[next_point] = double_point_weight;
649  else
650  this->weights[next_point] =
651  base_quadrature.weight(q_point) / n_copies;
652 
653  ++next_point;
654  }
655  }
656 
657 #if DEBUG
658  double sum_of_weights = 0;
659  for (unsigned int i = 0; i < this->size(); ++i)
660  sum_of_weights += this->weight(i);
661  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
662 #endif
663 }
664 
665 
666 // template <int dim>
667 // void
668 // QIterated<dim>::fill(Quadrature<dim>&, const Quadrature<1>&, unsigned int)
669 // {
670 // Assert(false, ExcNotImplemented());
671 // }
672 
673 
674 // construct higher dimensional quadrature formula by tensor product
675 // of lower dimensional iterated quadrature formulae
676 template <int dim>
678  const unsigned int N)
679  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
680  QIterated<1>(base_quadrature, N))
681 {}
682 
683 
684 
685 // explicit instantiations; note: we need them all for all dimensions
686 template class Quadrature<0>;
687 template class Quadrature<1>;
688 template class Quadrature<2>;
689 template class Quadrature<3>;
690 template class QAnisotropic<1>;
691 template class QAnisotropic<2>;
692 template class QAnisotropic<3>;
693 template class QIterated<1>;
694 template class QIterated<2>;
695 template class QIterated<3>;
696 
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
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:318
static ::ExceptionBase & ExcNotInitialized()
Definition: point.h:110
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:277
T fixed_power(const T t)
Definition: utilities.h:1045
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:346
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:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
std::size_t memory_consumption() const
Definition: quadrature.cc:306
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:31
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:297
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:677
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()