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>>>
114 "The incoming vector must contain at least one function."));
116 const int sign_of_first =
119 if (sign_of_first == 0)
122 for (
unsigned int j = 1; j <
functions.size(); ++j)
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);
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());
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();
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;
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)
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));
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 =
610 if (function_min > 0)
613 const double function_max =
617 if (function_max < 0)
652 this->weights.clear();
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);
726 return std::abs(bottom_value) < std::abs(top_value) ? bottom_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)
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 =
952 if (data && data->min_abs_dfdx >
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>
1180 const std::vector<std::reference_wrapper<
const Function<dim>>>
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))
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);
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;
1272 const double gradient_norm = normal.
norm();
1274 gradient_norm > 1
e-11,
1276 "The level set function has a gradient almost equal to 0."));
1277 normal *= 1. / gradient_norm;
1285 template <
int spacedim>
1299 namespace DiscreteQuadratureGeneratorImplementation
1303 "The set_active_cell function has to be called before calling this function.");
1307 "The reference cell of the incoming cell must be a hypercube.");
1335 template <
int dim,
typename VectorType = Vector<
double>>
1347 const VectorType &dof_values);
1360 set_subcell(
const std::vector<unsigned int> &mask,
1384 const unsigned int component = 0)
const override;
1395 const unsigned int component = 0)
const override;
1406 const unsigned int component = 0)
const override;
1465 std::vector<Polynomials::Polynomial<double>>
poly;
1480 template <
int dim,
typename VectorType>
1483 const VectorType &dof_values)
1484 : dof_handler(&dof_handler)
1485 , global_dof_values(&dof_values)
1494 template <
int dim,
typename VectorType>
1500 &cell->get_triangulation() == &dof_handler->get_triangulation(),
1502 "The incoming cell must belong to the triangulation associated with "
1503 "the DoFHandler passed to the constructor."));
1505 const auto dof_handler_cell =
1506 cell->as_dof_handler_iterator(*dof_handler);
1516 if (element != &dof_handler_cell->get_fe())
1519 element = &dof_handler_cell->get_fe();
1525 &dof_handler_cell->get_fe()))
1527 this->n_subdivisions_per_line = fe_q_iso_q1->get_degree();
1547 polynomials_are_hat_functions =
1548 (poly.size() == 2 && poly[0].value(0.) == 1. &&
1549 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1550 poly[1].value(1.) == 1.);
1554 element = &dof_handler_cell->get_fe();
1556 local_dof_indices.resize(element->dofs_per_cell);
1557 dof_handler_cell->get_dof_indices(local_dof_indices);
1559 local_dof_values.resize(element->dofs_per_cell);
1561 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1562 local_dof_values[i] =
1564 *global_dof_values, local_dof_indices[i]);
1569 template <
int dim,
typename VectorType>
1572 const std::vector<unsigned int> &mask,
1575 for (
unsigned int i = 0; i < local_dof_values_subcell.size(); ++i)
1576 local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]];
1578 this->subcell_box = subcell_box;
1583 template <
int dim,
typename VectorType>
1592 template <
int dim,
typename VectorType>
1596 return n_subdivisions_per_line;
1601 template <
int dim,
typename VectorType>
1607 return local_dof_values.size() > 0;
1612 template <
int dim,
typename VectorType>
1616 const unsigned int component)
const
1621 if (!poly.empty() && component == 0)
1626 this->is_fe_q_iso_q1() ? local_dof_values_subcell :
1628 this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(
point) :
point,
1629 polynomials_are_hat_functions,
1630 this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
1635 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1636 value += local_dof_values[i] *
1637 element->shape_value_component(i,
point, component);
1645 template <
int dim,
typename VectorType>
1649 const unsigned int component)
const
1654 if (!poly.empty() && component == 0)
1659 this->is_fe_q_iso_q1() ? local_dof_values_subcell :
1661 this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(
point) :
1663 polynomials_are_hat_functions,
1664 this->is_fe_q_iso_q1() ? std::vector<unsigned int>() :
1671 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1672 gradient += local_dof_values[i] *
1673 element->shape_grad_component(i,
point, component);
1681 template <
int dim,
typename VectorType>
1685 const unsigned int component)
const
1690 if (!poly.empty() && component == 0)
1695 this->is_fe_q_iso_q1() ? local_dof_values_subcell :
1697 this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(
point) :
point,
1698 this->is_fe_q_iso_q1() ? std::vector<unsigned int>() : renumber);
1703 for (
unsigned int i = 0; i < local_dof_indices.size(); ++i)
1705 local_dof_values[i] *
1706 element->shape_grad_grad_component(i,
point, component);
1717 const std::array<unsigned int, dim> &subcell_indices,
1718 const unsigned int n_subdivisions)
1721 for (
unsigned int d = 0;
d < dim; ++
d)
1725 const double split_box_side_length =
1727 corners.first[
d] = subcell_indices[
d] * split_box_side_length;
1728 corners.second[
d] = corners.first[
d] + split_box_side_length;
1737 std::vector<unsigned int>
1739 const std::array<unsigned int, dim> &subcell_indices,
1740 unsigned int dofs_per_line)
1745 std::vector<unsigned int> mask;
1746 mask.reserve(1 << dim);
1748 unsigned int stride = 1;
1749 for (
unsigned int d = 0;
d < dim; ++
d)
1751 const unsigned int previous_mask_size = 1U <<
d;
1753 for (
unsigned int i = 0; i < previous_mask_size; ++i)
1755 mask[i] += subcell_indices[
d] * stride;
1756 mask.push_back(stride + mask[i]);
1758 stride *= dofs_per_line;
1768 const unsigned int max_box_splits,
1769 const double lower_bound_implicit_function,
1770 const double min_distance_between_roots,
1771 const double limit_to_be_definite,
1772 const double root_finder_tolerance,
1773 const unsigned int max_root_finder_splits,
1775 : max_box_splits(max_box_splits)
1776 , lower_bound_implicit_function(lower_bound_implicit_function)
1777 , min_distance_between_roots(min_distance_between_roots)
1778 , limit_to_be_definite(limit_to_be_definite)
1779 , root_finder_tolerance(root_finder_tolerance)
1780 , max_root_finder_splits(max_root_finder_splits)
1781 , split_in_half(split_in_half)
1790 : q_generator(q_collection, additional_data)
1793 ExcMessage(
"Incoming hp::QCollection<1> is empty."));
1802 q_generator.clear_quadratures();
1812 clear_quadratures();
1813 generate_append(level_set, box);
1825 "The incoming function should be a scalar level set function,"
1826 " it should have one component."));
1829 std::vector<std::reference_wrapper<const Function<dim>>> level_sets;
1830 level_sets.push_back(level_set);
1832 const unsigned int n_box_splits = 0;
1833 q_generator.generate(level_sets, box, n_box_splits);
1839 q_generator.get_quadratures().indefinite.empty(),
1841 "Generation of quadrature rules failed. This can mean that the level "
1842 "set function is degenerate in some way, e.g. oscillating extremely "
1852 return q_generator.get_quadratures().negative;
1861 return q_generator.get_quadratures().positive;
1870 return q_generator.get_quadratures().surface;
1878 q_generator.set_1D_quadrature(q_index);
1887 : quadrature_generator(quadratures1D, additional_data)
1896 quadrature_generator.clear_quadratures();
1897 surface_quadrature.clear();
1906 const unsigned int face_index)
1908 clear_quadratures();
1909 generate_append(level_set, box, face_index);
1918 const unsigned int face_index)
1925 const unsigned int face_normal_direction =
1930 const double coordinate_value = vertex0(face_normal_direction);
1933 level_set, face_normal_direction, coordinate_value);
1939 quadrature_generator.generate_append(face_restriction, cross_section);
1944 &surface_quadrature_wrong_normal =
1945 quadrature_generator.get_surface_quadrature();
1947 std::vector<Tensor<1, dim>> normals;
1948 normals.reserve(surface_quadrature_wrong_normal.size());
1949 for (
unsigned int i = 0; i < surface_quadrature_wrong_normal.size(); ++i)
1952 surface_quadrature_wrong_normal.point(i),
1953 face_normal_direction,
1957 normal /= normal.
norm();
1958 normals.push_back(normal);
1961 surface_quadrature_wrong_normal.get_points(),
1962 surface_quadrature_wrong_normal.get_weights(),
1972 quadrature_generator.set_1D_quadrature(q_index);
1981 return quadrature_generator.get_inside_quadrature();
1989 return quadrature_generator.get_outside_quadrature();
1998 return surface_quadrature;
2007 (void)quadratures1D;
2008 (void)additional_data;
2024 const unsigned int face_index)
2035 const unsigned int face_index)
2043 const unsigned int n_points = 1;
2044 const double weight = 1;
2045 const std::vector<Point<0>> points(n_points);
2046 const std::vector<double> weights(n_points, weight);
2048 const double level_set_value = level_set.
value(vertex);
2049 if (level_set_value < 0)
2054 else if (level_set_value > 0)
2079 return inside_quadrature;
2086 return outside_quadrature;
2100 template <
typename VectorType>
2104 const VectorType &level_set,
2107 , reference_space_level_set(
2108 std::make_unique<
internal::DiscreteQuadratureGeneratorImplementation::
2109 RefSpaceFEFieldFunction<dim, VectorType>>(
2121 const unsigned int n_subdivisions =
2122 reference_space_level_set->n_subdivisions();
2124 const unsigned int dofs_per_line = n_subdivisions + 1;
2126 this->clear_quadratures();
2135 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2136 all_subcell_indices.reserve(
Utilities::pow(n_subdivisions, dim));
2137 all_subcell_indices.emplace_back();
2138 for (
unsigned d = 0;
d < dim; ++
d)
2140 const unsigned int old_size = all_subcell_indices.size();
2141 for (
unsigned int i = 1; i < n_subdivisions; ++i)
2143 for (
unsigned int j = 0; j < old_size; ++j)
2145 std::array<unsigned int, dim> next_entry =
2146 all_subcell_indices[j];
2148 all_subcell_indices.push_back(next_entry);
2154 for (
const auto &subcell_indices : all_subcell_indices)
2156 const std::vector<unsigned int> lexicographic_mask =
2157 internal::DiscreteQuadratureGeneratorImplementation::
2158 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2160 const auto subcell_box =
2161 internal::DiscreteQuadratureGeneratorImplementation::
2162 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2164 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2178 Assert(cell->reference_cell().is_hyper_cube(),
2179 internal::DiscreteQuadratureGeneratorImplementation::
2182 reference_space_level_set->set_active_cell(cell);
2184 if (reference_space_level_set->is_fe_q_iso_q1())
2185 generate_fe_q_iso_q1(unit_box);
2193 template <
typename VectorType>
2197 const VectorType &level_set,
2200 , reference_space_level_set(
2201 std::make_unique<
internal::DiscreteQuadratureGeneratorImplementation::
2202 RefSpaceFEFieldFunction<dim, VectorType>>(
2213 unsigned int face_index)
2215 const unsigned int n_subdivisions =
2216 reference_space_level_set->n_subdivisions();
2218 const unsigned int dofs_per_line = n_subdivisions + 1;
2220 this->clear_quadratures();
2222 const unsigned int coordinate_direction_face = face_index / 2;
2223 const unsigned int coordinate_orientation_face = face_index % 2;
2235 std::vector<std::array<unsigned int, dim>> all_subcell_indices;
2236 all_subcell_indices.reserve(
Utilities::pow(n_subdivisions, dim));
2237 all_subcell_indices.emplace_back();
2238 all_subcell_indices[0][coordinate_direction_face] =
2239 coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1;
2240 for (
unsigned d = 0;
d < dim - 1; ++
d)
2242 const unsigned int old_size = all_subcell_indices.size();
2243 for (
unsigned int i = 1; i < n_subdivisions; ++i)
2245 for (
unsigned int j = 0; j < old_size; ++j)
2247 std::array<unsigned int, dim> next_entry =
2248 all_subcell_indices[j];
2249 next_entry[(coordinate_direction_face +
d + 1) % dim] = i;
2250 all_subcell_indices.push_back(next_entry);
2255 for (
const auto &subcell_indices : all_subcell_indices)
2257 const std::vector<unsigned int> lexicographic_mask =
2258 internal::DiscreteQuadratureGeneratorImplementation::
2259 setup_lexicographic_mask<dim>(subcell_indices, dofs_per_line);
2261 const auto subcell_box =
2262 internal::DiscreteQuadratureGeneratorImplementation::
2263 create_subcell_box<dim>(unit_box, subcell_indices, n_subdivisions);
2265 reference_space_level_set->set_subcell(lexicographic_mask, subcell_box);
2268 *reference_space_level_set, subcell_box, face_index);
2278 const unsigned int face_index)
2280 Assert(cell->reference_cell().is_hyper_cube(),
2281 internal::DiscreteQuadratureGeneratorImplementation::
2284 reference_space_level_set->set_active_cell(cell);
2286 if (reference_space_level_set->is_fe_q_iso_q1())
2287 generate_fe_q_iso_q1(unit_box, face_index);
2294 #include "quadrature_generator.inst"
BoundingBox< 1, Number > bounds(const unsigned int direction) const
Point< spacedim, Number > center() const
Number lower_bound(const unsigned int direction) const
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
std::pair< Point< spacedim, Number >, Point< spacedim, Number > > & get_boundary_points()
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
unsigned int n_dofs_per_cell() 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
DiscreteFaceQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
void generate(const typename Triangulation< dim >::active_cell_iterator &cell, const unsigned int face_index)
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)
void generate_fe_q_iso_q1(const BoundingBox< dim > &unit_box)
DiscreteQuadratureGenerator(const hp::QCollection< 1 > &quadratures1D, const DoFHandler< dim > &dof_handler, const VectorType &level_set, const AdditionalData &additional_data=AdditionalData())
const Quadrature< dim - 1 > & get_inside_quadrature() const
void generate(const Function< dim > &level_set, const BoundingBox< dim > &box, const unsigned int face_index)
ImmersedSurfaceQuadrature< dim - 1, dim > surface_quadrature
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)
void set_active_cell(const typename Triangulation< dim >::active_cell_iterator &cell) override
bool is_fe_q_iso_q1() const override
Tensor< 1, dim > gradient(const Point< dim > &point, const unsigned int component=0) const override
double value(const Point< dim > &point, const unsigned int component=0) const override
std::vector< typename VectorType::value_type > local_dof_values_subcell
std::vector< types::global_dof_index > local_dof_indices
void set_subcell(const std::vector< unsigned int > &mask, const BoundingBox< dim > &subcell_box) override
unsigned int n_subdivisions_per_line
BoundingBox< dim > subcell_box
bool polynomials_are_hat_functions
RefSpaceFEFieldFunction(const DoFHandler< dim > &dof_handler, const VectorType &dof_values)
const SmartPointer< const DoFHandler< dim > > dof_handler
std::vector< Polynomials::Polynomial< double > > poly
SymmetricTensor< 2, dim > hessian(const Point< dim > &point, const unsigned int component=0) const override
unsigned int n_subdivisions() const override
SmartPointer< const FiniteElement< dim > > element
std::vector< unsigned int > renumber
std::vector< typename VectorType::value_type > local_dof_values
const SmartPointer< const VectorType > global_dof_values
void push_back(const Point< dim > &point, const double weight)
ExtendableQuadrature()=default
const SmartPointer< const hp::QCollection< 1 > > q_collection1D
const AdditionalQGeneratorData additional_data
QGeneratorBase(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
const QPartitioning< dim > & get_quadratures() const
QPartitioning< dim > q_partitioning
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)
QGenerator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
void set_1D_quadrature(const unsigned int q_index)
hp::QCollection< dim > tensor_products
void generate(const std::vector< std::reference_wrapper< const Function< dim >>> &level_sets, const BoundingBox< dim > &box, 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_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 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 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)
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)
UpThroughDimensionCreator(const hp::QCollection< 1 > &q_collection1D, const AdditionalQGeneratorData &additional_data)
double weight(const unsigned int i) const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
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
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcReferenceCellNotHypercube()
static ::ExceptionBase & ExcCellNotSet()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ExcMessage(std::string arg1)
Expression sign(const Expression &x)
BoundingBox< dim > create_subcell_box(const BoundingBox< dim > &unit_box, const std::array< unsigned int, dim > &subcell_indices, const unsigned int n_subdivisions)
std::vector< unsigned int > setup_lexicographic_mask(const std::array< unsigned int, dim > &subcell_indices, unsigned int dofs_per_line)
std::optional< HeightDirectionData > find_best_height_direction(const std::vector< FunctionBounds< dim >> &all_function_bounds)
std::pair< double, double > find_extreme_values(const std::vector< FunctionBounds< dim >> &all_function_bounds)
double lower_bound_on_abs(const std::pair< double, double > &function_bounds)
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)
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)
BoundingBox< dim > right_half(const BoundingBox< dim > &box, const unsigned int direction)
void map_quadrature_to_box(const Quadrature< dim > &unit_quadrature, const BoundingBox< dim > &box, ExtendableQuadrature< dim > &quadrature)
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)
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)
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)
bool is_indefinite(const std::pair< double, double > &function_bounds)
BoundingBox< dim > left_half(const BoundingBox< dim > &box, const unsigned int direction)
Definiteness pointwise_definiteness(const std::vector< std::reference_wrapper< const Function< dim >>> &functions, const Point< dim > &point)
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)
std::optional< unsigned int > direction_of_largest_extent(const BoundingBox< dim > &box)
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)
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)
bool one_positive_one_negative_definite(const std::vector< FunctionBounds< dim >> &all_function_bounds)
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)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
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)
int(&) functions(const void *v1, const void *v2)
SymmetricTensor< 2, dim, typename ProductTypeNoPoint< Number, Number2 >::type > evaluate_tensor_product_hessian(const std::vector< Polynomials::Polynomial< double >> &poly, const std::vector< 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 std::vector< 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 std::vector< Number > &values, const Point< dim, Number2 > &p, const bool d_linear=false, const std::vector< unsigned int > &renumber={})
static const unsigned int invalid_unsigned_int
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
constexpr DEAL_II_HOST SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)