34 , is_tensor_product_flag(false)
42 : quadrature_points(n_q,
Point<dim>())
44 , is_tensor_product_flag(dim == 1)
54 this->weights.clear();
55 if (weights.
size() > 0)
58 this->weights.insert(this->weights.
end(), weights.
begin(), weights.
end());
61 this->weights.resize(points.size(),
62 std::numeric_limits<double>::infinity());
64 quadrature_points.clear();
65 quadrature_points.insert(quadrature_points.end(),
69 is_tensor_product_flag = dim == 1;
76 const std::vector<double> &weights)
77 : quadrature_points(points)
79 , is_tensor_product_flag(dim == 1)
89 std::vector<double> &&weights)
90 : quadrature_points(
std::move(points))
91 , weights(
std::move(weights))
92 , is_tensor_product_flag(dim == 1)
102 : quadrature_points(points)
103 , weights(points.size(),
std::numeric_limits<double>::infinity())
104 , is_tensor_product_flag(dim == 1)
114 : quadrature_points(
std::vector<
Point<dim>>(1, point))
115 , weights(
std::vector<double>(1, 1.))
116 , is_tensor_product_flag(true)
119 for (
unsigned int i = 0; i < dim; ++i)
121 const std::vector<Point<1>> quad_vec_1d(1,
Point<1>(
point[i]));
131 : quadrature_points(
std::vector<
Point<1>>(1, point))
132 , weights(
std::vector<double>(1, 1.))
133 , is_tensor_product_flag(true)
140 : is_tensor_product_flag(false)
158 : quadrature_points(q1.size() * q2.size())
159 , weights(q1.size() * q2.size())
160 , is_tensor_product_flag(q1.is_tensor_product())
162 unsigned int present_index = 0;
163 for (
unsigned int i2 = 0; i2 < q2.
size(); ++i2)
164 for (
unsigned int i1 = 0; i1 < q1.
size(); ++i1)
168 for (
unsigned int d = 0; d < dim - 1; ++d)
181 for (
unsigned int i = 0; i <
size(); ++i)
191 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
192 for (
unsigned int i = 0; i < dim - 1; ++i)
194 (*tensor_basis)[dim - 1] = q2;
202 : quadrature_points(q2.size())
204 , is_tensor_product_flag(true)
206 unsigned int present_index = 0;
207 for (
unsigned int i2 = 0; i2 < q2.
size(); ++i2)
222 for (
unsigned int i = 0; i <
size(); ++i)
238 , is_tensor_product_flag(false)
257 , quadrature_points(
Utilities::fixed_power<dim>(q.size()))
258 , weights(
Utilities::fixed_power<dim>(q.size()))
259 , is_tensor_product_flag(true)
263 const unsigned int n0 = q.
size();
264 const unsigned int n1 = (dim > 1) ? n0 : 1;
265 const unsigned int n2 = (dim > 2) ? n0 : 1;
268 for (
unsigned int i2 = 0; i2 < n2; ++i2)
269 for (
unsigned int i1 = 0; i1 < n1; ++i1)
270 for (
unsigned int i0 = 0; i0 < n0; ++i0)
285 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
286 for (
unsigned int i = 0; i < dim; ++i)
295 , quadrature_points(q.quadrature_points)
297 , is_tensor_product_flag(q.is_tensor_product_flag)
301 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
313 if (dim > 1 && is_tensor_product_flag)
315 if (tensor_basis ==
nullptr)
317 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
346typename std::conditional_t<dim == 1,
347 std::array<Quadrature<1>, dim>,
348 const std::array<Quadrature<1>, dim> &>
351 Assert(this->is_tensor_product_flag ==
true,
352 ExcMessage(
"This function only makes sense if "
353 "this object represents a tensor product!"));
356 return *tensor_basis;
362std::array<Quadrature<1>, 1>
365 Assert(this->is_tensor_product_flag ==
true,
366 ExcMessage(
"This function only makes sense if "
367 "this object represents a tensor product!"));
369 return std::array<Quadrature<1>, 1>{{*
this}};
382 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
384 this->quadrature_points[k][0] = qx.
point(k1)[0];
385 this->weights[k++] = qx.
weight(k1);
388 this->is_tensor_product_flag =
true;
401 constexpr int dim_1 = dim == 2 ? 1 : 0;
404 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
405 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
407 this->quadrature_points[k][0] = qx.
point(k1)[0];
408 this->quadrature_points[k][dim_1] = qy.
point(k2)[0];
413 this->is_tensor_product_flag =
true;
414 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
415 (*this->tensor_basis)[0] = qx;
416 (*this->tensor_basis)[dim_1] = qy;
425 :
Quadrature<dim>(qx.size() * qy.size() * qz.size())
430 constexpr int dim_1 = dim == 3 ? 1 : 0;
431 constexpr int dim_2 = dim == 3 ? 2 : 0;
434 for (
unsigned int k3 = 0; k3 < qz.
size(); ++k3)
435 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
436 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
438 this->quadrature_points[k][0] = qx.
point(k1)[0];
439 this->quadrature_points[k][dim_1] = qy.
point(k2)[0];
440 this->quadrature_points[k][dim_2] = qz.
point(k3)[0];
445 this->is_tensor_product_flag =
true;
446 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
447 (*this->tensor_basis)[0] = qx;
448 (*this->tensor_basis)[dim_1] = qy;
449 (*this->tensor_basis)[dim_2] = qz;
458 namespace QIteratedImplementation
466 std::any_of(base_quadrature.
get_points().cbegin(),
468 [](
const Point<1> &p) { return p == Point<1>{0.}; });
469 const bool at_right =
470 std::any_of(base_quadrature.
get_points().cbegin(),
472 [](
const Point<1> &p) { return p == Point<1>{1.}; });
473 return (at_left && at_right);
476 std::vector<Point<1>>
477 create_equidistant_interval_points(
const unsigned int n_copies)
479 std::vector<Point<1>> support_points(n_copies + 1);
481 for (
unsigned int copy = 0; copy < n_copies; ++copy)
482 support_points[copy][0] =
483 static_cast<double>(copy) /
static_cast<double>(n_copies);
485 support_points[n_copies][0] = 1.0;
487 return support_points;
515 const std::vector<
Point<1>> &intervals)
517 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
518 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
519 base_quadrature.size() * (intervals.size() - 1))
524 const unsigned int n_copies = intervals.size() - 1;
526 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
530 unsigned int next_point = 0;
531 for (
unsigned int copy = 0; copy < n_copies; ++copy)
532 for (
unsigned int q_point = 0; q_point < base_quadrature.
size();
535 this->quadrature_points[next_point] =
537 (intervals[copy + 1][0] - intervals[copy][0]) +
539 this->weights[next_point] =
540 base_quadrature.
weight(q_point) *
541 (intervals[copy + 1][0] - intervals[copy][0]);
549 const unsigned int left_index =
550 std::distance(base_quadrature.
get_points().begin(),
551 std::find_if(base_quadrature.
get_points().cbegin(),
554 return p == Point<1>{0.};
557 const unsigned int right_index =
558 std::distance(base_quadrature.
get_points().begin(),
559 std::find_if(base_quadrature.
get_points().cbegin(),
562 return p == Point<1>{1.};
565 const unsigned double_point_offset =
566 left_index + (base_quadrature.size() - right_index);
568 for (
unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
569 for (
unsigned int q_point = 0; q_point < base_quadrature.size();
574 if ((copy > 0) && (base_quadrature.point(q_point) ==
Point<1>(0.0)))
576 Assert(this->quadrature_points[next_point - double_point_offset]
578 base_quadrature.point(q_point)[0] *
579 (intervals[copy + 1][0] - intervals[copy][0]) +
580 intervals[copy][0])) < 1e-10 ,
583 this->weights[next_point - double_point_offset] +=
584 base_quadrature.weight(q_point) *
585 (intervals[copy + 1][0] - intervals[copy][0]);
591 Point<1>(base_quadrature.point(q_point)[0] *
592 (intervals[copy + 1][0] - intervals[copy][0]) +
597 this->weights[next_point] =
598 base_quadrature.weight(q_point) *
599 (intervals[
copy + 1][0] - intervals[
copy][0]);
611 else if (
std::abs(i[0] - 1.0) < 1e-12)
615 double sum_of_weights = 0;
616 for (
unsigned int i = 0; i < this->size(); ++i)
617 sum_of_weights += this->weight(i);
626 const unsigned int n_copies)
629 internal::QIteratedImplementation::create_equidistant_interval_points(
642 const std::vector<
Point<1>> &intervals)
644 QIterated<1>(base_quadrature, intervals))
651 const unsigned int n_copies)
QAnisotropic(const Quadrature< 1 > &qx)
QIterated(const Quadrature< 1 > &base_quadrature, const unsigned int n_copies)
std::vector< Point< dim > > quadrature_points
void initialize(const ArrayView< const Point< dim > > &points, const ArrayView< const double > &weights={})
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Quadrature & operator=(const Quadrature< dim > &)
std::size_t memory_consumption() const
const Point< dim > & point(const unsigned int i) const
bool is_tensor_product_flag
Quadrature(const unsigned int n_quadrature_points=0)
double weight(const unsigned int i) const
bool operator==(const Quadrature< dim > &p) const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
std::vector< double > weights
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
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={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
void copy(const T *begin, const T *end, U *dest)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)