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)
79 , weights(points.size(), std::numeric_limits<
double>::infinity())
80 , is_tensor_product_flag(dim == 1)
91 , weights(std::vector<
double>(1, 1.))
92 , is_tensor_product_flag(true)
93 , tensor_basis(new std::array<
Quadrature<1>, dim>())
95 for (
unsigned int i = 0; i < dim; ++i)
97 const std::vector<Point<1>> quad_vec_1d(1,
Point<1>(
point[i]));
108 , weights(std::vector<
double>(1, 1.))
109 , is_tensor_product_flag(true)
116 : is_tensor_product_flag(false)
135 , weights(q1.size() * q2.size())
136 , is_tensor_product_flag(q1.is_tensor_product())
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)
144 for (
unsigned int d = 0;
d < dim - 1; ++
d)
157 for (
unsigned int i = 0; i <
size(); ++i)
167 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
168 for (
unsigned int i = 0; i < dim - 1; ++i)
170 (*tensor_basis)[dim - 1] = q2;
180 , is_tensor_product_flag(true)
182 unsigned int present_index = 0;
183 for (
unsigned int i2 = 0; i2 < q2.
size(); ++i2)
198 for (
unsigned int i = 0; i <
size(); ++i)
214 , is_tensor_product_flag(false)
235 , is_tensor_product_flag(true)
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;
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)
261 tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
262 for (
unsigned int i = 0; i < dim; ++i)
273 , is_tensor_product_flag(q.is_tensor_product_flag)
277 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
289 if (dim > 1 && is_tensor_product_flag)
291 if (tensor_basis ==
nullptr)
293 std::make_unique<std::array<Quadrature<1>, dim>>(*q.
tensor_basis);
322 typename std::conditional<dim == 1,
323 std::array<Quadrature<1>, dim>,
324 const std::array<Quadrature<1>, dim> &>::type
327 Assert(this->is_tensor_product_flag ==
true,
328 ExcMessage(
"This function only makes sense if "
329 "this object represents a tensor product!"));
332 return *tensor_basis;
338 std::array<Quadrature<1>, 1>
342 ExcMessage(
"This function only makes sense if "
343 "this object represents a tensor product!"));
345 return std::array<Quadrature<1>, 1>{{*
this}};
358 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
377 constexpr
int dim_1 = dim == 2 ? 1 : 0;
380 for (
unsigned int k2 = 0; k2 < qy.
size(); ++k2)
381 for (
unsigned int k1 = 0; k1 < qx.
size(); ++k1)
390 this->
tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
401 :
Quadrature<dim>(qx.size() * qy.size() * qz.size())
406 constexpr
int dim_1 = dim == 3 ? 1 : 0;
407 constexpr
int dim_2 = dim == 3 ? 2 : 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)
422 this->
tensor_basis = std::make_unique<std::array<Quadrature<1>, dim>>();
434 namespace QIteratedImplementation
442 std::any_of(base_quadrature.
get_points().cbegin(),
444 [](
const Point<1> &p) { return p == Point<1>{0.}; });
445 const bool at_right =
446 std::any_of(base_quadrature.
get_points().cbegin(),
448 [](
const Point<1> &p) { return p == Point<1>{1.}; });
449 return (at_left && at_right);
452 std::vector<Point<1>>
453 create_equidistant_interval_points(
const unsigned int n_copies)
455 std::vector<Point<1>> support_points(n_copies + 1);
458 support_points[
copy][0] =
459 static_cast<double>(
copy) /
static_cast<double>(n_copies);
461 support_points[n_copies][0] = 1.0;
463 return support_points;
491 const std::vector<
Point<1>> &intervals)
493 internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
494 (base_quadrature.size() - 1) * (intervals.size() - 1) + 1 :
495 base_quadrature.size() * (intervals.size() - 1))
500 const unsigned int n_copies = intervals.size() - 1;
502 if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
506 unsigned int next_point = 0;
508 for (
unsigned int q_point = 0; q_point < base_quadrature.
size();
513 (intervals[
copy + 1][0] - intervals[
copy][0]) +
516 base_quadrature.
weight(q_point) *
517 (intervals[
copy + 1][0] - intervals[
copy][0]);
525 const unsigned int left_index =
526 std::distance(base_quadrature.
get_points().begin(),
527 std::find_if(base_quadrature.
get_points().cbegin(),
530 return p == Point<1>{0.};
533 const unsigned int right_index =
534 std::distance(base_quadrature.
get_points().begin(),
535 std::find_if(base_quadrature.
get_points().cbegin(),
538 return p == Point<1>{1.};
541 const unsigned double_point_offset =
542 left_index + (base_quadrature.size() - right_index);
544 for (
unsigned int copy = 0, next_point = 0;
copy < n_copies; ++
copy)
545 for (
unsigned int q_point = 0; q_point < base_quadrature.size();
550 if ((
copy > 0) && (base_quadrature.point(q_point) ==
Point<1>(0.0)))
554 base_quadrature.point(q_point)(0) *
555 (intervals[
copy + 1][0] - intervals[
copy][0]) +
556 intervals[
copy][0])) < 1
e-10 ,
559 this->weights[next_point - double_point_offset] +=
560 base_quadrature.weight(q_point) *
561 (intervals[
copy + 1][0] - intervals[
copy][0]);
567 Point<1>(base_quadrature.point(q_point)(0) *
568 (intervals[
copy + 1][0] - intervals[
copy][0]) +
573 this->weights[next_point] =
574 base_quadrature.weight(q_point) *
575 (intervals[
copy + 1][0] - intervals[
copy][0]);
585 if (std::abs(i[0] - 0.0) < 1
e-12)
587 else if (std::abs(i[0] - 1.0) < 1
e-12)
591 double sum_of_weights = 0;
592 for (
unsigned int i = 0; i < this->size(); ++i)
593 sum_of_weights += this->weight(i);
602 const unsigned int n_copies)
605 internal::QIteratedImplementation::create_equidistant_interval_points(
618 const std::vector<
Point<1>> &intervals)
620 QIterated<1>(base_quadrature, intervals))
627 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
#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< T >::value, 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)