35 , is_tensor_product_flag(false)
45 , is_tensor_product_flag(dim == 1)
53 const std::vector<double> &
w)
58 is_tensor_product_flag = dim == 1;
65 const std::vector<double> &weights)
68 , is_tensor_product_flag(dim == 1)
78 std::vector<double> &&weights)
80 , weights(std::move(weights))
81 , is_tensor_product_flag(dim == 1)
92 , weights(points.size(), std::numeric_limits<
double>::infinity())
93 , is_tensor_product_flag(dim == 1)
104 , weights(std::vector<
double>(1, 1.))
105 , is_tensor_product_flag(true)
106 , tensor_basis(new std::array<
Quadrature<1>, dim>())
108 for (
unsigned int i = 0; i < dim; ++i)
110 const std::vector<Point<1>> quad_vec_1d(1,
Point<1>(
point[i]));
121 , weights(std::vector<
double>(1, 1.))
122 , is_tensor_product_flag(true)
129 : is_tensor_product_flag(false)
148 , weights(q1.size() * q2.size())
149 , is_tensor_product_flag(q1.is_tensor_product())
151 unsigned int present_index = 0;
152 for (
unsigned int i2 = 0; i2 < q2.
size(); ++i2)
153 for (
unsigned int i1 = 0; i1 < q1.
size(); ++i1)
157 for (
unsigned int d = 0;
d < dim - 1; ++
d)
170 for (
unsigned int i = 0; i < size(); ++i)
178 if (is_tensor_product_flag)
180 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
181 for (
unsigned int i = 0; i < dim - 1; ++i)
183 (*tensor_basis)[dim - 1] = q2;
193 , is_tensor_product_flag(true)
195 unsigned int present_index = 0;
196 for (
unsigned int i2 = 0; i2 < q2.
size(); ++i2)
211 for (
unsigned int i = 0; i <
size(); ++i)
227 , is_tensor_product_flag(false)
248 , is_tensor_product_flag(true)
252 const unsigned int n0 = q.
size();
253 const unsigned int n1 = (dim > 1) ? n0 : 1;
254 const unsigned int n2 = (dim > 2) ? n0 : 1;
257 for (
unsigned int i2 = 0; i2 < n2; ++i2)
258 for (
unsigned int i1 = 0; i1 < n1; ++i1)
259 for (
unsigned int i0 = 0; i0 < n0; ++i0)
274 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
275 for (
unsigned int i = 0; i < dim; ++i)
286 , is_tensor_product_flag(q.is_tensor_product_flag)
290 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
302 if (dim > 1 && is_tensor_product_flag)
304 if (tensor_basis ==
nullptr)
306 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
335 typename std::conditional<dim == 1,
336 std::array<Quadrature<1>, dim>,
337 const std::array<Quadrature<1>, dim> &>::type
340 Assert(this->is_tensor_product_flag ==
true,
341 ExcMessage(
"This function only makes sense if "
342 "this object represents a tensor product!"));
345 return *tensor_basis;
351 std::array<Quadrature<1>, 1>
355 ExcMessage(
"This function only makes sense if "
356 "this object represents a tensor product!"));
358 return std::array<Quadrature<1>, 1>{{*
this}};
371 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
390 constexpr
int dim_1 = dim == 2 ? 1 : 0;
393 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
394 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
403 this->
tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
414 :
Quadrature<dim>(qx.size() * qy.size() * qz.size())
419 constexpr
int dim_1 = dim == 3 ? 1 : 0;
420 constexpr
int dim_2 = dim == 3 ? 2 : 0;
423 for (
unsigned int k3 = 0; k3 < qz.
size(); ++k3)
424 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
425 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
435 this->
tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
447 namespace QIteratedImplementation
455 std::any_of(base_quadrature.
get_points().cbegin(),
457 [](
const Point<1> &p) { return p == Point<1>{0.}; });
458 const bool at_right =
459 std::any_of(base_quadrature.
get_points().cbegin(),
461 [](
const Point<1> &p) { return p == Point<1>{1.}; });
462 return (at_left && at_right);
465 std::vector<Point<1>>
466 create_equidistant_interval_points(
const unsigned int n_copies)
468 std::vector<Point<1>> support_points(n_copies + 1);
471 support_points[
copy][0] =
472 static_cast<double>(
copy) /
static_cast<double>(n_copies);
474 support_points[n_copies][0] = 1.0;
476 return support_points;
504 const std::vector<
Point<1>> &intervals)
506 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
507 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
508 base_quadrature.size() * (intervals.size() - 1))
513 const unsigned int n_copies = intervals.size() - 1;
515 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
519 unsigned int next_point = 0;
521 for (
unsigned int q_point = 0; q_point < base_quadrature.
size();
526 (intervals[
copy + 1][0] - intervals[
copy][0]) +
529 base_quadrature.
weight(q_point) *
530 (intervals[
copy + 1][0] - intervals[
copy][0]);
538 const unsigned int left_index =
539 std::distance(base_quadrature.
get_points().begin(),
540 std::find_if(base_quadrature.
get_points().cbegin(),
543 return p == Point<1>{0.};
546 const unsigned int right_index =
547 std::distance(base_quadrature.
get_points().begin(),
548 std::find_if(base_quadrature.
get_points().cbegin(),
551 return p == Point<1>{1.};
554 const unsigned double_point_offset =
555 left_index + (base_quadrature.size() - right_index);
557 for (
unsigned int copy = 0, next_point = 0;
copy < n_copies; ++
copy)
558 for (
unsigned int q_point = 0; q_point < base_quadrature.size();
563 if ((
copy > 0) && (base_quadrature.point(q_point) ==
Point<1>(0.0)))
567 base_quadrature.point(q_point)(0) *
568 (intervals[
copy + 1][0] - intervals[
copy][0]) +
569 intervals[
copy][0])) < 1
e-10 ,
572 this->weights[next_point - double_point_offset] +=
573 base_quadrature.weight(q_point) *
574 (intervals[
copy + 1][0] - intervals[
copy][0]);
580 Point<1>(base_quadrature.point(q_point)(0) *
581 (intervals[
copy + 1][0] - intervals[
copy][0]) +
586 this->weights[next_point] =
587 base_quadrature.weight(q_point) *
588 (intervals[
copy + 1][0] - intervals[
copy][0]);
598 if (std::abs(i[0] - 0.0) < 1
e-12)
600 else if (std::abs(i[0] - 1.0) < 1
e-12)
604 double sum_of_weights = 0;
605 for (
unsigned int i = 0; i < this->size(); ++i)
606 sum_of_weights += this->weight(i);
615 const unsigned int n_copies)
618 internal::QIteratedImplementation::create_equidistant_interval_points(
631 const std::vector<
Point<1>> &intervals)
633 QIterated<1>(base_quadrature, intervals))
640 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
std::unique_ptr< std::array< Quadrature< 1 >, dim > > tensor_basis
Quadrature & operator=(const Quadrature< dim > &)
const std::vector< Point< dim > > & get_points() const
std::size_t memory_consumption() const
void initialize(const std::vector< Point< dim >> &points, const std::vector< double > &weights)
bool is_tensor_product_flag
Quadrature(const unsigned int n_quadrature_points=0)
double weight(const unsigned int i) const
const std::array< Quadrature< 1 >, dim > & get_tensor_basis() const
bool operator==(const Quadrature< dim > &p) const
const Point< dim > & point(const unsigned int i) const
std::vector< double > weights
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
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcZero()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression fabs(const Expression &x)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
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 > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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)
T sum(const T &t, const MPI_Comm mpi_communicator)
void copy(const T *begin, const T *end, U *dest)