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>
135template <
int dim,
int spacedim>
136std::unique_ptr<Manifold<dim, spacedim>>
139 return std::make_unique<PolarManifold<dim, spacedim>>(*this);
144template <
int dim,
int spacedim>
161template <
int dim,
int spacedim>
166 Assert(spherical_point[0] >= 0.0,
167 ExcMessage(
"Negative radius for given point."));
168 const double rho = spherical_point[0];
169 const double theta = spherical_point[1];
181 const double phi = spherical_point[2];
195template <
int dim,
int spacedim>
201 const double rho = R.
norm();
210 p[1] = std::atan2(R[1], R[0]);
218 const double z = R[2];
219 p[2] = std::atan2(R[1], R[0]);
222 p[1] = std::atan2(
std::sqrt(R[0] * R[0] + R[1] * R[1]), z);
234template <
int dim,
int spacedim>
239 Assert(spherical_point[0] >= 0.0,
240 ExcMessage(
"Negative radius for given point."));
241 const double rho = spherical_point[0];
242 const double theta = spherical_point[1];
259 const double phi = spherical_point[2];
284 template <
int dim,
int spacedim>
286 spherical_face_is_horizontal(
296 constexpr unsigned int n_vertices =
298 std::array<double, n_vertices> sqr_distances_to_center;
299 std::array<double, n_vertices - 1> sqr_distances_to_first_vertex;
300 sqr_distances_to_center[0] =
301 (face->vertex(0) - manifold_center).norm_square();
302 for (
unsigned int i = 1; i < n_vertices; ++i)
304 sqr_distances_to_center[i] =
305 (face->vertex(i) - manifold_center).norm_square();
306 sqr_distances_to_first_vertex[i - 1] =
307 (face->vertex(i) - face->vertex(0)).norm_square();
309 const auto minmax_sqr_distance =
310 std::minmax_element(sqr_distances_to_center.begin(),
311 sqr_distances_to_center.end());
312 const auto min_sqr_distance_to_first_vertex =
313 std::min_element(sqr_distances_to_first_vertex.begin(),
314 sqr_distances_to_first_vertex.end());
316 return (*minmax_sqr_distance.second - *minmax_sqr_distance.first <
317 1.e-10 * *min_sqr_distance_to_first_vertex);
323template <
int dim,
int spacedim>
333 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
340 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
341 return normalized_spherical_normal;
357template <
int dim,
int spacedim>
366template <
int dim,
int spacedim>
367std::unique_ptr<Manifold<dim, spacedim>>
370 return std::make_unique<SphericalManifold<dim, spacedim>>(*this);
375template <
int dim,
int spacedim>
380 const double w)
const
382 const double tol = 1e-10;
384 if ((p1 - p2).norm_square() < tol * tol ||
std::abs(w) < tol)
396 const double r1 =
v1.norm();
397 const double r2 = v2.
norm();
399 Assert(r1 > tol && r2 > tol,
400 ExcMessage(
"p1 and p2 cannot coincide with the center."));
406 const double cosgamma = e1 * e2;
409 if (cosgamma < -1 + 8. * std::numeric_limits<double>::epsilon())
413 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
419 const double sigma = w * std::acos(cosgamma);
424 const double n_norm = n.
norm();
427 "Probably, this means v1==v2 or v2==0."));
441template <
int dim,
int spacedim>
447 const double tol = 1e-10;
454 const double r1 =
v1.norm();
455 const double r2 = v2.
norm();
465 const double cosgamma = e1 * e2;
467 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
468 ExcMessage(
"p1 and p2 cannot lie on the same diameter and be opposite "
469 "respect to the center."));
471 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
477 const double n_norm = n.
norm();
480 "Probably, this means v1==v2 or v2==0."));
486 const double gamma = std::acos(cosgamma);
487 return (r2 - r1) * e1 + r1 * gamma * n;
492template <
int dim,
int spacedim>
502 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
509 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
510 return normalized_spherical_normal;
544template <
int dim,
int spacedim>
555 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
560 for (
unsigned int vertex = 0;
561 vertex < GeometryInfo<spacedim>::vertices_per_face;
563 face_vertex_normals[vertex] = face->vertex(vertex) -
center;
571template <
int dim,
int spacedim>
581 do_get_new_points(surrounding_points,
make_array_view(weights), new_points);
588template <
int dim,
int spacedim>
612 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] -
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)
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)
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 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];
1171template <
int dim,
int spacedim>
1175 const double eccentricity)
1178 , direction(major_axis_direction)
1180 , cosh_u(1.0 / eccentricity)
1181 , sinh_u(
std::sqrt(cosh_u * cosh_u - 1.0))
1188 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1190 Assert(direction_norm != 0.0,
1192 "Invalid major axis direction vector: Null vector not allowed."));
1198template <
int dim,
int spacedim>
1199std::unique_ptr<Manifold<dim, spacedim>>
1202 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1207template <
int dim,
int spacedim>
1220template <
int dim,
int spacedim>
1234 const double cs =
std::cos(chart_point[1]);
1235 const double sn =
std::sin(chart_point[1]);
1238 const double x = chart_point[0] * cosh_u * cs;
1239 const double y = chart_point[0] * sinh_u * sn;
1241 const Point<2> p(direction[0] * x - direction[1] * y,
1242 direction[1] * x + direction[0] * y);
1248template <
int dim,
int spacedim>
1263 const double x0 = space_point[0] -
center[0];
1264 const double y0 = space_point[1] -
center[1];
1265 const double x = direction[0] * x0 + direction[1] * y0;
1266 const double y = -direction[1] * x0 + direction[0] * y0;
1268 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1274 double cos_eta = x / (pt0 * cosh_u);
1283 const double eta = std::acos(cos_eta);
1284 const double pt1 = (std::signbit(y) ? 2.0 *
numbers::PI - eta : eta);
1290template <
int dim,
int spacedim>
1306 const double cs =
std::cos(chart_point[1]);
1307 const double sn =
std::sin(chart_point[1]);
1309 dX[0][0] = cosh_u * cs;
1310 dX[0][1] = -chart_point[0] * cosh_u * sn;
1311 dX[1][0] = sinh_u * sn;
1312 dX[1][1] = chart_point[0] * sinh_u * cs;
1316 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1326template <
int dim,
int spacedim,
int chartdim>
1331 const double tolerance)
1334 , push_forward_function(&push_forward_function)
1335 , pull_back_function(&pull_back_function)
1336 , tolerance(tolerance)
1337 , owns_pointers(false)
1338 , finite_difference_step(0)
1346template <
int dim,
int spacedim,
int chartdim>
1351 const double tolerance)
1354 , push_forward_function(push_forward.release())
1355 , pull_back_function(pull_back.release())
1356 , tolerance(tolerance)
1357 , owns_pointers(true)
1358 , finite_difference_step(0)
1366template <
int dim,
int spacedim,
int chartdim>
1368 const std::string push_forward_expression,
1369 const std::string pull_back_expression,
1372 const std::string chart_vars,
1373 const std::string space_vars,
1374 const double tolerance,
1377 , const_map(const_map)
1378 , tolerance(tolerance)
1379 , owns_pointers(true)
1380 , push_forward_expression(push_forward_expression)
1381 , pull_back_expression(pull_back_expression)
1382 , chart_vars(chart_vars)
1383 , space_vars(space_vars)
1384 , finite_difference_step(h)
1396template <
int dim,
int spacedim,
int chartdim>
1399 if (owns_pointers ==
true)
1402 push_forward_function =
nullptr;
1406 pull_back_function =
nullptr;
1413template <
int dim,
int spacedim,
int chartdim>
1414std::unique_ptr<Manifold<dim, spacedim>>
1429 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1431 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1432 push_forward_expression,
1433 pull_back_expression,
1434 this->get_periodicity(),
1439 finite_difference_step);
1443 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1444 *push_forward_function,
1445 *pull_back_function,
1446 this->get_periodicity(),
1453template <
int dim,
int spacedim,
int chartdim>
1460 push_forward_function->vector_value(chart_point, pf);
1461 for (
unsigned int i = 0; i < spacedim; ++i)
1466 pull_back_function->vector_value(result, pb);
1467 for (
unsigned int i = 0; i < chartdim; ++i)
1469 (chart_point.
norm() > tolerance &&
1470 (
std::abs(pb[i] - chart_point[i]) < tolerance * chart_point.
norm())) ||
1471 (
std::abs(pb[i] - chart_point[i]) < tolerance),
1473 "The push forward is not the inverse of the pull back! Bailing out."));
1481template <
int dim,
int spacedim,
int chartdim>
1487 for (
unsigned int i = 0; i < spacedim; ++i)
1489 const auto gradient = push_forward_function->gradient(chart_point, i);
1490 for (
unsigned int j = 0; j < chartdim; ++j)
1491 DF[i][j] = gradient[j];
1498template <
int dim,
int spacedim,
int chartdim>
1505 pull_back_function->vector_value(space_point, pb);
1506 for (
unsigned int i = 0; i < chartdim; ++i)
1523 double phi = std::atan2(y, x);
1524 double theta = std::atan2(z,
std::sqrt(x * x + y * y) - R);
1527 Utilities::fixed_power<2>(x -
std::cos(phi) * R) + z * z) /
1529 return {phi, theta, w};
1538 double phi = chart_point[0];
1539 double theta = chart_point[1];
1540 double w = chart_point[2];
1556 ExcMessage(
"Outer radius R must be greater than the inner "
1564std::unique_ptr<Manifold<dim, 3>>
1567 return std::make_unique<TorusManifold<dim>>(R, r);
1578 double phi = chart_point[0];
1579 double theta = chart_point[1];
1580 double w = chart_point[2];
1587 DX[1][1] = r * w *
std::cos(theta);
1602template <
int dim,
int spacedim>
1613template <
int dim,
int spacedim>
1615 spacedim>::~TransfiniteInterpolationManifold()
1617 if (clear_signal.connected())
1618 clear_signal.disconnect();
1623template <
int dim,
int spacedim>
1624std::unique_ptr<Manifold<dim, spacedim>>
1630 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1635template <
int dim,
int spacedim>
1642 clear_signal.disconnect();
1644 this->triangulation =
nullptr;
1645 this->level_coarse = -1;
1648 coarse_cell_is_flat.resize(
triangulation.n_cells(level_coarse),
false);
1649 quadratic_approximation.clear();
1659 std::vector<Point<dim>> unit_points =
1661 std::vector<Point<spacedim>> real_points(unit_points.size());
1663 for (
const auto &cell :
triangulation.active_cell_iterators())
1665 bool cell_is_flat =
true;
1666 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1667 if (cell->line(l)->manifold_id() != cell->manifold_id() &&
1669 cell_is_flat =
false;
1671 for (
unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1672 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1674 cell_is_flat =
false;
1676 coarse_cell_is_flat.size());
1677 coarse_cell_is_flat[cell->index()] = cell_is_flat;
1680 for (
unsigned int i = 0; i < unit_points.size(); ++i)
1681 real_points[i] = push_forward(cell, unit_points[i]);
1682 quadratic_approximation.emplace_back(real_points, unit_points);
1691 template <
typename AccessorType>
1693 compute_transfinite_interpolation(
const AccessorType &cell,
1697 return cell.vertex(0) * (1. - chart_point[0]) +
1698 cell.vertex(1) * chart_point[0];
1702 template <
typename AccessorType>
1704 compute_transfinite_interpolation(
const AccessorType &cell,
1706 const bool cell_is_flat)
1708 const unsigned int dim = AccessorType::dimension;
1709 const unsigned int spacedim = AccessorType::space_dimension;
1717 const std::array<Point<spacedim>, 4>
vertices{
1718 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1723 std::array<double, 4> weights_vertices{
1724 {(1. - chart_point[0]) * (1. - chart_point[1]),
1725 chart_point[0] * (1. - chart_point[1]),
1726 (1. - chart_point[0]) * chart_point[1],
1727 chart_point[0] * chart_point[1]}};
1732 new_point += weights_vertices[v] *
vertices[v];
1744 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1747 const auto weights_view =
1749 const auto points_view =
make_array_view(points.begin(), points.end());
1751 for (
unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1754 const double my_weight =
1755 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1756 const double line_point = chart_point[1 - line / 2];
1763 cell.line(line)->manifold_id();
1764 if (line_manifold_id == my_manifold_id ||
1769 my_weight * (1. - line_point);
1772 my_weight * line_point;
1780 weights[0] = 1. - line_point;
1781 weights[1] = line_point;
1790 new_point -= weights_vertices[v] *
vertices[v];
1799 static constexpr unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1809 static constexpr unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1817 template <
typename AccessorType>
1819 compute_transfinite_interpolation(
const AccessorType &cell,
1821 const bool cell_is_flat)
1823 const unsigned int dim = AccessorType::dimension;
1824 const unsigned int spacedim = AccessorType::space_dimension;
1830 const std::array<Point<spacedim>, 8>
vertices{{cell.vertex(0),
1842 double linear_shapes[10];
1843 for (
unsigned int d = 0;
d < 3; ++
d)
1845 linear_shapes[2 *
d] = 1. - chart_point[
d];
1846 linear_shapes[2 *
d + 1] = chart_point[
d];
1850 for (
unsigned int d = 6;
d < 10; ++
d)
1851 linear_shapes[d] = linear_shapes[d - 6];
1853 std::array<double, 8> weights_vertices;
1854 for (
unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1855 for (
unsigned int i1 = 0; i1 < 2; ++i1)
1856 for (
unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1857 weights_vertices[v] =
1858 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1862 for (
unsigned int v = 0; v < 8; ++v)
1863 new_point += weights_vertices[v] *
vertices[v];
1869 std::array<double, GeometryInfo<3>::lines_per_cell> weights_lines;
1870 std::fill(weights_lines.begin(), weights_lines.end(), 0.0);
1873 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1876 const auto weights_view =
1878 const auto points_view =
make_array_view(points.begin(), points.end());
1880 for (
const unsigned int face :
GeometryInfo<3>::face_indices())
1882 const double my_weight = linear_shapes[face];
1883 const unsigned int face_even = face - face % 2;
1891 cell.face(face)->manifold_id();
1892 if (face_manifold_id == my_manifold_id ||
1895 for (
unsigned int line = 0;
1896 line < GeometryInfo<2>::lines_per_cell;
1899 const double line_weight =
1900 linear_shapes[face_even + 2 + line];
1901 weights_lines[face_to_cell_lines_3d[face][line]] +=
1902 my_weight * line_weight;
1908 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1909 linear_shapes[face_even + 2] *
1910 (linear_shapes[face_even + 4] * my_weight);
1911 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1912 linear_shapes[face_even + 3] *
1913 (linear_shapes[face_even + 4] * my_weight);
1914 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1915 linear_shapes[face_even + 2] *
1916 (linear_shapes[face_even + 5] * my_weight);
1917 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1918 linear_shapes[face_even + 3] *
1919 (linear_shapes[face_even + 5] * my_weight);
1924 points[v] =
vertices[face_to_cell_vertices_3d[face][v]];
1926 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1928 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1930 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1932 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1940 const auto weights_view_line =
1942 const auto points_view_line =
1944 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1947 const double line_point =
1948 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1949 double my_weight = 0.;
1951 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1954 const unsigned int subline = line - 8;
1956 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1958 my_weight -= weights_lines[line];
1964 cell.line(line)->manifold_id();
1965 if (line_manifold_id == my_manifold_id ||
1970 my_weight * (1. - line_point);
1973 my_weight * (line_point);
1981 weights[0] = 1. - line_point;
1982 weights[1] = line_point;
1991 new_point += weights_vertices[v] *
vertices[v];
1999template <
int dim,
int spacedim>
2011 ExcMessage(
"chart_point is not in unit interval"));
2013 return compute_transfinite_interpolation(*cell,
2015 coarse_cell_is_flat[cell->index()]);
2020template <
int dim,
int spacedim>
2029 for (
unsigned int d = 0; d < dim; ++d)
2032 const double step = chart_point[d] > 0.5 ? -1e-8 : 1e-8;
2035 modified[d] += step;
2037 compute_transfinite_interpolation(*cell,
2039 coarse_cell_is_flat[cell->index()]) -
2040 pushed_forward_chart_point;
2041 for (
unsigned int e = 0; e < spacedim; ++e)
2042 grad[e][d] = difference[e] / step;
2049template <
int dim,
int spacedim>
2057 for (
unsigned int d = 0; d < dim; ++d)
2061 Point<dim> chart_point = cell->reference_cell().closest_point(initial_guess);
2070 compute_transfinite_interpolation(*cell,
2072 coarse_cell_is_flat[cell->index()]);
2073 const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter());
2074 double residual_norm_square = residual.
norm_square();
2076 bool must_recompute_jacobian =
true;
2077 for (
unsigned int i = 0; i < 100; ++i)
2079 if (residual_norm_square < tolerance)
2086 for (
unsigned int d = 0; d < spacedim; ++d)
2087 for (
unsigned int e = 0; e < dim; ++e)
2088 update[e] += inv_grad[d][e] * residual[d];
2089 return chart_point + update;
2101 if (must_recompute_jacobian || i % 9 == 0)
2109 push_forward_gradient(cell,
2115 must_recompute_jacobian =
false;
2118 for (
unsigned int d = 0; d < spacedim; ++d)
2119 for (
unsigned int e = 0; e < dim; ++e)
2120 update[e] += inv_grad[d][e] * residual[d];
2134 while (alpha > 1e-4)
2136 Point<dim> guess = chart_point + alpha * update;
2138 point - compute_transfinite_interpolation(
2139 *cell, guess, coarse_cell_is_flat[cell->index()]);
2140 const double residual_norm_new = residual_guess.
norm_square();
2141 if (residual_norm_new < residual_norm_square)
2143 residual = residual_guess;
2144 residual_norm_square = residual_norm_new;
2145 chart_point += alpha * update;
2158 must_recompute_jacobian =
true;
2172 for (
unsigned int d = 0; d < spacedim; ++d)
2173 for (
unsigned int e = 0; e < dim; ++e)
2174 Jinv_deltaf[e] += inv_grad[d][e] * delta_f[d];
2181 if (
std::abs(delta_x * Jinv_deltaf) > 1e-12 && !must_recompute_jacobian)
2184 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
2186 for (
unsigned int d = 0; d < spacedim; ++d)
2187 for (
unsigned int e = 0; e < dim; ++e)
2188 jac_update[d] += delta_x[e] * inv_grad[d][e];
2189 for (
unsigned int d = 0; d < spacedim; ++d)
2190 for (
unsigned int e = 0; e < dim; ++e)
2191 inv_grad[d][e] += factor[e] * jac_update[d];
2199template <
int dim,
int spacedim>
2200std::array<unsigned int, 20>
2210 ExcMessage(
"The manifold was initialized with level " +
2211 std::to_string(level_coarse) +
" but there are now" +
2212 "active cells on a lower level. Coarsening the mesh is " +
2213 "currently not supported"));
2223 boost::container::small_vector<std::pair<double, unsigned int>, 200>
2224 distances_and_cells;
2225 for (; cell != endc; ++cell)
2228 if (&cell->get_manifold() !=
this)
2235 vertices[vertex_n] = cell->vertex(vertex_n);
2245 double radius_square = 0.;
2249 bool inside_circle =
true;
2250 for (
unsigned int i = 0; i < points.size(); ++i)
2251 if ((
center - points[i]).norm_square() > radius_square * 1.5)
2253 inside_circle =
false;
2256 if (inside_circle ==
false)
2260 double current_distance = 0;
2261 for (
unsigned int i = 0; i < points.size(); ++i)
2264 quadratic_approximation[cell->index()].compute(points[i]);
2267 distances_and_cells.push_back(
2268 std::make_pair(current_distance, cell->index()));
2273 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2274 std::array<unsigned int, 20> cells;
2276 for (
unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2278 cells[i] = distances_and_cells[i].
second;
2285template <
int dim,
int spacedim>
2291 Assert(surrounding_points.size() == chart_points.size(),
2292 ExcMessage(
"The chart points array view must be as large as the "
2293 "surrounding points array view."));
2295 std::array<unsigned int, 20> nearby_cells =
2296 get_possible_cells_around_points(surrounding_points);
2320 auto guess_chart_point_structdim_2 = [&](
const unsigned int i) ->
Point<dim> {
2321 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2322 ExcMessage(
"This function assumes that there are eight surrounding "
2323 "points around a two-dimensional object. It also assumes "
2324 "that the first three chart points have already been "
2334 return chart_points[1] + (chart_points[2] - chart_points[0]);
2336 return 0.5 * (chart_points[0] + chart_points[2]);
2338 return 0.5 * (chart_points[1] + chart_points[3]);
2340 return 0.5 * (chart_points[0] + chart_points[1]);
2342 return 0.5 * (chart_points[2] + chart_points[3]);
2371 auto guess_chart_point_structdim_3 = [&](
const unsigned int i) ->
Point<dim> {
2372 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2373 ExcMessage(
"This function assumes that there are eight surrounding "
2374 "points around a three-dimensional object. It also "
2375 "assumes that the first five chart points have already "
2377 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2381 bool use_structdim_2_guesses =
false;
2382 bool use_structdim_3_guesses =
false;
2387 if (surrounding_points.size() == 8)
2390 surrounding_points[6] - surrounding_points[0];
2392 surrounding_points[7] - surrounding_points[2];
2395 const double cosine = scalar_product(v06, v27) /
2400 use_structdim_2_guesses =
true;
2401 else if (spacedim == 3)
2404 use_structdim_3_guesses =
true;
2407 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2408 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2413 auto compute_chart_point =
2415 const unsigned int point_index) {
2421 bool used_quadratic_approximation =
false;
2424 if (point_index == 3 && surrounding_points.size() >= 8)
2425 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2426 else if (use_structdim_2_guesses && 3 < point_index)
2427 guess = guess_chart_point_structdim_2(point_index);
2428 else if (use_structdim_3_guesses && 4 < point_index)
2429 guess = guess_chart_point_structdim_3(point_index);
2430 else if (dim == 3 && point_index > 7 && surrounding_points.size() == 26)
2432 if (point_index < 20)
2435 point_index - 8, 0)] +
2437 point_index - 8, 1)]);
2441 point_index - 20, 0)] +
2443 point_index - 20, 1)] +
2445 point_index - 20, 2)] +
2447 point_index - 20, 3)]);
2451 guess = quadratic_approximation[cell->index()].compute(
2452 surrounding_points[point_index]);
2453 used_quadratic_approximation =
true;
2455 chart_points[point_index] =
2456 pull_back(cell, surrounding_points[point_index], guess);
2461 if (chart_points[point_index][0] ==
2463 !used_quadratic_approximation)
2465 guess = quadratic_approximation[cell->index()].compute(
2466 surrounding_points[point_index]);
2467 chart_points[point_index] =
2468 pull_back(cell, surrounding_points[point_index], guess);
2471 if (chart_points[point_index][0] ==
2474 for (
unsigned int d = 0; d < dim; ++d)
2476 chart_points[point_index] =
2477 pull_back(cell, surrounding_points[point_index], guess);
2482 for (
unsigned int c = 0; c < nearby_cells.size(); ++c)
2486 bool inside_unit_cell =
true;
2487 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2489 compute_chart_point(cell, i);
2496 inside_unit_cell =
false;
2500 if (inside_unit_cell ==
true)
2507 if (c == nearby_cells.size() - 1 ||
2512 std::ostringstream message;
2513 for (
unsigned int b = 0; b <= c; ++b)
2517 message <<
"Looking at cell " << cell->id()
2518 <<
" with vertices: " << std::endl;
2520 message << cell->vertex(v) <<
" ";
2521 message << std::endl;
2522 message <<
"Transformation to chart coordinates: " << std::endl;
2523 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2525 compute_chart_point(cell, i);
2526 message << surrounding_points[i] <<
" -> " << chart_points[i]
2546template <
int dim,
int spacedim>
2552 boost::container::small_vector<Point<dim>, 100> chart_points(
2553 surrounding_points.size());
2556 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2559 chart_manifold.get_new_point(chart_points_view, weights);
2561 return push_forward(cell, p_chart);
2566template <
int dim,
int spacedim>
2576 boost::container::small_vector<Point<dim>, 100> chart_points(
2577 surrounding_points.size());
2580 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2582 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2584 chart_manifold.get_new_points(chart_points_view,
2587 new_points_on_chart.end()));
2589 for (
unsigned int row = 0; row < weights.size(0); ++row)
2590 new_points[row] = push_forward(cell, new_points_on_chart[row]);
2596#include "manifold_lib.inst"
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, 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
Tensor< 1, spacedim > direction
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()
SmartPointer< const Function< spacedim >, FunctionManifold< dim, spacedim, chartdim > > pull_back_function
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
SmartPointer< const Function< chartdim >, FunctionManifold< dim, spacedim, chartdim > > push_forward_function
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
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
TorusManifold(const double R, const double r)
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
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
#define DEAL_II_NAMESPACE_CLOSE
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
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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)
static constexpr double PI
static const unsigned int invalid_unsigned_int
const 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 > &)
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 > &)