413 constexpr int dim_1 = dim == 2 ? 1 : 0;
416 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
417 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
419 this->quadrature_points[k][0] = qx.
point(k1)[0];
420 this->quadrature_points[k][dim_1] = qy.
point(k2)[0];
425 this->is_tensor_product_flag =
true;
426 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
427 (*this->tensor_basis)[0] = qx;
428 (*this->tensor_basis)[dim_1] = qy;
437 :
Quadrature<dim>(qx.size() * qy.size() * qz.size())
442 constexpr int dim_1 = dim == 3 ? 1 : 0;
443 constexpr int dim_2 = dim == 3 ? 2 : 0;
446 for (
unsigned int k3 = 0; k3 < qz.
size(); ++k3)
447 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
448 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
450 this->quadrature_points[k][0] = qx.
point(k1)[0];
451 this->quadrature_points[k][dim_1] = qy.
point(k2)[0];
452 this->quadrature_points[k][dim_2] = qz.
point(k3)[0];
457 this->is_tensor_product_flag =
true;
458 this->tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
459 (*this->tensor_basis)[0] = qx;
460 (*this->tensor_basis)[dim_1] = qy;
461 (*this->tensor_basis)[dim_2] = qz;
478 std::any_of(base_quadrature.
get_points().cbegin(),
480 [](
const Point<1> &p) { return p == Point<1>{0.}; });
481 const bool at_right =
482 std::any_of(base_quadrature.
get_points().cbegin(),
484 [](
const Point<1> &p) { return p == Point<1>{1.}; });
485 return (at_left && at_right);
488 std::vector<Point<1>>
489 create_equidistant_interval_points(
const unsigned int n_copies)
491 std::vector<Point<1>> support_points(n_copies + 1);
493 for (
unsigned int copy = 0; copy < n_copies; ++copy)
494 support_points[copy][0] =
495 static_cast<double>(copy) /
static_cast<double>(n_copies);
497 support_points[n_copies][0] = 1.0;
499 return support_points;
527 const std::vector<
Point<1>> &intervals)
529 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
530 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
531 base_quadrature.size() * (intervals.size() - 1))
536 const unsigned int n_copies = intervals.size() - 1;
538 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
542 unsigned int next_point = 0;
543 for (
unsigned int copy = 0; copy < n_copies; ++copy)
544 for (
unsigned int q_point = 0; q_point < base_quadrature.
size();
547 this->quadrature_points[next_point] =
549 (intervals[copy + 1][0] - intervals[copy][0]) +
551 this->weights[next_point] =
552 base_quadrature.
weight(q_point) *
553 (intervals[copy + 1][0] - intervals[copy][0]);
561 const unsigned int left_index =
562 std::distance(base_quadrature.
get_points().begin(),
563 std::find_if(base_quadrature.
get_points().cbegin(),
566 return p == Point<1>{0.};
569 const unsigned int right_index =
570 std::distance(base_quadrature.
get_points().begin(),
571 std::find_if(base_quadrature.
get_points().cbegin(),
574 return p == Point<1>{1.};
577 const unsigned double_point_offset =
578 left_index + (base_quadrature.size() - right_index);
580 for (
unsigned int copy = 0, next_point = 0; copy < n_copies; ++copy)
581 for (
unsigned int q_point = 0; q_point < base_quadrature.size();
586 if ((copy > 0) && (base_quadrature.point(q_point) ==
Point<1>(0.0)))
588 Assert(this->quadrature_points[next_point - double_point_offset]
590 base_quadrature.point(q_point)[0] *
591 (intervals[copy + 1][0] - intervals[copy][0]) +
592 intervals[copy][0])) < 1e-10 ,
595 this->weights[next_point - double_point_offset] +=
596 base_quadrature.weight(q_point) *
597 (intervals[copy + 1][0] - intervals[copy][0]);
603 Point<1>(base_quadrature.point(q_point)[0] *
604 (intervals[copy + 1][0] - intervals[copy][0]) +
609 this->weights[next_point] =
610 base_quadrature.weight(q_point) *
611 (intervals[
copy + 1][0] - intervals[
copy][0]);
620 for (
auto &i : this->quadrature_points)
623 else if (
std::abs(i[0] - 1.0) < 1e-12)
627 double sum_of_weights = 0;
628 for (
unsigned int i = 0; i < this->size(); ++i)
629 sum_of_weights += this->weight(i);