38#include <boost/math/special_functions/relative_difference.hpp>
39#include <boost/math/special_functions/sign.hpp>
40#include <boost/math/tools/roots.hpp>
51 namespace QuadratureGeneratorImplementation
60 const unsigned int component_in_dim,
64 ExcMessage(
"Interval start must be less than interval end."));
66 const double length = end - start;
67 for (
unsigned int j = 0; j < quadrature1D.
size(); ++j)
69 const double x = start + length * quadrature1D.
point(j)[0];
71 point, component_in_dim, x),
72 length * weight * quadrature1D.
weight(j));
88 const unsigned int component_in_dim,
91 for (
unsigned int j = 0; j < lower.
size(); ++j)
108 const std::vector<std::reference_wrapper<
const Function<dim>>>
112 Assert(functions.size() > 0,
114 "The incoming vector must contain at least one function."));
116 const int sign_of_first =
117 boost::math::sign(functions[0].get().value(point));
119 if (sign_of_first == 0)
122 for (
unsigned int j = 1; j < functions.size(); ++j)
124 const int sign = boost::math::sign(functions[j].get().value(point));
126 if (sign != sign_of_first)
130 if (sign_of_first < 0)
154 std::pair<double, double> &value_bounds)
156 const ReferenceCell &cube = ReferenceCells::get_hypercube<dim>();
157 for (
unsigned int i = 0; i < cube.
n_vertices(); ++i)
159 const double vertex_value = function.
value(box.
vertex(i));
161 value_bounds.first =
std::min(value_bounds.first, vertex_value);
162 value_bounds.second =
std::max(value_bounds.second, vertex_value);
181 const std::vector<std::reference_wrapper<
const Function<dim>>>
186 all_function_bounds.clear();
187 all_function_bounds.reserve(functions.size());
191 FunctionTools::taylor_estimate_function_bounds<dim>(
195 all_function_bounds.push_back(bounds);
202 std::pair<double, double>
207 std::pair<double, double> extremes = bounds[0].value;
208 for (
unsigned int i = 1; i < bounds.size(); ++i)
210 extremes.first =
std::min(extremes.first, bounds[i].value.first);
211 extremes.second =
std::max(extremes.second, bounds[i].value.second);
226 if (function_bounds.first > 0)
228 if (function_bounds.second < 0)
258 Assert(function_bounds.first <= function_bounds.second,
259 ExcMessage(
"Function bounds reversed, max < min."));
261 return 0.5 * (
std::abs(function_bounds.second + function_bounds.first) -
262 (function_bounds.second - function_bounds.first));
276 std::optional<HeightDirectionData>
282 std::optional<std::array<double, dim>> min_lower_abs_grad;
290 if (!min_lower_abs_grad)
292 min_lower_abs_grad.emplace();
293 for (
unsigned int d = 0; d < dim; ++d)
295 (*min_lower_abs_grad)[d] =
301 for (
unsigned int d = 0; d < dim; ++d)
303 (*min_lower_abs_grad)[d] =
311 if (min_lower_abs_grad)
313 const auto max_element =
314 std::max_element(min_lower_abs_grad->begin(),
315 min_lower_abs_grad->end());
319 std::distance(min_lower_abs_grad->begin(), max_element);
320 data.min_abs_dfdx = *max_element;
325 return std::optional<HeightDirectionData>();
340 if (all_function_bounds.size() != 2)
347 if (bounds0.
value.first > 0 && bounds1.
value.second < 0)
349 if (bounds1.
value.first > 0 && bounds0.
value.second < 0)
370 for (
unsigned int i = 0; i < unit_quadrature.
size(); ++i)
373 const double weight = unit_quadrature.
weight(i) * box.
volume();
394 const std::vector<std::reference_wrapper<
const Function<dim>>>
397 const unsigned int direction,
402 restrictions.clear();
403 restrictions.reserve(2 * functions.size());
408 for (
const auto &function : functions)
411 function, direction, bottom));
413 function, direction, top));
430 const std::vector<std::reference_wrapper<
const Function<dim>>>
433 const unsigned int open_direction,
438 restrictions.clear();
439 restrictions.reserve(functions.size());
440 for (
const auto &function : functions)
443 function, open_direction, point));
467 const std::vector<double> &roots,
470 const unsigned int height_function_direction,
471 const std::vector<std::reference_wrapper<
const Function<1>>>
477 const int n_roots = roots.size();
480 for (
int i = -1; i < n_roots; ++i)
483 const double start = i < 0 ? interval.
lower_bound(0) : roots[i];
485 i + 1 < n_roots ? roots[i + 1] : interval.
upper_bound(0);
487 const double length = end - start;
497 const Point<1> center(start + 0.5 * length);
508 height_function_direction,
517 const double tolerance,
518 const unsigned int max_recursion_depth,
519 const unsigned int max_iterations)
520 : tolerance(tolerance)
521 , max_recursion_depth(max_recursion_depth)
522 , max_iterations(max_iterations)
535 const std::vector<std::reference_wrapper<
const Function<1>>> &functions,
537 std::vector<double> &roots)
541 const unsigned int recursion_depth = 0;
542 find_roots(function, interval, recursion_depth, roots);
545 std::sort(roots.begin(), roots.end());
547 const auto roots_are_equal = [
this](
const double &a,
const double &b) {
550 roots.erase(unique(roots.begin(), roots.end(), roots_are_equal),
559 const unsigned int recursion_depth,
560 std::vector<double> &roots)
563 const double left_value = function.
value(interval.
vertex(0));
564 const double right_value = function.
value(interval.
vertex(1));
567 if (boost::math::sign(left_value) != boost::math::sign(right_value))
569 const auto lambda = [&function](
const double x) {
573 const auto stopping_criteria = [
this](
const double &a,
580 const std::pair<double, double> root_bracket =
581 boost::math::tools::toms748_solve(lambda,
589 const double root = .5 * (root_bracket.first + root_bracket.second);
590 roots.push_back(root);
597 std::pair<double, double> value_bounds;
598 std::array<std::pair<double, double>, 1> gradient_bounds;
599 FunctionTools::taylor_estimate_function_bounds<1>(function,
606 const double function_min =
607 std::min({left_value, right_value, value_bounds.first});
610 if (function_min > 0)
613 const double function_max =
614 std::max({left_value, right_value, value_bounds.second});
617 if (function_max < 0)
651 this->quadrature_points.clear();
652 this->weights.clear();
662 this->quadrature_points.push_back(point);
663 this->weights.push_back(weight);
673 switch (definiteness)
708 const unsigned int direction,
716 const double bottom_value = level_set.
value(bottom_point);
722 const double top_value = level_set.
value(top_point);
732 template <
int dim,
int spacedim>
736 : q_collection1D(&q_collection1D)
737 , additional_data(additional_data)
739 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
740 additional_data.max_root_finder_splits))
747 template <
int dim,
int spacedim>
750 const std::vector<std::reference_wrapper<
const Function<dim>>>
754 const unsigned int height_function_direction,
757 const Quadrature<1> &quadrature1D = (*q_collection1D)[q_index];
759 for (
unsigned int q = 0; q < low_dim_quadrature.
size(); ++q)
761 const Point<dim - 1> &point = low_dim_quadrature.
point(q);
762 const double weight = low_dim_quadrature.
weight(q);
765 height_function_direction,
769 const std::vector<std::reference_wrapper<const Function<1>>>
770 restrictions(point_restrictions.begin(),
771 point_restrictions.end());
774 box.
bounds(height_function_direction);
777 root_finder.find_roots(restrictions, bounds_in_direction, roots);
784 height_function_direction,
790 create_surface_point(point,
794 height_function_direction,
798 point_restrictions.clear();
803 template <
int dim,
int spacedim>
808 const std::vector<std::reference_wrapper<
const Function<dim>>>
811 const unsigned int height_function_direction,
821 if (roots.size() == 1)
824 point, height_function_direction, roots[0]);
840 point, height_function_direction, box, level_set);
845 normal *= 1. / normal.
norm();
849 const double surface_weight =
850 weight * gradient.norm() /
851 std::abs(gradient[height_function_direction]);
852 surface_quadrature.
push_back(surface_point, surface_weight, normal);
857 template <
int dim,
int spacedim>
860 unsigned int q_index)
863 this->q_index = q_index;
868 template <
int dim,
int spacedim>
872 : additional_data(additional_data)
873 , q_collection1D(&q_collection1D)
880 template <
int dim,
int spacedim>
885 , low_dim_algorithm(q_collection1D, additional_data)
886 , up_through_dimension_creator(q_collection1D, additional_data)
894 template <
int dim,
int spacedim>
898 q_partitioning.clear();
903 template <
int dim,
int spacedim>
907 return q_partitioning;
912 template <
int dim,
int spacedim>
915 const std::vector<std::reference_wrapper<
const Function<dim>>>
918 const unsigned int n_box_splits)
920 std::vector<FunctionBounds<dim>> all_function_bounds;
923 const std::pair<double, double> extreme_values =
926 if (extreme_values.first > this->additional_data.limit_to_be_definite)
930 this->q_partitioning.positive);
932 else if (extreme_values.second <
933 -(this->additional_data.limit_to_be_definite))
937 this->q_partitioning.negative);
943 this->q_partitioning.indefinite);
947 const std::optional<HeightDirectionData>
data =
953 this->additional_data.lower_bound_implicit_function)
955 create_low_dim_quadratures(
data->direction,
959 create_high_dim_quadratures(
data->direction, level_sets, box);
961 else if (n_box_splits < this->additional_data.max_box_splits)
963 split_box_and_recurse(level_sets, box,
data, n_box_splits);
969 use_midpoint_method(level_sets, box);
982 std::optional<unsigned int>
986 std::array<std::pair<double, unsigned int>, dim> side_lengths;
987 for (
unsigned int i = 0; i < dim; ++i)
990 side_lengths[i].second = i;
993 std::sort(side_lengths.begin(), side_lengths.end());
997 if (boost::math::epsilon_difference(side_lengths[dim - 1].
first,
998 side_lengths[dim - 2].
first) < 100)
999 return std::optional<unsigned int>();
1001 return side_lengths.back().second;
1021 const std::optional<HeightDirectionData> &height_direction_data)
1023 const std::optional<unsigned int> direction =
1031 if (height_direction_data)
1032 return height_direction_data->direction;
1052 corners.second[direction] -= .5 * box.
side_length(direction);
1071 corners.first[direction] += .5 * box.
side_length(direction);
1078 template <
int dim,
int spacedim>
1081 const std::vector<std::reference_wrapper<
const Function<dim>>>
1084 const std::optional<HeightDirectionData> &direction_data,
1085 const unsigned int n_box_splits)
1087 if (this->additional_data.split_in_half)
1089 const unsigned int direction =
1095 generate(level_sets, left_box, n_box_splits + 1);
1096 generate(level_sets, right_box, n_box_splits + 1);
1100 for (
unsigned int i = 0;
1101 i < GeometryInfo<dim>::max_children_per_cell;
1104 generate(level_sets, box.
child(i), n_box_splits + 1);
1111 template <
int dim,
int spacedim>
1114 const unsigned int height_function_direction,
1115 const std::vector<std::reference_wrapper<
const Function<dim>>>
1118 const unsigned int n_box_splits)
1124 height_function_direction,
1128 const std::vector<std::reference_wrapper<
const Function<dim - 1>>>
1129 restrictions(face_restrictions.begin(), face_restrictions.end());
1134 low_dim_algorithm.clear_quadratures();
1135 low_dim_algorithm.generate(restrictions, cross_section, n_box_splits);
1140 template <
int dim,
int spacedim>
1143 const unsigned int height_function_direction,
1144 const std::vector<std::reference_wrapper<
const Function<dim>>>
1149 low_dim_algorithm.get_quadratures();
1152 (*this->q_collection1D)[this->q_index];
1158 height_function_direction,
1159 this->q_partitioning.negative);
1165 height_function_direction,
1166 this->q_partitioning.positive);
1168 up_through_dimension_creator.generate(level_sets,
1170 low_dim_quadratures.indefinite,
1171 height_function_direction,
1172 this->q_partitioning);
1177 template <
int dim,
int spacedim>
1189 this->q_partitioning.quadrature_by_definiteness(definiteness);
1196 template <
int dim,
int spacedim>
1202 this->q_index = q_index;
1203 low_dim_algorithm.set_1D_quadrature(q_index);
1204 up_through_dimension_creator.set_1D_quadrature(q_index);
1209 template <
int spacedim>
1215 RootFinder::AdditionalData(additional_data.root_finder_tolerance,
1216 additional_data.max_root_finder_splits))
1219 ExcMessage(
"Incoming quadrature collection is empty."));
1224 template <
int spacedim>
1227 const std::vector<std::reference_wrapper<
const Function<1>>>
1230 const unsigned int n_box_splits)
1235 root_finder.find_roots(level_sets, box, roots);
1238 (*this->q_collection1D)[this->q_index];
1247 this->additional_data,
1248 this->q_partitioning);
1251 this->create_surface_points(level_sets);
1256 template <
int spacedim>
1259 const std::vector<std::reference_wrapper<
const Function<1>>>
1264 for (
const double root : roots)
1268 const double weight = 1;
1271 Tensor<1, 1> normal = level_sets[0].get().gradient(point);
1272 const double gradient_norm = normal.
norm();
1274 gradient_norm > 1e-11,
1276 "The level set function has a gradient almost equal to 0."));
1277 normal *= 1. / gradient_norm;
1279 this->q_partitioning.surface.push_back(point, weight, normal);
1285 template <
int spacedim>
1290 this->q_index = q_index;
1299 namespace DiscreteQuadratureGeneratorImplementation
1303 "The set_active_cell function has to be called before calling this function.");
1306 ExcReferenceCellNotHypercube,
1307 "The reference cell of the incoming cell must be a hypercube.");
1335 template <
int dim,
typename Number>
1360 set_subcell(
const std::vector<unsigned int> &mask,
1367 is_fe_q_iso_q1()
const override;
1373 n_subdivisions()
const override;
1384 const unsigned int component = 0)
const override;
1395 const unsigned int component = 0)
const override;
1406 const unsigned int component = 0)
const override;
1413 cell_is_set()
const;
1465 std::vector<Polynomials::Polynomial<double>>
poly;
1485 template <
int dim,
typename Number>
1489 : dof_handler(&dof_handler)
1490 , global_dof_values(&dof_values)
1491 , n_subdivisions_per_line(
numbers::invalid_unsigned_int)
1499 template <
int dim,
typename Number>
1505 &cell->get_triangulation() == &dof_handler->get_triangulation(),
1507 "The incoming cell must belong to the triangulation associated with "
1508 "the DoFHandler passed to the constructor."));
1510 const auto dof_handler_cell =
1511 cell->as_dof_handler_iterator(*dof_handler);
1519 if (element != &dof_handler_cell->get_fe())
1522 element = &dof_handler_cell->get_fe();
1528 &dof_handler_cell->get_fe()))
1530 this->n_subdivisions_per_line = fe_q_iso_q1->get_degree();
1533 .value_or_initialize(
1534 []() {
return std::make_unique<FE_Q<dim>>(1); })
1536 local_dof_values_subcell.resize(
1537 fe_q1.value()->n_dofs_per_cell());
1555 polynomials_are_hat_functions =
1556 (poly.size() == 2 && poly[0].value(0.) == 1. &&
1557 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1558 poly[1].value(1.) == 1.);
1562 element = &dof_handler_cell->get_fe();
1564 local_dof_indices.resize(element->dofs_per_cell);
1565 dof_handler_cell->get_dof_indices(local_dof_indices);
1567 local_dof_values.resize(element->dofs_per_cell);
1569 global_dof_values->extract_subvector_to(local_dof_indices,
1575 template <
int dim,
typename Number>
1578 const std::vector<unsigned int> &mask,
1581 for (
unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
1582 local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]];
1584 this->subcell_box = subcell_box;
1589 template <
int dim,
typename Number>
1598 template <
int dim,
typename Number>
1602 return n_subdivisions_per_line;
1607 template <
int dim,
typename Number>
1613 return local_dof_values.size() > 0;
1618 template <
int dim,
typename Number>
1622 const unsigned int component)
const
1627 if (!poly.empty() && component == 0)
1630 return this->is_fe_q_iso_q1() ?
1634 subcell_box.real_to_unit(point),
1635 polynomials_are_hat_functions) :
1640 polynomials_are_hat_functions,
1646 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1647 value += local_dof_values[i] *
1648 element->shape_value_component(i, point, component);
1656 template <
int dim,
typename Number>
1660 const unsigned int component)
const
1665 if (!poly.empty() && component == 0)
1668 return (this->is_fe_q_iso_q1() ?
1670 evaluate_tensor_product_value_and_gradient(
1673 subcell_box.real_to_unit(point),
1674 polynomials_are_hat_functions) :
1680 polynomials_are_hat_functions,
1687 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1688 gradient += local_dof_values[i] *
1689 element->shape_grad_component(i, point, component);
1697 template <
int dim,
typename Number>
1701 const unsigned int component)
const
1706 if (!poly.empty() && component == 0)
1709 return this->is_fe_q_iso_q1() ?
1713 subcell_box.real_to_unit(point)) :
1723 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1725 local_dof_values[i] *
1726 element->shape_grad_grad_component(i, point, component);
1737 const std::array<unsigned int, dim> &subcell_indices,
1738 const unsigned int n_subdivisions)
1741 for (
unsigned int d = 0; d < dim; ++d)
1745 const double split_box_side_length =
1747 corners.first[d] = subcell_indices[d] * split_box_side_length;
1748 corners.second[d] = corners.first[d] + split_box_side_length;
1757 std::vector<unsigned int>
1759 const std::array<unsigned int, dim> &subcell_indices,
1760 unsigned int dofs_per_line)
1765 std::vector<unsigned int> mask;
1766 mask.reserve(1 << dim);
1768 unsigned int stride = 1;
1769 for (
unsigned int d = 0; d < dim; ++d)
1771 const unsigned int previous_mask_size = 1U << d;
1773 for (
unsigned int i = 0; i < previous_mask_size; ++i)
1775 mask[i] += subcell_indices[d] * stride;
1776 mask.push_back(stride + mask[i]);
1778 stride *= dofs_per_line;
1788 const unsigned int max_box_splits,
1789 const double lower_bound_implicit_function,
1790 const double min_distance_between_roots,
1791 const double limit_to_be_definite,
1792 const double root_finder_tolerance,
1793 const unsigned int max_root_finder_splits,
1795 : max_box_splits(max_box_splits)
1796 , lower_bound_implicit_function(lower_bound_implicit_function)
1797 , min_distance_between_roots(min_distance_between_roots)
1798 , limit_to_be_definite(limit_to_be_definite)
1799 , root_finder_tolerance(root_finder_tolerance)
1800 , max_root_finder_splits(max_root_finder_splits)
1801 , split_in_half(split_in_half)
1810 : q_generator(q_collection, additional_data)
1813 ExcMessage(
"Incoming hp::QCollection<1> is empty."));
1822 q_generator.clear_quadratures();
1832 clear_quadratures();
1833 generate_append(level_set, box);
1845 "The incoming function should be a scalar level set function,"
1846 " it should have one component."));
1849 std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
1850 level_sets.push_back(level_set);
1852 const unsigned int n_box_splits = 0;
1853 q_generator.generate(level_sets, box, n_box_splits);
1859 q_generator.get_quadratures().indefinite.empty(),
1861 "Generation of quadrature rules failed. This can mean that the level "
1862 "set function is degenerate in some way, e.g. oscillating extremely "
1872 return q_generator.get_quadratures().negative;
1881 return q_generator.get_quadratures().positive;
1890 return q_generator.get_quadratures().surface;
1898 q_generator.set_1D_quadrature(q_index);
1907 : quadrature_generator(quadratures1D, additional_data)
1916 quadrature_generator.clear_quadratures();
1917 surface_quadrature.clear();
1926 const unsigned int face_index)
1928 clear_quadratures();
1929 generate_append(level_set, box, face_index);
1938 const unsigned int face_index)
1945 const unsigned int face_normal_direction =
1950 const double coordinate_value = vertex0[face_normal_direction];
1953 level_set, face_normal_direction, coordinate_value);
1959 quadrature_generator.generate_append(face_restriction, cross_section);
1964 &surface_quadrature_wrong_normal =
1965 quadrature_generator.get_surface_quadrature();
1967 std::vector<Tensor<1, dim>> normals;
1968 normals.reserve(surface_quadrature_wrong_normal.size());
1969 for (
unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1972 surface_quadrature_wrong_normal.point(i),
1973 face_normal_direction,
1977 normal /= normal.
norm();
1978 normals.push_back(normal);
1981 surface_quadrature_wrong_normal.get_points(),
1982 surface_quadrature_wrong_normal.get_weights(),
1992 quadrature_generator.set_1D_quadrature(q_index);
2001 return quadrature_generator.get_inside_quadrature();
2009 return quadrature_generator.get_outside_quadrature();
2018 return surface_quadrature;
2027 (void)quadratures1D;
2028 (void)additional_data;
2044 const unsigned int face_index)
2046 clear_quadratures();
2047 generate_append(level_set, box, face_index);
2055 const unsigned int face_index)
2063 const unsigned int n_points = 1;
2064 const double weight = 1;
2065 const std::vector<Point<0>> points(n_points);
2066 const std::vector<double> weights(n_points, weight);
2068 const double level_set_value = level_set.
value(vertex);
2069 if (level_set_value < 0)
2074 else if (level_set_value > 0)
2099 return inside_quadrature;
2106 return outside_quadrature;
2114 return surface_quadrature;
2120 template <
typename Number>
2127 , reference_space_level_set(
2128 std::make_unique<
internal::DiscreteQuadratureGeneratorImplementation::
2129 RefSpaceFEFieldFunction<dim, Number>>(dof_handler,
2140 const unsigned int n_subdivisions =
2141 reference_space_level_set->n_subdivisions();
2143 const unsigned int dofs_per_line = n_subdivisions + 1;
2145 this->clear_quadratures();
2154 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2155 all_subcell_indices.reserve(
Utilities::pow(n_subdivisions, dim));
2156 all_subcell_indices.emplace_back();
2157 for (
unsigned d = 0; d < dim; ++d)
2159 const unsigned int old_size = all_subcell_indices.size();
2160 for (
unsigned int i = 1; i < n_subdivisions; ++i)
2162 for (
unsigned int j = 0; j < old_size; ++j)
2164 std::array<unsigned int, dim> next_entry =
2165 all_subcell_indices[j];
2167 all_subcell_indices.push_back(next_entry);
2173 for (
const auto &subcell_indices : all_subcell_indices)
2175 const std::vector<unsigned int> lexicographic_mask =
2176 internal::DiscreteQuadratureGeneratorImplementation::
2177 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2179 const auto subcell_box =
2180 internal::DiscreteQuadratureGeneratorImplementation::
2181 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2183 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2197 Assert(cell->reference_cell().is_hyper_cube(),
2198 internal::DiscreteQuadratureGeneratorImplementation::
2199 ExcReferenceCellNotHypercube());
2201 reference_space_level_set->set_active_cell(cell);
2203 if (reference_space_level_set->is_fe_q_iso_q1())
2204 generate_fe_q_iso_q1(unit_box);
2212 template <
typename Number>
2219 , reference_space_level_set(
2220 std::make_unique<
internal::DiscreteQuadratureGeneratorImplementation::
2221 RefSpaceFEFieldFunction<dim, Number>>(dof_handler,
2231 unsigned int face_index)
2233 const unsigned int n_subdivisions =
2234 reference_space_level_set->n_subdivisions();
2236 const unsigned int dofs_per_line = n_subdivisions + 1;
2238 this->clear_quadratures();
2240 const unsigned int coordinate_direction_face = face_index / 2;
2241 const unsigned int coordinate_orientation_face = face_index % 2;
2253 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2254 all_subcell_indices.reserve(
Utilities::pow(n_subdivisions, dim));
2255 all_subcell_indices.emplace_back();
2256 all_subcell_indices[0][coordinate_direction_face] =
2257 coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1;
2258 for (
unsigned d = 0; d < dim - 1; ++d)
2260 const unsigned int old_size = all_subcell_indices.size();
2261 for (
unsigned int i = 1; i < n_subdivisions; ++i)
2263 for (
unsigned int j = 0; j < old_size; ++j)
2265 std::array<unsigned int, dim> next_entry =
2266 all_subcell_indices[j];
2267 next_entry[(coordinate_direction_face + d + 1) % dim] = i;
2268 all_subcell_indices.push_back(next_entry);
2273 for (
const auto &subcell_indices : all_subcell_indices)
2275 const std::vector<unsigned int> lexicographic_mask =
2276 internal::DiscreteQuadratureGeneratorImplementation::
2277 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2279 const auto subcell_box =
2280 internal::DiscreteQuadratureGeneratorImplementation::
2281 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2283 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2286 *reference_space_level_set, subcell_box, face_index);
2296 const unsigned int face_index)
2298 Assert(cell->reference_cell().is_hyper_cube(),
2299 internal::DiscreteQuadratureGeneratorImplementation::
2300 ExcReferenceCellNotHypercube());
2302 reference_space_level_set->set_active_cell(cell);
2304 if (reference_space_level_set->is_fe_q_iso_q1())
2305 generate_fe_q_iso_q1(unit_box, face_index);
2312#include "non_matching/quadrature_generator.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
BoundingBox< 1, Number > bounds(const unsigned int direction) const
Point< spacedim, Number > center() const
Number lower_bound(const unsigned int direction) const
std::pair< Point< spacedim, Number >, Point< spacedim, Number > > & get_boundary_points()
Number side_length(const unsigned int direction) const
BoundingBox< spacedim, Number > child(const unsigned int index) const
Point< spacedim, Number > vertex(const unsigned int index) const
BoundingBox< spacedim - 1, Number > cross_section(const unsigned int direction) const
Number upper_bound(const unsigned int direction) const
Point< spacedim, Number > unit_to_real(const Point< spacedim, Number > &point) const
virtual const FiniteElement< dim, spacedim > & base_element(const unsigned int index) const
unsigned int n_base_elements() const
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
const unsigned int n_components
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void generate(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index)
DiscreteFaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const ReadVector< Number > &level_set, const AdditionalData &additional_data=AdditionalData())
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box, unsigned int face_index)
void generate(const typename Triangulation< dim >::active_cell_iterator &cell)
DiscreteQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const ReadVector< Number > &level_set, const AdditionalData &additional_data=AdditionalData())
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box)
const Quadrature< dim - 1 > & get_inside_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
void set_1D_quadrature(const unsigned int q_index)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
FaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const ImmersedSurfaceQuadrature< dim - 1, dim > & get_surface_quadrature() const
const Quadrature< dim - 1 > & get_outside_quadrature() const
void push_back(const Point< dim > &point, const double weight, const Tensor< 1, spacedim > &normal)
void generate_append(const Function< dim > &level_set, const BoundingBox< dim > &box)
void set_1D_quadrature(const unsigned int q_index)
QuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim > & get_inside_quadrature() const
const Quadrature< dim > & get_outside_quadrature() const
const ImmersedSurfaceQuadrature< dim > & get_surface_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box)
Lazy< std::unique_ptr< FE_Q< dim > > > fe_q1
std::vector< unsigned int > renumber
std::vector< Number > local_dof_values_subcell
unsigned int n_subdivisions() const override
Tensor< 1, dim > gradient(const Point< dim > &point, const unsigned int component=0) const override
std::vector< types::global_dof_index > local_dof_indices
bool is_fe_q_iso_q1() const override
BoundingBox< dim > subcell_box
const ObserverPointer< const DoFHandler< dim > > dof_handler
unsigned int n_subdivisions_per_line
double value(const Point< dim > &point, const unsigned int component=0) const override
const ObserverPointer< const ReadVector< Number > > global_dof_values
SymmetricTensor< 2, dim > hessian(const Point< dim > &point, const unsigned int component=0) const override
void set_subcell(const std::vector< unsigned int > &mask, const BoundingBox< dim > &subcell_box) override
ObserverPointer< const FiniteElement< dim > > element
void set_active_cell(const typename Triangulation< dim >::active_cell_iterator &cell) override
std::vector< Polynomials::Polynomial< double > > poly
bool polynomials_are_hat_functions
std::vector< Number > local_dof_values
void push_back(const Point< dim > &point, const double weight)
ExtendableQuadrature()=default
QGeneratorBase(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
const ObserverPointer< const hp::QCollection< 1 > > q_collection1D
const QPartitioning< dim > & get_quadratures() const
QGenerator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
void set_1D_quadrature(const unsigned int q_index)
void create_high_dim_quadratures(const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box)
void generate(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int n_box_splits)
hp::QCollection< dim > tensor_products
void split_box_and_recurse(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const std::optional< HeightDirectionData > &direction_data, const unsigned int n_box_splits)
void use_midpoint_method(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box)
void create_low_dim_quadratures(const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int n_box_splits)
ExtendableQuadrature< dim > & quadrature_by_definiteness(const Definiteness definiteness)
ImmersedSurfaceQuadrature< dim > surface
const AdditionalData additional_data
void find_roots(const std::vector< std::reference_wrapper< const Function< 1 > > > &functions, const BoundingBox< 1 > &interval, std::vector< double > &roots)
RootFinder(const AdditionalData &data=AdditionalData())
void generate(const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const Quadrature< dim - 1 > &low_dim_quadrature, const unsigned int height_function_direction, QPartitioning< dim > &q_partitioning)
void set_1D_quadrature(const unsigned int q_index)
void create_surface_point(const Point< dim - 1 > &point, const double weight, const std::vector< std::reference_wrapper< const Function< dim > > > &level_sets, const BoundingBox< dim > &box, const unsigned int height_function_direction, ImmersedSurfaceQuadrature< dim > &surface_quadrature)
UpThroughDimensionCreator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
virtual size_type size() const =0
unsigned int n_vertices() const
numbers::NumberTraits< Number >::real_type norm() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcCellNotSet()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
std::vector< unsigned int > setup_lexicographic_mask(const std::array< unsigned int, dim > &subcell_indices, unsigned int dofs_per_line)
BoundingBox< dim > create_subcell_box(const BoundingBox< dim > &unit_box, const std::array< unsigned int, dim > &subcell_indices, const unsigned int n_subdivisions)
void distribute_points_between_roots(const Quadrature< 1 > &quadrature1D, const BoundingBox< 1 > &interval, const std::vector< double > &roots, const Point< dim - 1 > &point, const double weight, const unsigned int height_function_direction, const std::vector< std::reference_wrapper< const Function< 1 > > > &level_sets, const AdditionalQGeneratorData &additional_data, QPartitioning< dim > &q_partitioning)
double lower_bound_on_abs(const std::pair< double, double > &function_bounds)
void estimate_function_bounds(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const BoundingBox< dim > &box, std::vector< FunctionBounds< dim > > &all_function_bounds)
BoundingBox< dim > left_half(const BoundingBox< dim > &box, const unsigned int direction)
std::optional< unsigned int > direction_of_largest_extent(const BoundingBox< dim > &box)
void add_tensor_product(const Quadrature< dim - 1 > &lower, const Quadrature< 1 > &quadrature1D, const double start, const double end, const unsigned int component_in_dim, ExtendableQuadrature< dim > &quadrature)
bool one_positive_one_negative_definite(const std::vector< FunctionBounds< dim > > &all_function_bounds)
void map_quadrature_to_box(const Quadrature< dim > &unit_quadrature, const BoundingBox< dim > &box, ExtendableQuadrature< dim > &quadrature)
std::pair< double, double > find_extreme_values(const std::vector< FunctionBounds< dim > > &all_function_bounds)
std::optional< HeightDirectionData > find_best_height_direction(const std::vector< FunctionBounds< dim > > &all_function_bounds)
void tensor_point_with_1D_quadrature(const Point< dim - 1 > &point, const double weight, const Quadrature< 1 > &quadrature1D, const double start, const double end, const unsigned int component_in_dim, ExtendableQuadrature< dim > &quadrature)
void restrict_to_point(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const Point< dim - 1 > &point, const unsigned int open_direction, std::vector< Functions::PointRestriction< dim - 1 > > &restrictions)
BoundingBox< dim > right_half(const BoundingBox< dim > &box, const unsigned int direction)
bool is_indefinite(const std::pair< double, double > &function_bounds)
void take_min_max_at_vertices(const Function< dim > &function, const BoundingBox< dim > &box, std::pair< double, double > &value_bounds)
unsigned int compute_split_direction(const BoundingBox< dim > &box, const std::optional< HeightDirectionData > &height_direction_data)
void restrict_to_top_and_bottom(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const BoundingBox< dim > &box, const unsigned int direction, std::vector< Functions::CoordinateRestriction< dim - 1 > > &restrictions)
Point< dim > face_projection_closest_zero_contour(const Point< dim - 1 > &point, const unsigned int direction, const BoundingBox< dim > &box, const Function< dim > &level_set)
Definiteness pointwise_definiteness(const std::vector< std::reference_wrapper< const Function< dim > > > &functions, const Point< dim > &point)
constexpr T pow(const T base, const int iexp)
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const std::vector< unsigned int > &renumber={})
std::pair< typename ProductTypeNoPoint< Number, Number2 >::type, Tensor< 1, dim, typename ProductTypeNoPoint< Number, Number2 >::type > > evaluate_tensor_product_value_and_gradient(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
Point< dim+1 > create_higher_dim_point(const Point< dim > &point, const unsigned int component_in_dim_plus_1, const double coordinate_value)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value(const std::vector< Polynomials::Polynomial< double > > &poly, const ArrayView< const Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
constexpr unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
double min_distance_between_roots
AdditionalQGeneratorData(const unsigned int max_box_splits=4, const double lower_bound_implicit_function=1e-11, const double min_distance_between_roots=1e-12, const double limit_to_be_definite=1e-11, const double root_finder_tolerance=1e-12, const unsigned int max_root_finder_splits=2, bool split_in_half=true)
std::array< std::pair< double, double >, dim > gradient
std::pair< double, double > value
AdditionalData(const double tolerance=1e-12, const unsigned int max_recursion_depth=2, const unsigned int max_iterations=500)
unsigned int max_iterations
unsigned int max_recursion_depth
void reinit(const Quadrature< dim_q > &quad, const FiniteElement< dim, spacedim > &fe_dim, const unsigned int base_element=0)
std::vector< unsigned int > lexicographic_numbering
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)