30#include <boost/container/small_vector.hpp>
53 const double theta = dir.
norm();
63 return tmp / tmp.
norm();
80 template <
int spacedim>
92 ExcMessage(
"The direction parameter must not be zero!"));
99 normal[0] = (vector[1] + vector[2]) / vector[0];
106 normal[1] = (vector[0] + vector[2]) / vector[1];
112 normal[2] = (vector[0] + vector[1]) / vector[2];
115 normal /= normal.
norm();
126template <
int dim,
int spacedim>
136template <
int dim,
int spacedim>
137std::unique_ptr<Manifold<dim, spacedim>>
140 return std::make_unique<PolarManifold<dim, spacedim>>(*this);
145template <
int dim,
int spacedim>
162template <
int dim,
int spacedim>
167 Assert(spherical_point[0] >= 0.0,
168 ExcMessage(
"Negative radius for given point."));
169 const double rho = spherical_point[0];
170 const double theta = spherical_point[1];
182 const double phi = spherical_point[2];
196template <
int dim,
int spacedim>
202 const double rho = R.
norm();
211 p[1] = std::atan2(R[1], R[0]);
219 const double z = R[2];
220 p[2] = std::atan2(R[1], R[0]);
223 p[1] = std::atan2(
std::sqrt(R[0] * R[0] + R[1] * R[1]), z);
235template <
int dim,
int spacedim>
240 Assert(spherical_point[0] >= 0.0,
241 ExcMessage(
"Negative radius for given point."));
242 const double rho = spherical_point[0];
243 const double theta = spherical_point[1];
260 const double phi = spherical_point[2];
285 template <
int dim,
int spacedim>
287 spherical_face_is_horizontal(
297 constexpr unsigned int n_vertices =
299 std::array<double, n_vertices> sqr_distances_to_center;
300 std::array<double, n_vertices - 1> sqr_distances_to_first_vertex;
301 sqr_distances_to_center[0] =
302 (face->vertex(0) - manifold_center).norm_square();
303 for (
unsigned int i = 1; i < n_vertices; ++i)
305 sqr_distances_to_center[i] =
306 (face->vertex(i) - manifold_center).norm_square();
307 sqr_distances_to_first_vertex[i - 1] =
308 (face->vertex(i) - face->vertex(0)).norm_square();
310 const auto minmax_sqr_distance =
311 std::minmax_element(sqr_distances_to_center.begin(),
312 sqr_distances_to_center.end());
313 const auto min_sqr_distance_to_first_vertex =
314 std::min_element(sqr_distances_to_first_vertex.begin(),
315 sqr_distances_to_first_vertex.end());
317 return (*minmax_sqr_distance.second - *minmax_sqr_distance.first <
318 1.e-10 * *min_sqr_distance_to_first_vertex);
324template <
int dim,
int spacedim>
334 if (spherical_face_is_horizontal<dim, spacedim>(face, center))
341 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
342 return normalized_spherical_normal;
358template <
int dim,
int spacedim>
363 , polar_manifold(center)
368template <
int dim,
int spacedim>
369std::unique_ptr<Manifold<dim, spacedim>>
372 return std::make_unique<SphericalManifold<dim, spacedim>>(*this);
377template <
int dim,
int spacedim>
382 const double w)
const
384 const double tol = 1e-10;
386 if ((p1 - p2).norm_square() < tol * tol ||
std::abs(w) < tol)
398 const double r1 =
v1.norm();
399 const double r2 = v2.
norm();
401 Assert(r1 > tol && r2 > tol,
402 ExcMessage(
"p1 and p2 cannot coincide with the center."));
408 const double cosgamma = e1 * e2;
411 if (cosgamma < -1 + 8. * std::numeric_limits<double>::epsilon())
415 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
421 const double sigma = w *
std::acos(cosgamma);
426 const double n_norm = n.
norm();
429 "Probably, this means v1==v2 or v2==0."));
443template <
int dim,
int spacedim>
449 [[maybe_unused]]
const double tol = 1e-10;
455 const double r1 =
v1.norm();
456 const double r2 = v2.
norm();
466 const double cosgamma = e1 * e2;
468 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
469 ExcMessage(
"p1 and p2 cannot lie on the same diameter and be opposite "
470 "respect to the center."));
472 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
478 const double n_norm = n.
norm();
481 "Probably, this means v1==v2 or v2==0."));
487 const double gamma =
std::acos(cosgamma);
488 return (r2 - r1) * e1 + r1 * gamma * n;
493template <
int dim,
int spacedim>
503 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
510 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
511 return normalized_spherical_normal;
545template <
int dim,
int spacedim>
556 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
561 for (
unsigned int vertex = 0;
562 vertex < GeometryInfo<spacedim>::vertices_per_face;
564 face_vertex_normals[vertex] = face->vertex(vertex) - p_center;
572template <
int dim,
int spacedim>
582 do_get_new_points(surrounding_points,
make_array_view(weights), new_points);
589template <
int dim,
int spacedim>
598 do_get_new_points(vertices,
613 template <
int spacedim>
636 Point<3> candidate = candidate_point;
637 const unsigned int n_merged_points = directions.size();
638 const double tolerance = 1e-10;
639 const int max_iterations = 10;
644 for (
unsigned int i = 0; i < n_merged_points; ++i)
646 const double squared_distance =
647 (candidate - directions[i]).norm_square();
648 if (squared_distance < tolerance * tolerance)
654 if (n_merged_points == 2)
656 static const ::SphericalManifold<3, 3> unit_manifold;
660 unit_manifold.get_intermediate_point(
Point<3>(directions[0]),
673 for (
unsigned int i = 0; i < max_iterations; ++i)
679 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
687 for (
unsigned int i = 0; i < n_merged_points; ++i)
693 const double sintheta =
std::sqrt(sinthetaSq);
694 if (sintheta < tolerance)
696 Hessian[0][0] += weights[i];
697 Hessian[1][1] += weights[i];
701 const double costheta = (directions[i]) * candidate;
702 const double theta = std::atan2(sintheta, costheta);
703 const double sincthetaInv = theta / sintheta;
705 const double cosphi = vPerp * Clocalx;
706 const double sinphi = vPerp * Clocaly;
708 gradlocal[0] = cosphi;
709 gradlocal[1] = sinphi;
710 gradient += (weights[i] * sincthetaInv) * gradlocal;
712 const double wt = weights[i] / sinthetaSq;
713 const double sinphiSq = sinphi * sinphi;
714 const double cosphiSq = cosphi * cosphi;
715 const double tt = sincthetaInv * costheta;
716 const double offdiag =
717 cosphi * sinphi * wt * (1.0 - tt);
718 Hessian[0][0] += wt * (cosphiSq + tt * sinphiSq);
719 Hessian[0][1] += offdiag;
720 Hessian[1][0] += offdiag;
721 Hessian[1][1] += wt * (sinphiSq + tt * cosphiSq);
731 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
735 const Point<3> candidateOld = candidate;
740 if ((candidate - candidateOld).norm_square() <
741 tolerance * tolerance)
753template <
int dim,
int spacedim>
761 new_points.size() * surrounding_points.size());
762 const unsigned int weight_rows = new_points.size();
763 const unsigned int weight_columns = surrounding_points.size();
765 if (surrounding_points.size() == 2)
767 for (
unsigned int row = 0; row < weight_rows; ++row)
770 surrounding_points[0],
771 surrounding_points[1],
772 weights[row * weight_columns + 1]);
776 boost::container::small_vector<std::pair<double, Tensor<1, spacedim>>, 100>
777 new_candidates(new_points.size());
778 boost::container::small_vector<Tensor<1, spacedim>, 100> directions(
780 boost::container::small_vector<double, 100> distances(
781 surrounding_points.size(), 0.0);
782 double max_distance = 0.;
783 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
785 directions[i] = surrounding_points[i] - p_center;
786 distances[i] = directions[i].norm();
788 if (distances[i] != 0.)
789 directions[i] /= distances[i];
792 ExcMessage(
"One of the vertices coincides with the center. "
793 "This is not allowed!"));
797 for (
unsigned int k = 0; k < i; ++k)
799 const double squared_distance =
800 (directions[i] - directions[k]).norm_square();
801 max_distance =
std::max(max_distance, squared_distance);
807 const double tolerance = 1e-10;
808 boost::container::small_vector<bool, 100> accurate_point_was_found(
809 new_points.size(),
false);
814 for (
unsigned int row = 0; row < weight_rows; ++row)
816 new_candidates[row] =
817 guess_new_point(array_directions,
824 if (new_candidates[row].
first == 0.0)
826 new_points[row] = p_center;
827 accurate_point_was_found[row] =
true;
834 new_points[row] = polar_manifold.get_new_point(
842 if constexpr (spacedim < 3)
850 if (max_distance < 2e-2)
852 for (
unsigned int row = 0; row < weight_rows; ++row)
854 p_center + new_candidates[row].
first * new_candidates[row].
second;
864 boost::container::small_vector<double, 1000> merged_weights(
866 boost::container::small_vector<Tensor<1, spacedim>, 100>
868 boost::container::small_vector<double, 100> merged_distances(
869 surrounding_points.size(), 0.0);
871 unsigned int n_unique_directions = 0;
872 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
874 bool found_duplicate =
false;
878 for (
unsigned int j = 0; j < n_unique_directions; ++j)
880 const double squared_distance =
881 (directions[i] - directions[j]).norm_square();
882 if (!found_duplicate && squared_distance < 1e-28)
884 found_duplicate =
true;
885 for (
unsigned int row = 0; row < weight_rows; ++row)
886 merged_weights[row * weight_columns + j] +=
887 weights[row * weight_columns + i];
891 if (found_duplicate ==
false)
893 merged_directions[n_unique_directions] = directions[i];
894 merged_distances[n_unique_directions] = distances[i];
895 for (
unsigned int row = 0; row < weight_rows; ++row)
896 merged_weights[row * weight_columns + n_unique_directions] =
897 weights[row * weight_columns + i];
899 ++n_unique_directions;
905 boost::container::small_vector<unsigned int, 100> merged_weights_index(
907 for (
unsigned int row = 0; row < weight_rows; ++row)
909 for (
unsigned int existing_row = 0; existing_row < row;
912 bool identical_weights =
true;
914 for (
unsigned int weight_index = 0;
915 weight_index < n_unique_directions;
918 merged_weights[row * weight_columns + weight_index] -
919 merged_weights[existing_row * weight_columns +
920 weight_index]) > tolerance)
922 identical_weights =
false;
926 if (identical_weights)
928 merged_weights_index[row] = existing_row;
938 merged_directions.begin() + n_unique_directions);
941 merged_distances.begin() + n_unique_directions);
943 for (
unsigned int row = 0; row < weight_rows; ++row)
944 if (!accurate_point_was_found[row])
949 &merged_weights[row * weight_columns], n_unique_directions);
950 new_candidates[row].second =
951 internal::SphericalManifold::do_get_new_point(
952 array_merged_directions,
953 array_merged_distances,
954 array_merged_weights,
958 new_candidates[row].second =
959 new_candidates[merged_weights_index[row]].second;
962 p_center + new_candidates[row].first * new_candidates[row].second;
969template <
int dim,
int spacedim>
970std::pair<double, Tensor<1, spacedim>>
976 const double tolerance = 1e-10;
981 double total_weights = 0.;
982 for (
unsigned int i = 0; i < directions.size(); ++i)
985 if (
std::abs(1 - weights[i]) < tolerance)
986 return std::make_pair(distances[i], directions[i]);
988 rho += distances[i] * weights[i];
989 candidate += directions[i] * weights[i];
990 total_weights += weights[i];
994 const double norm = candidate.
norm();
998 rho /= total_weights;
1000 return std::make_pair(rho, candidate);
1008template <
int dim,
int spacedim>
1010 const double tolerance)
1018 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1023template <
int dim,
int spacedim>
1027 const double tolerance)
1030 , direction(direction / direction.norm())
1031 , point_on_axis(point_on_axis)
1032 , tolerance(tolerance)
1033 , dxn(cross_product_3d(this->direction, normal_direction))
1038 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1043template <
int dim,
int spacedim>
1044std::unique_ptr<Manifold<dim, spacedim>>
1047 return std::make_unique<CylindricalManifold<dim, spacedim>>(*this);
1052template <
int dim,
int spacedim>
1059 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1063 double average_length = 0.;
1064 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
1066 middle += surrounding_points[i] * weights[i];
1067 average_length += surrounding_points[i].square() * weights[i];
1069 middle -= point_on_axis;
1070 const double lambda = middle * direction;
1072 if ((middle - direction * lambda).square() < tolerance * average_length)
1073 return point_on_axis + direction * lambda;
1081template <
int dim,
int spacedim>
1087 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1091 const double lambda = normalized_point * direction;
1092 const Point<spacedim> projection = point_on_axis + direction * lambda;
1107template <
int dim,
int spacedim>
1113 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1116 const double sine_r =
std::sin(chart_point[1]) * chart_point[0];
1117 const double cosine_r =
std::cos(chart_point[1]) * chart_point[0];
1119 normal_direction * cosine_r + dxn * sine_r;
1122 return point_on_axis + direction * chart_point[2] + intermediate;
1127template <
int dim,
int spacedim>
1133 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1138 const double sine =
std::sin(chart_point[1]);
1139 const double cosine =
std::cos(chart_point[1]);
1141 normal_direction * cosine + dxn * sine;
1144 constexpr int s0 = 0 % spacedim;
1145 constexpr int s1 = 1 % spacedim;
1146 constexpr int s2 = 2 % spacedim;
1149 derivatives[s0][s0] = intermediate[s0];
1150 derivatives[s1][s0] = intermediate[s1];
1151 derivatives[s2][s0] = intermediate[s2];
1154 derivatives[s0][s1] = -normal_direction[s0] * sine + dxn[s0] * cosine;
1155 derivatives[s1][s1] = -normal_direction[s1] * sine + dxn[s1] * cosine;
1156 derivatives[s2][s1] = -normal_direction[s2] * sine + dxn[s2] * cosine;
1159 derivatives[s0][s2] = direction[s0];
1160 derivatives[s1][s2] = direction[s1];
1161 derivatives[s2][s2] = direction[s2];
1174 const double norm = t.
norm();
1175 Assert(norm > 0.0,
ExcMessage(
"The major axis must have a positive norm."));
1185template <
int dim,
int spacedim>
1189 const double eccentricity)
1192 , direction(check_and_normalize(major_axis_direction))
1194 , eccentricity(eccentricity)
1195 , cosh_u(1.0 / eccentricity)
1196 , sinh_u(
std::sqrt(cosh_u * cosh_u - 1.0))
1203 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1208template <
int dim,
int spacedim>
1209std::unique_ptr<Manifold<dim, spacedim>>
1212 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1217template <
int dim,
int spacedim>
1230template <
int dim,
int spacedim>
1244 const double cs =
std::cos(chart_point[1]);
1245 const double sn =
std::sin(chart_point[1]);
1248 const double x = chart_point[0] * cosh_u * cs;
1249 const double y = chart_point[0] * sinh_u * sn;
1251 const Point<2> p(direction[0] * x - direction[1] * y,
1252 direction[1] * x + direction[0] * y);
1258template <
int dim,
int spacedim>
1273 const double x0 = space_point[0] - center[0];
1274 const double y0 = space_point[1] - center[1];
1275 const double x = direction[0] * x0 + direction[1] * y0;
1276 const double y = -direction[1] * x0 + direction[0] * y0;
1278 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1284 double cos_eta = x / (pt0 * cosh_u);
1294 const double pt1 = (std::signbit(y) ? 2.0 *
numbers::PI - eta : eta);
1300template <
int dim,
int spacedim>
1316 const double cs =
std::cos(chart_point[1]);
1317 const double sn =
std::sin(chart_point[1]);
1319 dX[0][0] = cosh_u * cs;
1320 dX[0][1] = -chart_point[0] * cosh_u * sn;
1321 dX[1][0] = sinh_u * sn;
1322 dX[1][1] = chart_point[0] * sinh_u * cs;
1326 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1336template <
int dim,
int spacedim,
int chartdim>
1341 const double tolerance)
1344 , push_forward_function(&push_forward_function)
1345 , pull_back_function(&pull_back_function)
1346 , tolerance(tolerance)
1347 , owns_pointers(false)
1348 , finite_difference_step(0)
1356template <
int dim,
int spacedim,
int chartdim>
1361 const double tolerance)
1364 , push_forward_function(push_forward.
release())
1365 , pull_back_function(pull_back.
release())
1366 , tolerance(tolerance)
1367 , owns_pointers(true)
1368 , finite_difference_step(0)
1376template <
int dim,
int spacedim,
int chartdim>
1378 const std::string push_forward_expression,
1379 const std::string pull_back_expression,
1382 const std::string chart_vars,
1383 const std::string space_vars,
1384 const double tolerance,
1387 , const_map(const_map)
1388 , tolerance(tolerance)
1389 , owns_pointers(true)
1390 , push_forward_expression(push_forward_expression)
1391 , pull_back_expression(pull_back_expression)
1392 , chart_vars(chart_vars)
1393 , space_vars(space_vars)
1394 , finite_difference_step(h)
1406template <
int dim,
int spacedim,
int chartdim>
1409 if (owns_pointers ==
true)
1412 push_forward_function =
nullptr;
1416 pull_back_function =
nullptr;
1423template <
int dim,
int spacedim,
int chartdim>
1424std::unique_ptr<Manifold<dim, spacedim>>
1439 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1441 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1442 push_forward_expression,
1443 pull_back_expression,
1444 this->get_periodicity(),
1449 finite_difference_step);
1453 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1454 *push_forward_function,
1455 *pull_back_function,
1456 this->get_periodicity(),
1463template <
int dim,
int spacedim,
int chartdim>
1470 push_forward_function->vector_value(chart_point, pf);
1471 for (
unsigned int i = 0; i < spacedim; ++i)
1477 pull_back_function->vector_value(result, pb);
1478 for (
unsigned int i = 0; i < chartdim; ++i)
1480 (chart_point.
norm() > tolerance &&
1481 (
std::abs(pb[i] - chart_point[i]) <
1482 tolerance * chart_point.
norm())) ||
1483 (
std::abs(pb[i] - chart_point[i]) < tolerance),
1485 "The push forward is not the inverse of the pull back! Bailing out."));
1493template <
int dim,
int spacedim,
int chartdim>
1499 for (
unsigned int i = 0; i < spacedim; ++i)
1501 const auto gradient = push_forward_function->gradient(chart_point, i);
1502 for (
unsigned int j = 0; j < chartdim; ++j)
1503 DF[i][j] = gradient[j];
1510template <
int dim,
int spacedim,
int chartdim>
1517 pull_back_function->vector_value(space_point, pb);
1518 for (
unsigned int i = 0; i < chartdim; ++i)
1535 double phi = std::atan2(y, x);
1536 double theta = std::atan2(z,
std::sqrt(x * x + y * y) - centerline_radius);
1539 Utilities::fixed_power<2>(x -
std::cos(phi) * centerline_radius) +
1542 return {phi, theta, w};
1551 double phi = chart_point[0];
1552 double theta = chart_point[1];
1553 double w = chart_point[2];
1555 return {
std::cos(phi) * centerline_radius +
1557 inner_radius * w *
std::sin(theta),
1558 std::sin(phi) * centerline_radius +
1566 const double inner_radius)
1568 , centerline_radius(centerline_radius)
1569 , inner_radius(inner_radius)
1572 ExcMessage(
"The centerline radius must be greater than the "
1580std::unique_ptr<Manifold<dim, 3>>
1583 return std::make_unique<TorusManifold<dim>>(centerline_radius, inner_radius);
1594 double phi = chart_point[0];
1595 double theta = chart_point[1];
1596 double w = chart_point[2];
1598 DX[0][0] = -
std::sin(phi) * centerline_radius -
1604 DX[1][1] = inner_radius * w *
std::cos(theta);
1605 DX[1][2] = inner_radius *
std::sin(theta);
1607 DX[2][0] =
std::cos(phi) * centerline_radius +
1620template <
int dim,
int spacedim>
1631template <
int dim,
int spacedim>
1633 spacedim>::~TransfiniteInterpolationManifold()
1635 if (clear_signal.connected())
1636 clear_signal.disconnect();
1633 spacedim>::~TransfiniteInterpolationManifold() {
…}
1641template <
int dim,
int spacedim>
1642std::unique_ptr<Manifold<dim, spacedim>>
1648 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1653template <
int dim,
int spacedim>
1660 clear_signal.disconnect();
1662 this->triangulation =
nullptr;
1663 this->level_coarse = -1;
1666 coarse_cell_is_flat.resize(
triangulation.n_cells(level_coarse),
false);
1667 quadratic_approximation.clear();
1677 const std::vector<Point<dim>> unit_points =
1679 std::vector<Point<spacedim>> real_points(unit_points.size());
1681 for (
const auto &cell :
triangulation.active_cell_iterators())
1683 bool cell_is_flat =
true;
1684 for (
const auto l : cell->line_indices())
1685 if (cell->line(l)->manifold_id() != cell->manifold_id() &&
1687 cell_is_flat =
false;
1688 if constexpr (dim > 2)
1689 for (
const auto q : cell->face_indices())
1690 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1692 cell_is_flat =
false;
1694 coarse_cell_is_flat.size());
1695 coarse_cell_is_flat[cell->index()] = cell_is_flat;
1698 for (
unsigned int i = 0; i < unit_points.size(); ++i)
1699 real_points[i] = push_forward(cell, unit_points[i]);
1700 quadratic_approximation.emplace_back(real_points, unit_points);
1709 template <
typename AccessorType>
1711 compute_transfinite_interpolation(
const AccessorType &cell,
1715 return cell.vertex(0) * (1. - chart_point[0]) +
1716 cell.vertex(1) * chart_point[0];
1720 template <
typename AccessorType>
1722 compute_transfinite_interpolation(
const AccessorType &cell,
1724 const bool cell_is_flat)
1726 const unsigned int dim = AccessorType::dimension;
1727 const unsigned int spacedim = AccessorType::space_dimension;
1735 const std::array<Point<spacedim>, 4> vertices{
1736 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1741 std::array<double, 4> weights_vertices{
1742 {(1. - chart_point[0]) * (1. - chart_point[1]),
1743 chart_point[0] * (1. - chart_point[1]),
1744 (1. - chart_point[0]) * chart_point[1],
1745 chart_point[0] * chart_point[1]}};
1750 new_point += weights_vertices[v] * vertices[v];
1762 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1765 const auto weights_view =
1767 const auto points_view =
make_array_view(points.begin(), points.end());
1769 for (
unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1772 const double my_weight =
1773 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1774 const double line_point = chart_point[1 - line / 2];
1781 cell.line(line)->manifold_id();
1782 if (line_manifold_id == my_manifold_id ||
1787 my_weight * (1. - line_point);
1790 my_weight * line_point;
1798 weights[0] = 1. - line_point;
1799 weights[1] = line_point;
1808 new_point -= weights_vertices[v] * vertices[v];
1817 static constexpr unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1827 static constexpr unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1835 template <
typename AccessorType>
1837 compute_transfinite_interpolation(
const AccessorType &cell,
1839 const bool cell_is_flat)
1841 const unsigned int dim = AccessorType::dimension;
1842 const unsigned int spacedim = AccessorType::space_dimension;
1848 const std::array<Point<spacedim>, 8> vertices{{cell.vertex(0),
1860 double linear_shapes[10];
1861 for (
unsigned int d = 0;
d < 3; ++
d)
1863 linear_shapes[2 *
d] = 1. - chart_point[
d];
1864 linear_shapes[2 *
d + 1] = chart_point[
d];
1868 for (
unsigned int d = 6;
d < 10; ++
d)
1869 linear_shapes[d] = linear_shapes[d - 6];
1871 std::array<double, 8> weights_vertices;
1872 for (
unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1873 for (
unsigned int i1 = 0; i1 < 2; ++i1)
1874 for (
unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1875 weights_vertices[v] =
1876 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1880 for (
unsigned int v = 0; v < 8; ++v)
1881 new_point += weights_vertices[v] * vertices[v];
1887 std::array<double, GeometryInfo<3>::lines_per_cell> weights_lines;
1888 std::fill(weights_lines.begin(), weights_lines.end(), 0.0);
1891 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1894 const auto weights_view =
1896 const auto points_view =
make_array_view(points.begin(), points.end());
1898 for (
const unsigned int face :
GeometryInfo<3>::face_indices())
1900 const double my_weight = linear_shapes[face];
1901 const unsigned int face_even = face - face % 2;
1909 cell.face(face)->manifold_id();
1910 if (face_manifold_id == my_manifold_id ||
1913 for (
unsigned int line = 0;
1914 line < GeometryInfo<2>::lines_per_cell;
1917 const double line_weight =
1918 linear_shapes[face_even + 2 + line];
1919 weights_lines[face_to_cell_lines_3d[face][line]] +=
1920 my_weight * line_weight;
1926 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1927 linear_shapes[face_even + 2] *
1928 (linear_shapes[face_even + 4] * my_weight);
1929 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1930 linear_shapes[face_even + 3] *
1931 (linear_shapes[face_even + 4] * my_weight);
1932 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1933 linear_shapes[face_even + 2] *
1934 (linear_shapes[face_even + 5] * my_weight);
1935 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1936 linear_shapes[face_even + 3] *
1937 (linear_shapes[face_even + 5] * my_weight);
1942 points[v] = vertices[face_to_cell_vertices_3d[face][v]];
1944 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1946 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1948 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1950 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1958 const auto weights_view_line =
1960 const auto points_view_line =
1962 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1965 const double line_point =
1966 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1967 double my_weight = 0.;
1969 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1972 const unsigned int subline = line - 8;
1974 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1976 my_weight -= weights_lines[line];
1982 cell.line(line)->manifold_id();
1983 if (line_manifold_id == my_manifold_id ||
1988 my_weight * (1. - line_point);
1991 my_weight * (line_point);
1999 weights[0] = 1. - line_point;
2000 weights[1] = line_point;
2001 new_point -= my_weight * tria.
get_manifold(line_manifold_id)
2009 new_point += weights_vertices[v] * vertices[v];
2017template <
int dim,
int spacedim>
2029 ExcMessage(
"chart_point is not in unit interval"));
2031 return compute_transfinite_interpolation(*cell,
2033 coarse_cell_is_flat[cell->index()]);
2038template <
int dim,
int spacedim>
2047 for (
unsigned int d = 0; d < dim; ++d)
2050 const double step = chart_point[d] > 0.5 ? -1e-8 : 1e-8;
2053 modified[d] += step;
2055 compute_transfinite_interpolation(*cell,
2057 coarse_cell_is_flat[cell->index()]) -
2058 pushed_forward_chart_point;
2059 for (
unsigned int e = 0; e < spacedim; ++e)
2060 grad[e][d] = difference[e] / step;
2067template <
int dim,
int spacedim>
2075 for (
unsigned int d = 0; d < dim; ++d)
2079 Point<dim> chart_point = cell->reference_cell().closest_point(initial_guess);
2088 compute_transfinite_interpolation(*cell,
2090 coarse_cell_is_flat[cell->index()]);
2091 const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter());
2092 double residual_norm_square = residual.
norm_square();
2094 bool must_recompute_jacobian =
true;
2095 for (
unsigned int i = 0; i < 100; ++i)
2097 if (residual_norm_square < tolerance)
2104 for (
unsigned int d = 0; d < spacedim; ++d)
2105 for (
unsigned int e = 0; e < dim; ++e)
2106 update[e] += inv_grad[d][e] * residual[d];
2107 return chart_point + update;
2119 if (must_recompute_jacobian || i % 9 == 0)
2127 push_forward_gradient(cell,
2133 must_recompute_jacobian =
false;
2136 for (
unsigned int d = 0; d < spacedim; ++d)
2137 for (
unsigned int e = 0; e < dim; ++e)
2138 update[e] += inv_grad[d][e] * residual[d];
2152 while (alpha > 1e-4)
2154 Point<dim> guess = chart_point + alpha * update;
2156 point - compute_transfinite_interpolation(
2157 *cell, guess, coarse_cell_is_flat[cell->index()]);
2158 const double residual_norm_new = residual_guess.
norm_square();
2159 if (residual_norm_new < residual_norm_square)
2161 residual = residual_guess;
2162 residual_norm_square = residual_norm_new;
2163 chart_point += alpha * update;
2176 must_recompute_jacobian =
true;
2190 for (
unsigned int d = 0; d < spacedim; ++d)
2191 for (
unsigned int e = 0; e < dim; ++e)
2192 Jinv_deltaf[e] += inv_grad[d][e] * delta_f[d];
2199 if (
std::abs(delta_x * Jinv_deltaf) > 1e-12 && !must_recompute_jacobian)
2202 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
2204 for (
unsigned int d = 0; d < spacedim; ++d)
2205 for (
unsigned int e = 0; e < dim; ++e)
2206 jac_update[d] += delta_x[e] * inv_grad[d][e];
2207 for (
unsigned int d = 0; d < spacedim; ++d)
2208 for (
unsigned int e = 0; e < dim; ++e)
2209 inv_grad[d][e] += factor[e] * jac_update[d];
2217template <
int dim,
int spacedim>
2218std::array<unsigned int, 20>
2228 ExcMessage(
"The manifold was initialized with level " +
2229 std::to_string(level_coarse) +
" but there are now" +
2230 "active cells on a lower level. Coarsening the mesh is " +
2231 "currently not supported"));
2241 boost::container::small_vector<std::pair<double, unsigned int>, 200>
2242 distances_and_cells;
2243 for (; cell != endc; ++cell)
2246 if (&cell->get_manifold() !=
this)
2253 vertices[vertex_n] = cell->vertex(vertex_n);
2261 center += vertices[v];
2263 double radius_square = 0.;
2266 std::max(radius_square, (center - vertices[v]).norm_square());
2267 bool inside_circle =
true;
2268 for (
unsigned int i = 0; i < points.size(); ++i)
2269 if ((center - points[i]).norm_square() > radius_square * 1.5)
2271 inside_circle =
false;
2274 if (inside_circle ==
false)
2278 double current_distance = 0;
2279 for (
unsigned int i = 0; i < points.size(); ++i)
2282 quadratic_approximation[cell->index()].compute(points[i]);
2285 distances_and_cells.push_back(
2286 std::make_pair(current_distance, cell->index()));
2291 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2292 std::array<unsigned int, 20> cells;
2294 for (
unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2296 cells[i] = distances_and_cells[i].
second;
2303template <
int dim,
int spacedim>
2309 Assert(surrounding_points.size() == chart_points.size(),
2310 ExcMessage(
"The chart points array view must be as large as the "
2311 "surrounding points array view."));
2313 std::array<unsigned int, 20> nearby_cells =
2314 get_possible_cells_around_points(surrounding_points);
2338 auto guess_chart_point_structdim_2 = [&](
const unsigned int i) ->
Point<dim> {
2339 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2340 ExcMessage(
"This function assumes that there are eight surrounding "
2341 "points around a two-dimensional object. It also assumes "
2342 "that the first three chart points have already been "
2352 return chart_points[1] + (chart_points[2] - chart_points[0]);
2354 return 0.5 * (chart_points[0] + chart_points[2]);
2356 return 0.5 * (chart_points[1] + chart_points[3]);
2358 return 0.5 * (chart_points[0] + chart_points[1]);
2360 return 0.5 * (chart_points[2] + chart_points[3]);
2389 auto guess_chart_point_structdim_3 = [&](
const unsigned int i) ->
Point<dim> {
2390 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2391 ExcMessage(
"This function assumes that there are eight surrounding "
2392 "points around a three-dimensional object. It also "
2393 "assumes that the first five chart points have already "
2395 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2399 bool use_structdim_2_guesses =
false;
2400 bool use_structdim_3_guesses =
false;
2405 if (surrounding_points.size() == 8)
2408 surrounding_points[6] - surrounding_points[0];
2410 surrounding_points[7] - surrounding_points[2];
2413 const double cosine = scalar_product(v06, v27) /
2418 use_structdim_2_guesses =
true;
2419 else if (spacedim == 3)
2422 use_structdim_3_guesses =
true;
2425 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2426 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2431 auto compute_chart_point =
2433 const unsigned int point_index) {
2439 bool used_quadratic_approximation =
false;
2442 if (point_index == 3 && surrounding_points.size() >= 8)
2443 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2444 else if (use_structdim_2_guesses && 3 < point_index)
2445 guess = guess_chart_point_structdim_2(point_index);
2446 else if (use_structdim_3_guesses && 4 < point_index)
2447 guess = guess_chart_point_structdim_3(point_index);
2448 else if (dim == 3 && point_index > 7 && surrounding_points.size() == 26)
2450 if (point_index < 20)
2453 point_index - 8, 0)] +
2455 point_index - 8, 1)]);
2459 point_index - 20, 0)] +
2461 point_index - 20, 1)] +
2463 point_index - 20, 2)] +
2465 point_index - 20, 3)]);
2469 guess = quadratic_approximation[cell->index()].compute(
2470 surrounding_points[point_index]);
2471 used_quadratic_approximation =
true;
2473 chart_points[point_index] =
2474 pull_back(cell, surrounding_points[point_index], guess);
2479 if (chart_points[point_index][0] ==
2481 !used_quadratic_approximation)
2483 guess = quadratic_approximation[cell->index()].compute(
2484 surrounding_points[point_index]);
2485 chart_points[point_index] =
2486 pull_back(cell, surrounding_points[point_index], guess);
2489 if (chart_points[point_index][0] ==
2492 for (
unsigned int d = 0; d < dim; ++d)
2494 chart_points[point_index] =
2495 pull_back(cell, surrounding_points[point_index], guess);
2500 for (
unsigned int c = 0; c < nearby_cells.size(); ++c)
2504 bool inside_unit_cell =
true;
2505 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2507 compute_chart_point(cell, i);
2514 inside_unit_cell =
false;
2518 if (inside_unit_cell ==
true)
2525 if (c == nearby_cells.size() - 1 ||
2530 std::ostringstream message;
2531 for (
unsigned int b = 0; b <= c; ++b)
2535 message <<
"Looking at cell " << cell->id()
2536 <<
" with vertices: " << std::endl;
2538 message << cell->vertex(v) <<
" ";
2539 message << std::endl;
2540 message <<
"Transformation to chart coordinates: " << std::endl;
2541 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2543 compute_chart_point(cell, i);
2544 message << surrounding_points[i] <<
" -> " << chart_points[i]
2564template <
int dim,
int spacedim>
2570 boost::container::small_vector<Point<dim>, 100> chart_points(
2571 surrounding_points.size());
2574 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2577 chart_manifold.get_new_point(chart_points_view, weights);
2579 return push_forward(cell, p_chart);
2584template <
int dim,
int spacedim>
2594 boost::container::small_vector<Point<dim>, 100> chart_points(
2595 surrounding_points.size());
2598 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2600 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2602 chart_manifold.get_new_points(chart_points_view,
2605 new_points_on_chart.end()));
2607 for (
unsigned int row = 0; row < weights.size(0); ++row)
2608 new_points[row] = push_forward(cell, new_points_on_chart[row]);
2614#include "grid/manifold_lib.inst"
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
virtual Point< spacedim > push_forward(const Point< 3 > &chart_point) const override
virtual DerivativeForm< 1, 3, spacedim > push_forward_gradient(const Point< 3 > &chart_point) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
CylindricalManifold(const unsigned int axis=0, const double tolerance=1e-10)
virtual Point< 3 > pull_back(const Point< spacedim > &space_point) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > push_forward(const Point< spacedim > &chart_point) const override
EllipticalManifold(const Point< spacedim > ¢er, const Tensor< 1, spacedim > &major_axis_direction, const double eccentricity)
virtual DerivativeForm< 1, spacedim, spacedim > push_forward_gradient(const Point< spacedim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
static Tensor< 1, spacedim > get_periodicity()
const FunctionParser< spacedim >::ConstMap const_map
const std::string chart_vars
const std::string pull_back_expression
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const override
const std::string space_vars
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const override
FunctionManifold(const Function< chartdim > &push_forward_function, const Function< spacedim > &pull_back_function, const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >(), const double tolerance=1e-10)
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
const std::string push_forward_expression
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const override
ObserverPointer< const Function< chartdim >, FunctionManifold< dim, spacedim, chartdim > > push_forward_function
ObserverPointer< const Function< spacedim >, FunctionManifold< dim, spacedim, chartdim > > pull_back_function
virtual ~FunctionManifold() override
virtual void initialize(const std::string &vars, const std::vector< std::string > &expressions, const ConstMap &constants, const bool time_dependent=false) override
std::map< std::string, double > ConstMap
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const
Abstract base class for mapping classes.
PolarManifold(const Point< spacedim > center=Point< spacedim >())
static Tensor< 1, spacedim > get_periodicity()
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual DerivativeForm< 1, spacedim, spacedim > push_forward_gradient(const Point< spacedim > &chart_point) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
virtual Point< spacedim > push_forward(const Point< spacedim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &vertices, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
std::pair< double, Tensor< 1, spacedim > > guess_new_point(const ArrayView< const Tensor< 1, spacedim > > &directions, const ArrayView< const double > &distances, const ArrayView< const double > &weights) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
SphericalManifold(const Point< spacedim > center=Point< spacedim >())
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
void do_get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights, ArrayView< Point< spacedim > > new_points) const
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
virtual Point< 3 > pull_back(const Point< 3 > &p) const override
virtual DerivativeForm< 1, 3, 3 > push_forward_gradient(const Point< 3 > &chart_point) const override
virtual Point< 3 > push_forward(const Point< 3 > &chart_point) const override
TorusManifold(const double centerline_radius, const double inner_radius)
virtual std::unique_ptr< Manifold< dim, 3 > > clone() const override
DerivativeForm< 1, dim, spacedim > push_forward_gradient(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &chart_point, const Point< spacedim > &pushed_forward_chart_point) const
void initialize(const Triangulation< dim, spacedim > &triangulation)
Triangulation< dim, spacedim >::cell_iterator compute_chart_points(const ArrayView< const Point< spacedim > > &surrounding_points, ArrayView< Point< dim > > chart_points) const
std::array< unsigned int, 20 > get_possible_cells_around_points(const ArrayView< const Point< spacedim > > &surrounding_points) const
Point< dim > pull_back(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_guess) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Point< spacedim > push_forward(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &chart_point) const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Triangulation< dim, spacedim > & get_triangulation()
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number signed_angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b, const Tensor< 1, spacedim, Number > &axis)
Tensor< 1, 3 > projected_direction(const Tensor< 1, 3 > &u, const Tensor< 1, 3 > &v)
Tensor< 1, 3 > apply_exponential_map(const Tensor< 1, 3 > &u, const Tensor< 1, 3 > &dir)
static constexpr double invalid_pull_back_coordinate
Point< spacedim > compute_normal(const Tensor< 1, spacedim > &, bool=false)
constexpr unsigned int invalid_unsigned_int
constexpr types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > acos(const ::VectorizedArray< Number, width > &x)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static double distance_to_unit_cell(const Point< dim > &p)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
boost::signals2::signal< void()> clear
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)