Reference documentation for deal.II version 9.4.0
\(\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
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
31template <>
32Quadrature<0>::Quadrature(const unsigned int n_q)
34 , weights(n_q, 0)
35 , is_tensor_product_flag(false)
36{}
37#endif
38
39
40
41template <int dim>
42Quadrature<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
50template <int dim>
51void
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
63template <int dim>
64Quadrature<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
76template <int dim>
77Quadrature<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
88template <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
105template <>
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
114template <>
116 : is_tensor_product_flag(false)
117{
118 Assert(false, ExcImpossibleInDim(0));
119}
120
121
122
123template <>
124Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
125{
126 Assert(false, ExcImpossibleInDim(0));
127}
128#endif // DOXYGEN
129
130
131
132template <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
165 if (is_tensor_product_flag)
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
176template <>
177Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2)
178 : quadrature_points(q2.size())
179 , weights(q2.size())
180 , is_tensor_product_flag(true)
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
209template <>
211 : Subscriptor()
213 , weights(1, 1.)
214 , is_tensor_product_flag(false)
215{}
216
217
218template <>
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
230template <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
268template <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)
277 std::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
278}
279
280
281
282template <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
302template <int dim>
303bool
305{
306 return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
307}
308
309
310
311template <int dim>
312std::size_t
314{
317}
318
319
320
321template <int dim>
322typename 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
337template <>
338std::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//---------------------------------------------------------------------------
352template <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
369template <int dim>
371 const Quadrature<1> &qy)
372 : Quadrature<dim>(qx.size() * qy.size())
373{
374 Assert(dim == 2, ExcImpossibleInDim(dim));
375}
376
377
378
379template <>
381 : Quadrature<2>(qx.size() * qy.size())
382{
383 unsigned int k = 0;
384 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
385 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
386 {
387 this->quadrature_points[k](0) = qx.point(k1)(0);
388 this->quadrature_points[k](1) = qy.point(k2)(0);
389 this->weights[k++] = qx.weight(k1) * qy.weight(k2);
390 }
391 Assert(k == this->size(), ExcInternalError());
392 this->is_tensor_product_flag = true;
393 const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
394 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 2>>(q_array);
395}
396
397
398
399template <int dim>
401 const Quadrature<1> &qy,
402 const Quadrature<1> &qz)
403 : Quadrature<dim>(qx.size() * qy.size() * qz.size())
404{
405 Assert(dim == 3, ExcImpossibleInDim(dim));
406}
407
408
409
410template <>
412 const Quadrature<1> &qy,
413 const Quadrature<1> &qz)
414 : Quadrature<3>(qx.size() * qy.size() * qz.size())
415{
416 unsigned int k = 0;
417 for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
418 for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
419 for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
420 {
421 this->quadrature_points[k](0) = qx.point(k1)(0);
422 this->quadrature_points[k](1) = qy.point(k2)(0);
423 this->quadrature_points[k](2) = qz.point(k3)(0);
424 this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
425 }
426 Assert(k == this->size(), ExcInternalError());
427 this->is_tensor_product_flag = true;
428 const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
429 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, 3>>(q_array);
430}
431
432
433
434// ------------------------------------------------------------ //
435
436namespace internal
437{
438 namespace QIteratedImplementation
439 {
440 namespace
441 {
442 bool
443 uses_both_endpoints(const Quadrature<1> &base_quadrature)
444 {
445 const bool at_left =
446 std::any_of(base_quadrature.get_points().cbegin(),
447 base_quadrature.get_points().cend(),
448 [](const Point<1> &p) { return p == Point<1>{0.}; });
449 const bool at_right =
450 std::any_of(base_quadrature.get_points().cbegin(),
451 base_quadrature.get_points().cend(),
452 [](const Point<1> &p) { return p == Point<1>{1.}; });
453 return (at_left && at_right);
454 }
455
456 std::vector<Point<1>>
457 create_equidistant_interval_points(const unsigned int n_copies)
458 {
459 std::vector<Point<1>> support_points(n_copies + 1);
460
461 for (unsigned int copy = 0; copy < n_copies; ++copy)
462 support_points[copy][0] =
463 static_cast<double>(copy) / static_cast<double>(n_copies);
464
465 support_points[n_copies][0] = 1.0;
466
467 return support_points;
468 }
469 } // namespace
470 } // namespace QIteratedImplementation
471} // namespace internal
472
473
474
475template <>
476QIterated<0>::QIterated(const Quadrature<1> &, const std::vector<Point<1>> &)
477 : Quadrature<0>()
478{
479 Assert(false, ExcNotImplemented());
480}
481
482
483
484template <>
485QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
486 : Quadrature<0>()
487{
488 Assert(false, ExcNotImplemented());
489}
490
491
492
493template <>
495 const std::vector<Point<1>> &intervals)
496 : Quadrature<1>(
497 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
498 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
499 base_quadrature.size() * (intervals.size() - 1))
500{
501 Assert(base_quadrature.size() > 0, ExcNotInitialized());
502 Assert(intervals.size() > 1, ExcZero());
503
504 const unsigned int n_copies = intervals.size() - 1;
505
506 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
507 // we don't have to skip some points in order to get a reasonable quadrature
508 // formula
509 {
510 unsigned int next_point = 0;
511 for (unsigned int copy = 0; copy < n_copies; ++copy)
512 for (unsigned int q_point = 0; q_point < base_quadrature.size();
513 ++q_point)
514 {
515 this->quadrature_points[next_point] =
516 Point<1>(base_quadrature.point(q_point)(0) *
517 (intervals[copy + 1][0] - intervals[copy][0]) +
518 intervals[copy][0]);
519 this->weights[next_point] =
520 base_quadrature.weight(q_point) *
521 (intervals[copy + 1][0] - intervals[copy][0]);
522
523 ++next_point;
524 }
525 }
526 else
527 // skip doubly available points
528 {
529 const unsigned int left_index =
530 std::distance(base_quadrature.get_points().begin(),
531 std::find_if(base_quadrature.get_points().cbegin(),
532 base_quadrature.get_points().cend(),
533 [](const Point<1> &p) {
534 return p == Point<1>{0.};
535 }));
536
537 const unsigned int right_index =
538 std::distance(base_quadrature.get_points().begin(),
539 std::find_if(base_quadrature.get_points().cbegin(),
540 base_quadrature.get_points().cend(),
541 [](const Point<1> &p) {
542 return p == Point<1>{1.};
543 }));
544
545 const unsigned double_point_offset =
546 left_index + (base_quadrature.size() - right_index);
547
548 for (unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
549 for (unsigned int q_point = 0; q_point < base_quadrature.size();
550 ++q_point)
551 {
552 // skip the left point of this copy since we have already entered it
553 // the last time
554 if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
555 {
556 Assert(this->quadrature_points[next_point - double_point_offset]
557 .distance(Point<1>(
558 base_quadrature.point(q_point)(0) *
559 (intervals[copy + 1][0] - intervals[copy][0]) +
560 intervals[copy][0])) < 1e-10 /*tolerance*/,
562
563 this->weights[next_point - double_point_offset] +=
564 base_quadrature.weight(q_point) *
565 (intervals[copy + 1][0] - intervals[copy][0]);
566
567 continue;
568 }
569
570 this->quadrature_points[next_point] =
571 Point<1>(base_quadrature.point(q_point)(0) *
572 (intervals[copy + 1][0] - intervals[copy][0]) +
573 intervals[copy][0]);
574
575 // if this is the rightmost point of one of the non-last copies:
576 // give it the double weight
577 this->weights[next_point] =
578 base_quadrature.weight(q_point) *
579 (intervals[copy + 1][0] - intervals[copy][0]);
580
581 ++next_point;
582 }
583 }
584
585 // make sure that there is no rounding error for 0.0 and 1.0, since there
586 // are multiple asserts in the library checking for equality without
587 // tolorances
588 for (auto &i : this->quadrature_points)
589 if (std::abs(i[0] - 0.0) < 1e-12)
590 i[0] = 0.0;
591 else if (std::abs(i[0] - 1.0) < 1e-12)
592 i[0] = 1.0;
593
594#if DEBUG
595 double sum_of_weights = 0;
596 for (unsigned int i = 0; i < this->size(); ++i)
597 sum_of_weights += this->weight(i);
598 Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
599#endif
600}
601
602
603
604template <>
606 const unsigned int n_copies)
607 : QIterated<1>(
608 base_quadrature,
609 internal::QIteratedImplementation::create_equidistant_interval_points(
610 n_copies))
611{
612 Assert(base_quadrature.size() > 0, ExcNotInitialized());
613 Assert(n_copies > 0, ExcZero());
614}
615
616
617
618// construct higher dimensional quadrature formula by tensor product
619// of lower dimensional iterated quadrature formulae
620template <int dim>
622 const std::vector<Point<1>> &intervals)
623 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, intervals),
624 QIterated<1>(base_quadrature, intervals))
625{}
626
627
628
629template <int dim>
631 const unsigned int n_copies)
632 : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, n_copies),
633 QIterated<1>(base_quadrature, n_copies))
634{}
635
636
637
638// explicit instantiations; note: we need them all for all dimensions
639template class Quadrature<0>;
640template class Quadrature<1>;
641template class Quadrature<2>;
642template class Quadrature<3>;
643template class QAnisotropic<1>;
644template class QAnisotropic<2>;
645template class QAnisotropic<3>;
646template class QIterated<1>;
647template class QIterated<2>;
648template class QIterated<3>;
649
Definition: point.h:111
QAnisotropic(const Quadrature< 1 > &qx)
Definition: quadrature.cc:353
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
Definition: quadrature.cc:630
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:290
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Definition: quadrature.h:311
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:284
std::size_t memory_consumption() const
Definition: quadrature.cc:313
const Point< dim > & point(const unsigned int i) const
bool is_tensor_product_flag
Definition: quadrature.h:305
Quadrature(const unsigned int n_quadrature_points=0)
Definition: quadrature.cc:42
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
Definition: quadrature.cc:325
double weight(const unsigned int i) const
bool operator==(const Quadrature< dim > &p) const
Definition: quadrature.cc:304
void initialize(const std::vector< Point< dim > > &points, const std::vector< double > &weights)
Definition: quadrature.cc:52
std::vector< double > weights
Definition: quadrature.h:296
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression fabs(const Expression &x)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
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
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)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T fixed_power(const T t)
Definition: utilities.h:1123
void copy(const T *begin, const T *end, U *dest)
STL namespace.
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)