32 #include <boost/container/small_vector.hpp>
56 const double theta = dir.
norm();
65 std::cos(theta) * u + std::sin(theta) * dirUnit;
66 return tmp / tmp.
norm();
83 template <
int spacedim>
95 ExcMessage(
"The direction parameter must not be zero!"));
97 if (std::abs(vector[0]) >= std::abs(vector[1]) &&
98 std::abs(vector[0]) >= std::abs(vector[2]))
102 normal[0] = (vector[1] + vector[2]) / vector[0];
104 else if (std::abs(vector[1]) >= std::abs(vector[0]) &&
105 std::abs(vector[1]) >= std::abs(vector[2]))
109 normal[1] = (vector[0] + vector[2]) / vector[1];
115 normal[2] = (vector[0] + vector[1]) / vector[2];
118 normal /= normal.
norm();
129 template <
int dim,
int spacedim>
138 template <
int dim,
int spacedim>
139 std::unique_ptr<Manifold<dim, spacedim>>
142 return std::make_unique<PolarManifold<dim, spacedim>>(*this);
147 template <
int dim,
int spacedim>
164 template <
int dim,
int spacedim>
169 Assert(spherical_point[0] >= 0.0,
170 ExcMessage(
"Negative radius for given point."));
171 const double rho = spherical_point[0];
172 const double theta = spherical_point[1];
179 p[0] = rho * std::cos(theta);
180 p[1] = rho * std::sin(theta);
184 const double phi = spherical_point[2];
185 p[0] = rho * std::sin(theta) * std::cos(phi);
186 p[1] = rho * std::sin(theta) * std::sin(phi);
187 p[2] = rho * std::cos(theta);
198 template <
int dim,
int spacedim>
204 const double rho = R.
norm();
221 const double z = R[2];
225 p[1] =
std::atan2(std::sqrt(R[0] * R[0] + R[1] * R[1]), z);
237 template <
int dim,
int spacedim>
242 Assert(spherical_point[0] >= 0.0,
243 ExcMessage(
"Negative radius for given point."));
244 const double rho = spherical_point[0];
245 const double theta = spherical_point[1];
253 DX[0][0] = std::cos(theta);
254 DX[0][1] = -rho * std::sin(theta);
255 DX[1][0] = std::sin(theta);
256 DX[1][1] = rho * std::cos(theta);
262 const double phi = spherical_point[2];
263 DX[0][0] = std::sin(theta) * std::cos(phi);
264 DX[0][1] = rho * std::cos(theta) * std::cos(phi);
265 DX[0][2] = -rho * std::sin(theta) * std::sin(phi);
267 DX[1][0] = std::sin(theta) * std::sin(phi);
268 DX[1][1] = rho * std::cos(theta) * std::sin(phi);
269 DX[1][2] = rho * std::sin(theta) * std::cos(phi);
271 DX[2][0] = std::cos(theta);
272 DX[2][1] = -rho * std::sin(theta);
287 template <
int dim,
int spacedim>
289 spherical_face_is_horizontal(
299 constexpr
unsigned int n_vertices =
301 std::array<double, n_vertices> sqr_distances_to_center;
302 std::array<
double, n_vertices - 1> sqr_distances_to_first_vertex;
303 sqr_distances_to_center[0] =
304 (face->vertex(0) - manifold_center).norm_square();
305 for (
unsigned int i = 1; i < n_vertices; ++i)
307 sqr_distances_to_center[i] =
308 (face->vertex(i) - manifold_center).norm_square();
309 sqr_distances_to_first_vertex[i - 1] =
310 (face->vertex(i) - face->vertex(0)).norm_square();
312 const auto minmax_sqr_distance =
313 std::minmax_element(sqr_distances_to_center.begin(),
314 sqr_distances_to_center.end());
315 const auto min_sqr_distance_to_first_vertex =
316 std::min_element(sqr_distances_to_first_vertex.begin(),
317 sqr_distances_to_first_vertex.end());
319 return (*minmax_sqr_distance.second - *minmax_sqr_distance.first <
320 1.e-10 * *min_sqr_distance_to_first_vertex);
326 template <
int dim,
int spacedim>
336 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
343 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
344 return normalized_spherical_normal;
360 template <
int dim,
int spacedim>
369 template <
int dim,
int spacedim>
370 std::unique_ptr<Manifold<dim, spacedim>>
373 return std::make_unique<SphericalManifold<dim, spacedim>>(*this);
378 template <
int dim,
int spacedim>
383 const double w)
const
385 const double tol = 1
e-10;
387 if ((p1 - p2).norm_square() < tol * tol || std::abs(
w) < tol)
389 else if (std::abs(
w - 1.0) < tol)
399 const double r1 =
v1.norm();
400 const double r2 = v2.
norm();
402 Assert(r1 > tol && r2 > tol,
403 ExcMessage(
"p1 and p2 cannot coincide with the center."));
409 const double cosgamma = e1 * e2;
427 const double n_norm = n.
norm();
430 "Probably, this means v1==v2 or v2==0."));
444 template <
int dim,
int spacedim>
450 const double tol = 1
e-10;
457 const double r1 =
v1.norm();
458 const double r2 = v2.
norm();
468 const double cosgamma = e1 * e2;
471 ExcMessage(
"p1 and p2 cannot lie on the same diameter and be opposite "
472 "respect to the center."));
480 const double n_norm = n.
norm();
483 "Probably, this means v1==v2 or v2==0."));
490 return (r2 - r1) * e1 + r1 *
gamma * n;
495 template <
int dim,
int spacedim>
505 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
512 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
513 return normalized_spherical_normal;
547 template <
int dim,
int spacedim>
558 if (spherical_face_is_horizontal<dim, spacedim>(face,
center))
563 for (
unsigned int vertex = 0;
564 vertex < GeometryInfo<spacedim>::vertices_per_face;
566 face_vertex_normals[vertex] = face->vertex(vertex) -
center;
574 template <
int dim,
int spacedim>
584 get_new_points(surrounding_points,
make_array_view(weights), new_points);
591 template <
int dim,
int spacedim>
609 template <
int dim,
int spacedim>
617 new_points.size() * surrounding_points.size());
618 const unsigned int weight_rows = new_points.size();
619 const unsigned int weight_columns = surrounding_points.size();
621 if (surrounding_points.size() == 2)
623 for (
unsigned int row = 0; row < weight_rows; ++row)
625 get_intermediate_point(surrounding_points[0],
626 surrounding_points[1],
627 weights[row * weight_columns + 1]);
631 boost::container::small_vector<std::pair<double, Tensor<1, spacedim>>, 100>
632 new_candidates(new_points.size());
633 boost::container::small_vector<Tensor<1, spacedim>, 100> directions(
635 boost::container::small_vector<double, 100> distances(
636 surrounding_points.size(), 0.0);
637 double max_distance = 0.;
638 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
640 directions[i] = surrounding_points[i] -
center;
641 distances[i] = directions[i].
norm();
643 if (distances[i] != 0.)
644 directions[i] /= distances[i];
647 ExcMessage(
"One of the vertices coincides with the center. "
648 "This is not allowed!"));
652 for (
unsigned int k = 0; k < i; ++k)
654 const double squared_distance =
655 (directions[i] - directions[k]).norm_square();
656 max_distance =
std::max(max_distance, squared_distance);
662 const double tolerance = 1
e-10;
663 boost::container::small_vector<bool, 100> accurate_point_was_found(
664 new_points.size(),
false);
669 for (
unsigned int row = 0; row < weight_rows; ++row)
671 new_candidates[row] =
672 guess_new_point(array_directions,
679 if (new_candidates[row].
first == 0.0)
682 accurate_point_was_found[row] =
true;
689 new_points[row] = polar_manifold.get_new_point(
704 if (max_distance < 2
e-2)
706 for (
unsigned int row = 0; row < weight_rows; ++row)
718 boost::container::small_vector<double, 1000> merged_weights(weights.
size());
719 boost::container::small_vector<Tensor<1, spacedim>, 100> merged_directions(
721 boost::container::small_vector<double, 100> merged_distances(
722 surrounding_points.size(), 0.0);
724 unsigned int n_unique_directions = 0;
725 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
727 bool found_duplicate =
false;
731 for (
unsigned int j = 0; j < n_unique_directions; ++j)
733 const double squared_distance =
734 (directions[i] - directions[j]).norm_square();
735 if (!found_duplicate && squared_distance < 1
e-28)
737 found_duplicate =
true;
738 for (
unsigned int row = 0; row < weight_rows; ++row)
739 merged_weights[row * weight_columns + j] +=
740 weights[row * weight_columns + i];
744 if (found_duplicate ==
false)
746 merged_directions[n_unique_directions] = directions[i];
747 merged_distances[n_unique_directions] = distances[i];
748 for (
unsigned int row = 0; row < weight_rows; ++row)
749 merged_weights[row * weight_columns + n_unique_directions] =
750 weights[row * weight_columns + i];
752 ++n_unique_directions;
758 boost::container::small_vector<unsigned int, 100> merged_weights_index(
760 for (
unsigned int row = 0; row < weight_rows; ++row)
762 for (
unsigned int existing_row = 0; existing_row < row; ++existing_row)
764 bool identical_weights =
true;
766 for (
unsigned int weight_index = 0;
767 weight_index < n_unique_directions;
769 if (std::abs(merged_weights[row * weight_columns + weight_index] -
770 merged_weights[existing_row * weight_columns +
771 weight_index]) > tolerance)
773 identical_weights =
false;
777 if (identical_weights)
779 merged_weights_index[row] = existing_row;
789 merged_directions.begin() + n_unique_directions);
792 merged_distances.begin() + n_unique_directions);
794 for (
unsigned int row = 0; row < weight_rows; ++row)
795 if (!accurate_point_was_found[row])
800 &merged_weights[row * weight_columns], n_unique_directions);
801 new_candidates[row].second =
802 get_new_point(array_merged_directions,
803 array_merged_distances,
804 array_merged_weights,
808 new_candidates[row].second =
809 new_candidates[merged_weights_index[row]].second;
812 center + new_candidates[row].first * new_candidates[row].second;
818 template <
int dim,
int spacedim>
819 std::pair<double, Tensor<1, spacedim>>
825 const double tolerance = 1
e-10;
830 double total_weights = 0.;
831 for (
unsigned int i = 0; i < directions.size(); ++i)
834 if (std::abs(1 - weights[i]) < tolerance)
835 return std::make_pair(distances[i], directions[i]);
837 rho += distances[i] * weights[i];
838 candidate += directions[i] * weights[i];
839 total_weights += weights[i];
843 const double norm = candidate.
norm();
847 rho /= total_weights;
849 return std::make_pair(rho, candidate);
855 template <
int spacedim>
878 Point<3> candidate = candidate_point;
879 const unsigned int n_merged_points = directions.size();
880 const double tolerance = 1
e-10;
881 const int max_iterations = 10;
886 for (
unsigned int i = 0; i < n_merged_points; ++i)
888 const double squared_distance =
889 (candidate - directions[i]).norm_square();
890 if (squared_distance < tolerance * tolerance)
896 if (n_merged_points == 2)
899 Assert(std::abs(weights[0] + weights[1] - 1.0) < 1
e-13,
915 for (
unsigned int i = 0; i < max_iterations; ++i)
921 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
929 for (
unsigned int i = 0; i < n_merged_points; ++i)
930 if (std::abs(weights[i]) > 1.e-15)
934 const double sintheta =
std::sqrt(sinthetaSq);
935 if (sintheta < tolerance)
937 Hessian[0][0] += weights[i];
938 Hessian[1][1] += weights[i];
942 const double costheta = (directions[i]) * candidate;
943 const double theta =
std::atan2(sintheta, costheta);
944 const double sincthetaInv = theta / sintheta;
946 const double cosphi = vPerp * Clocalx;
947 const double sinphi = vPerp * Clocaly;
949 gradlocal[0] = cosphi;
950 gradlocal[1] = sinphi;
951 gradient += (weights[i] * sincthetaInv) * gradlocal;
953 const double wt = weights[i] / sinthetaSq;
954 const double sinphiSq = sinphi * sinphi;
955 const double cosphiSq = cosphi * cosphi;
956 const double tt = sincthetaInv * costheta;
957 const double offdiag = cosphi * sinphi * wt * (1.0 - tt);
958 Hessian[0][0] += wt * (cosphiSq + tt * sinphiSq);
959 Hessian[0][1] += offdiag;
960 Hessian[1][0] += offdiag;
961 Hessian[1][1] += wt * (sinphiSq + tt * cosphiSq);
971 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
974 const Point<3> candidateOld = candidate;
979 if ((candidate - candidateOld).norm_square() < tolerance * tolerance)
989 template <
int dim,
int spacedim>
1009 const Point<3> & candidate_point)
const
1011 return do_get_new_point(directions, distances, weights, candidate_point);
1022 const Point<3> & candidate_point)
const
1024 return do_get_new_point(directions, distances, weights, candidate_point);
1035 const Point<3> & candidate_point)
const
1037 return do_get_new_point(directions, distances, weights, candidate_point);
1045 template <
int dim,
int spacedim>
1047 const double tolerance)
1055 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1060 template <
int dim,
int spacedim>
1064 const double tolerance)
1067 , direction(direction / direction.norm())
1068 , point_on_axis(point_on_axis)
1069 , tolerance(tolerance)
1070 , dxn(cross_product_3d(this->direction, normal_direction))
1075 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1080 template <
int dim,
int spacedim>
1081 std::unique_ptr<Manifold<dim, spacedim>>
1084 return std::make_unique<CylindricalManifold<dim, spacedim>>(*this);
1089 template <
int dim,
int spacedim>
1096 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1100 double average_length = 0.;
1101 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
1103 middle += surrounding_points[i] * weights[i];
1104 average_length += surrounding_points[i].square() * weights[i];
1106 middle -= point_on_axis;
1107 const double lambda = middle * direction;
1109 if ((middle - direction * lambda).square() < tolerance * average_length)
1110 return point_on_axis + direction * lambda;
1118 template <
int dim,
int spacedim>
1124 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1128 const double lambda = normalized_point * direction;
1129 const Point<spacedim> projection = point_on_axis + direction * lambda;
1144 template <
int dim,
int spacedim>
1150 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1153 const double sine_r = std::sin(chart_point(1)) * chart_point(0);
1154 const double cosine_r = std::cos(chart_point(1)) * chart_point(0);
1156 normal_direction * cosine_r + dxn * sine_r;
1159 return point_on_axis + direction * chart_point(2) + intermediate;
1164 template <
int dim,
int spacedim>
1170 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1175 const double sine = std::sin(chart_point(1));
1176 const double cosine = std::cos(chart_point(1));
1178 normal_direction * cosine + dxn * sine;
1181 constexpr
int s0 = 0 % spacedim;
1182 constexpr
int s1 = 1 % spacedim;
1183 constexpr
int s2 = 2 % spacedim;
1186 derivatives[s0][s0] = intermediate[s0];
1187 derivatives[s1][s0] = intermediate[s1];
1188 derivatives[s2][s0] = intermediate[s2];
1191 derivatives[s0][s1] = -normal_direction[s0] * sine + dxn[s0] * cosine;
1192 derivatives[s1][s1] = -normal_direction[s1] * sine + dxn[s1] * cosine;
1193 derivatives[s2][s1] = -normal_direction[s2] * sine + dxn[s2] * cosine;
1196 derivatives[s0][s2] = direction[s0];
1197 derivatives[s1][s2] = direction[s1];
1198 derivatives[s2][s2] = direction[s2];
1208 template <
int dim,
int spacedim>
1212 const double eccentricity)
1215 , direction(major_axis_direction)
1217 , cosh_u(1.0 / eccentricity)
1218 , sinh_u(std::
sqrt(cosh_u * cosh_u - 1.0))
1225 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1227 Assert(direction_norm != 0.0,
1229 "Invalid major axis direction vector: Null vector not allowed."));
1235 template <
int dim,
int spacedim>
1236 std::unique_ptr<Manifold<dim, spacedim>>
1239 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1244 template <
int dim,
int spacedim>
1257 template <
int dim,
int spacedim>
1271 const double cs = std::cos(chart_point[1]);
1272 const double sn = std::sin(chart_point[1]);
1275 const double x = chart_point[0] * cosh_u * cs;
1276 const double y = chart_point[0] * sinh_u * sn;
1278 const Point<2> p(direction[0] * x - direction[1] * y,
1279 direction[1] * x + direction[0] * y);
1285 template <
int dim,
int spacedim>
1300 const double x0 = space_point[0] -
center[0];
1301 const double y0 = space_point[1] -
center[1];
1302 const double x = direction[0] * x0 + direction[1] * y0;
1303 const double y = -direction[1] * x0 + direction[0] * y0;
1305 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1311 double cos_eta = x / (pt0 * cosh_u);
1321 const double pt1 = (std::signbit(y) ? 2.0 *
numbers::PI - eta : eta);
1327 template <
int dim,
int spacedim>
1343 const double cs = std::cos(chart_point[1]);
1344 const double sn = std::sin(chart_point[1]);
1346 dX[0][0] = cosh_u * cs;
1347 dX[0][1] = -chart_point[0] * cosh_u * sn;
1348 dX[1][0] = sinh_u * sn;
1349 dX[1][1] = chart_point[0] * sinh_u * cs;
1353 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1363 template <
int dim,
int spacedim,
int chartdim>
1368 const double tolerance)
1371 , push_forward_function(&push_forward_function)
1372 , pull_back_function(&pull_back_function)
1373 , tolerance(tolerance)
1374 , owns_pointers(false)
1375 , finite_difference_step(0)
1383 template <
int dim,
int spacedim,
int chartdim>
1388 const double tolerance)
1392 , pull_back_function(
pull_back.release())
1393 , tolerance(tolerance)
1394 , owns_pointers(true)
1395 , finite_difference_step(0)
1403 template <
int dim,
int spacedim,
int chartdim>
1405 const std::string push_forward_expression,
1406 const std::string pull_back_expression,
1409 const std::string chart_vars,
1410 const std::string space_vars,
1411 const double tolerance,
1414 , const_map(const_map)
1415 , tolerance(tolerance)
1416 , owns_pointers(true)
1417 , push_forward_expression(push_forward_expression)
1418 , pull_back_expression(pull_back_expression)
1419 , chart_vars(chart_vars)
1420 , space_vars(space_vars)
1421 , finite_difference_step(h)
1433 template <
int dim,
int spacedim,
int chartdim>
1436 if (owns_pointers ==
true)
1439 push_forward_function =
nullptr;
1443 pull_back_function =
nullptr;
1450 template <
int dim,
int spacedim,
int chartdim>
1451 std::unique_ptr<Manifold<dim, spacedim>>
1466 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1468 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1469 push_forward_expression,
1470 pull_back_expression,
1471 this->get_periodicity(),
1476 finite_difference_step);
1480 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1481 *push_forward_function,
1482 *pull_back_function,
1483 this->get_periodicity(),
1490 template <
int dim,
int spacedim,
int chartdim>
1497 push_forward_function->vector_value(chart_point, pf);
1498 for (
unsigned int i = 0; i < spacedim; ++i)
1503 pull_back_function->vector_value(result, pb);
1504 for (
unsigned int i = 0; i < chartdim; ++i)
1506 (chart_point.
norm() > tolerance &&
1507 (std::abs(pb[i] - chart_point[i]) < tolerance * chart_point.
norm())) ||
1508 (std::abs(pb[i] - chart_point[i]) < tolerance),
1510 "The push forward is not the inverse of the pull back! Bailing out."));
1518 template <
int dim,
int spacedim,
int chartdim>
1524 for (
unsigned int i = 0; i < spacedim; ++i)
1526 const auto gradient = push_forward_function->gradient(chart_point, i);
1527 for (
unsigned int j = 0; j < chartdim; ++j)
1528 DF[i][j] = gradient[j];
1535 template <
int dim,
int spacedim,
int chartdim>
1542 pull_back_function->vector_value(space_point, pb);
1543 for (
unsigned int i = 0; i < chartdim; ++i)
1561 double theta =
std::atan2(z, std::sqrt(x * x + y * y) - R);
1562 double w = std::sqrt(std::pow(y - std::sin(phi) * R, 2.0) +
1563 std::pow(x - std::cos(phi) * R, 2.0) + z * z) /
1565 return {phi, theta,
w};
1574 double phi = chart_point(0);
1575 double theta = chart_point(1);
1576 double w = chart_point(2);
1578 return {std::cos(phi) * R + r *
w * std::cos(theta) * std::cos(phi),
1579 r *
w * std::sin(theta),
1580 std::sin(phi) * R + r *
w * std::cos(theta) * std::sin(phi)};
1592 ExcMessage(
"Outer radius R must be greater than the inner "
1600 std::unique_ptr<Manifold<dim, 3>>
1603 return std::make_unique<TorusManifold<dim>>(R, r);
1614 double phi = chart_point(0);
1615 double theta = chart_point(1);
1616 double w = chart_point(2);
1618 DX[0][0] = -std::sin(phi) * R - r *
w * std::cos(theta) * std::sin(phi);
1619 DX[0][1] = -r *
w * std::sin(theta) * std::cos(phi);
1620 DX[0][2] = r * std::cos(theta) * std::cos(phi);
1623 DX[1][1] = r *
w * std::cos(theta);
1624 DX[1][2] = r * std::sin(theta);
1626 DX[2][0] = std::cos(phi) * R + r *
w * std::cos(theta) * std::cos(phi);
1627 DX[2][1] = -r *
w * std::sin(theta) * std::sin(phi);
1628 DX[2][2] = r * std::cos(theta) * std::sin(phi);
1638 template <
int dim,
int spacedim>
1649 template <
int dim,
int spacedim>
1651 spacedim>::~TransfiniteInterpolationManifold()
1653 if (clear_signal.connected())
1654 clear_signal.disconnect();
1659 template <
int dim,
int spacedim>
1660 std::unique_ptr<Manifold<dim, spacedim>>
1666 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1671 template <
int dim,
int spacedim>
1678 clear_signal.disconnect();
1679 clear_signal =
triangulation.signals.clear.connect([&]() ->
void {
1680 this->triangulation =
nullptr;
1681 this->level_coarse = -1;
1684 coarse_cell_is_flat.resize(
triangulation.n_cells(level_coarse),
false);
1685 quadratic_approximation.clear();
1695 std::vector<Point<dim>> unit_points =
1697 std::vector<Point<spacedim>> real_points(unit_points.size());
1699 for (
const auto &cell :
triangulation.active_cell_iterators())
1701 bool cell_is_flat =
true;
1702 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++
l)
1703 if (cell->line(
l)->manifold_id() != cell->manifold_id() &&
1705 cell_is_flat =
false;
1707 for (
unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1708 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1710 cell_is_flat =
false;
1712 coarse_cell_is_flat.size());
1713 coarse_cell_is_flat[cell->index()] = cell_is_flat;
1716 for (
unsigned int i = 0; i < unit_points.size(); ++i)
1718 quadratic_approximation.emplace_back(real_points, unit_points);
1727 template <
typename AccessorType>
1729 compute_transfinite_interpolation(
const AccessorType &cell,
1733 return cell.vertex(0) * (1. - chart_point[0]) +
1734 cell.vertex(1) * chart_point[0];
1738 template <
typename AccessorType>
1740 compute_transfinite_interpolation(
const AccessorType &cell,
1742 const bool cell_is_flat)
1744 const unsigned int dim = AccessorType::dimension;
1745 const unsigned int spacedim = AccessorType::space_dimension;
1753 const std::array<Point<spacedim>, 4>
vertices{
1754 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1759 std::array<double, 4> weights_vertices{
1760 {(1. - chart_point[0]) * (1. - chart_point[1]),
1761 chart_point[0] * (1. - chart_point[1]),
1762 (1. - chart_point[0]) * chart_point[1],
1763 chart_point[0] * chart_point[1]}};
1768 new_point += weights_vertices[v] *
vertices[v];
1780 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1783 const auto weights_view =
1785 const auto points_view =
make_array_view(points.begin(), points.end());
1787 for (
unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1790 const double my_weight =
1791 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1792 const double line_point = chart_point[1 - line / 2];
1799 cell.line(line)->manifold_id();
1800 if (line_manifold_id == my_manifold_id ||
1805 my_weight * (1. - line_point);
1808 my_weight * line_point;
1816 weights[0] = 1. - line_point;
1817 weights[1] = line_point;
1820 .get_new_point(points_view, weights_view);
1826 new_point -= weights_vertices[v] *
vertices[v];
1835 static constexpr
unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1845 static constexpr
unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1853 template <
typename AccessorType>
1855 compute_transfinite_interpolation(
const AccessorType &cell,
1857 const bool cell_is_flat)
1859 const unsigned int dim = AccessorType::dimension;
1860 const unsigned int spacedim = AccessorType::space_dimension;
1866 const std::array<Point<spacedim>, 8>
vertices{{cell.vertex(0),
1878 double linear_shapes[10];
1879 for (
unsigned int d = 0;
d < 3; ++
d)
1881 linear_shapes[2 *
d] = 1. - chart_point[
d];
1882 linear_shapes[2 *
d + 1] = chart_point[
d];
1886 for (
unsigned int d = 6;
d < 10; ++
d)
1887 linear_shapes[
d] = linear_shapes[
d - 6];
1889 std::array<double, 8> weights_vertices;
1890 for (
unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1891 for (
unsigned int i1 = 0; i1 < 2; ++i1)
1892 for (
unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1893 weights_vertices[v] =
1894 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1898 for (
unsigned int v = 0; v < 8; ++v)
1899 new_point += weights_vertices[v] *
vertices[v];
1905 std::array<double, GeometryInfo<3>::lines_per_cell> weights_lines;
1906 std::fill(weights_lines.begin(), weights_lines.end(), 0.0);
1909 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1912 const auto weights_view =
1914 const auto points_view =
make_array_view(points.begin(), points.end());
1918 const double my_weight = linear_shapes[face];
1919 const unsigned int face_even = face - face % 2;
1921 if (std::abs(my_weight) < 1
e-13)
1927 cell.face(face)->manifold_id();
1928 if (face_manifold_id == my_manifold_id ||
1931 for (
unsigned int line = 0;
1932 line < GeometryInfo<2>::lines_per_cell;
1935 const double line_weight =
1936 linear_shapes[face_even + 2 + line];
1937 weights_lines[face_to_cell_lines_3d[face][line]] +=
1938 my_weight * line_weight;
1944 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1945 linear_shapes[face_even + 2] *
1946 (linear_shapes[face_even + 4] * my_weight);
1947 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1948 linear_shapes[face_even + 3] *
1949 (linear_shapes[face_even + 4] * my_weight);
1950 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1951 linear_shapes[face_even + 2] *
1952 (linear_shapes[face_even + 5] * my_weight);
1953 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1954 linear_shapes[face_even + 3] *
1955 (linear_shapes[face_even + 5] * my_weight);
1960 points[v] =
vertices[face_to_cell_vertices_3d[face][v]];
1962 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1964 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1966 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1968 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1971 .get_new_point(points_view, weights_view);
1976 const auto weights_view_line =
1978 const auto points_view_line =
1980 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1983 const double line_point =
1984 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1985 double my_weight = 0.;
1987 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1990 const unsigned int subline = line - 8;
1992 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1994 my_weight -= weights_lines[line];
1996 if (std::abs(my_weight) < 1
e-13)
2000 cell.line(line)->manifold_id();
2001 if (line_manifold_id == my_manifold_id ||
2006 my_weight * (1. - line_point);
2009 my_weight * (line_point);
2017 weights[0] = 1. - line_point;
2018 weights[1] = line_point;
2020 .get_new_point(points_view_line,
2027 new_point += weights_vertices[v] *
vertices[v];
2035 template <
int dim,
int spacedim>
2047 ExcMessage(
"chart_point is not in unit interval"));
2049 return compute_transfinite_interpolation(*cell,
2051 coarse_cell_is_flat[cell->index()]);
2056 template <
int dim,
int spacedim>
2065 for (
unsigned int d = 0;
d < dim; ++
d)
2068 const double step = chart_point[
d] > 0.5 ? -1
e-8 : 1
e-8;
2071 modified[
d] += step;
2073 compute_transfinite_interpolation(*cell,
2075 coarse_cell_is_flat[cell->index()]) -
2076 pushed_forward_chart_point;
2077 for (
unsigned int e = 0;
e < spacedim; ++
e)
2078 grad[
e][
d] = difference[
e] / step;
2085 template <
int dim,
int spacedim>
2093 for (
unsigned int d = 0;
d < dim; ++
d)
2107 compute_transfinite_interpolation(*cell,
2109 coarse_cell_is_flat[cell->index()]);
2110 const double tolerance = 1
e-21 * Utilities::fixed_power<2>(cell->diameter());
2111 double residual_norm_square = residual.
norm_square();
2113 bool must_recompute_jacobian =
true;
2114 for (
unsigned int i = 0; i < 100; ++i)
2116 if (residual_norm_square < tolerance)
2123 for (
unsigned int d = 0;
d < spacedim; ++
d)
2124 for (
unsigned int e = 0;
e < dim; ++
e)
2125 update[
e] += inv_grad[
d][
e] * residual[
d];
2126 return chart_point + update;
2138 if (must_recompute_jacobian || i % 9 == 0)
2146 push_forward_gradient(cell,
2152 must_recompute_jacobian =
false;
2155 for (
unsigned int d = 0;
d < spacedim; ++
d)
2156 for (
unsigned int e = 0;
e < dim; ++
e)
2157 update[
e] += inv_grad[
d][
e] * residual[
d];
2171 while (alpha > 1
e-4)
2173 Point<dim> guess = chart_point + alpha * update;
2175 point - compute_transfinite_interpolation(
2176 *cell, guess, coarse_cell_is_flat[cell->index()]);
2177 const double residual_norm_new = residual_guess.
norm_square();
2178 if (residual_norm_new < residual_norm_square)
2180 residual = residual_guess;
2181 residual_norm_square = residual_norm_new;
2182 chart_point += alpha * update;
2195 must_recompute_jacobian =
true;
2209 for (
unsigned int d = 0;
d < spacedim; ++
d)
2210 for (
unsigned int e = 0;
e < dim; ++
e)
2211 Jinv_deltaf[
e] += inv_grad[
d][
e] * delta_f[
d];
2218 if (std::abs(delta_x * Jinv_deltaf) > 1
e-12 && !must_recompute_jacobian)
2221 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
2223 for (
unsigned int d = 0;
d < spacedim; ++
d)
2224 for (
unsigned int e = 0;
e < dim; ++
e)
2225 jac_update[
d] += delta_x[
e] * inv_grad[
d][
e];
2226 for (
unsigned int d = 0;
d < spacedim; ++
d)
2227 for (
unsigned int e = 0;
e < dim; ++
e)
2228 inv_grad[
d][
e] += factor[
e] * jac_update[
d];
2236 template <
int dim,
int spacedim>
2237 std::array<unsigned int, 20>
2247 ExcMessage(
"The manifold was initialized with level " +
2249 "active cells on a lower level. Coarsening the mesh is " +
2250 "currently not supported"));
2260 boost::container::small_vector<std::pair<double, unsigned int>, 200>
2261 distances_and_cells;
2262 for (; cell != endc; ++cell)
2272 vertices[vertex_n] = cell->vertex(vertex_n);
2282 double radius_square = 0.;
2286 bool inside_circle =
true;
2287 for (
unsigned int i = 0; i < points.size(); ++i)
2288 if ((
center - points[i]).norm_square() > radius_square * 1.5)
2290 inside_circle =
false;
2293 if (inside_circle ==
false)
2297 double current_distance = 0;
2298 for (
unsigned int i = 0; i < points.size(); ++i)
2301 quadratic_approximation[cell->index()].compute(points[i]);
2304 distances_and_cells.push_back(
2305 std::make_pair(current_distance, cell->index()));
2310 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2311 std::array<unsigned int, 20> cells;
2313 for (
unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2315 cells[i] = distances_and_cells[i].second;
2322 template <
int dim,
int spacedim>
2328 Assert(surrounding_points.size() == chart_points.size(),
2329 ExcMessage(
"The chart points array view must be as large as the "
2330 "surrounding points array view."));
2332 std::array<unsigned int, 20> nearby_cells =
2333 get_possible_cells_around_points(surrounding_points);
2357 auto guess_chart_point_structdim_2 = [&](
const unsigned int i) ->
Point<dim> {
2358 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2359 ExcMessage(
"This function assumes that there are eight surrounding "
2360 "points around a two-dimensional object. It also assumes "
2361 "that the first three chart points have already been "
2371 return chart_points[1] + (chart_points[2] - chart_points[0]);
2373 return 0.5 * (chart_points[0] + chart_points[2]);
2375 return 0.5 * (chart_points[1] + chart_points[3]);
2377 return 0.5 * (chart_points[0] + chart_points[1]);
2379 return 0.5 * (chart_points[2] + chart_points[3]);
2408 auto guess_chart_point_structdim_3 = [&](
const unsigned int i) ->
Point<dim> {
2409 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2410 ExcMessage(
"This function assumes that there are eight surrounding "
2411 "points around a three-dimensional object. It also "
2412 "assumes that the first five chart points have already "
2414 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2418 bool use_structdim_2_guesses =
false;
2419 bool use_structdim_3_guesses =
false;
2424 if (surrounding_points.size() == 8)
2427 surrounding_points[6] - surrounding_points[0];
2429 surrounding_points[7] - surrounding_points[2];
2437 use_structdim_2_guesses =
true;
2438 else if (spacedim == 3)
2441 use_structdim_3_guesses =
true;
2444 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2445 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2450 auto compute_chart_point =
2452 const unsigned int point_index) {
2458 bool used_quadratic_approximation =
false;
2461 if (point_index == 3 && surrounding_points.size() >= 8)
2462 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2463 else if (use_structdim_2_guesses && 3 < point_index)
2464 guess = guess_chart_point_structdim_2(point_index);
2465 else if (use_structdim_3_guesses && 4 < point_index)
2466 guess = guess_chart_point_structdim_3(point_index);
2467 else if (dim == 3 && point_index > 7 && surrounding_points.size() == 26)
2469 if (point_index < 20)
2472 point_index - 8, 0)] +
2474 point_index - 8, 1)]);
2478 point_index - 20, 0)] +
2480 point_index - 20, 1)] +
2482 point_index - 20, 2)] +
2484 point_index - 20, 3)]);
2488 guess = quadratic_approximation[cell->index()].compute(
2489 surrounding_points[point_index]);
2490 used_quadratic_approximation =
true;
2492 chart_points[point_index] =
2493 pull_back(cell, surrounding_points[point_index], guess);
2498 if (chart_points[point_index][0] ==
2500 !used_quadratic_approximation)
2502 guess = quadratic_approximation[cell->index()].compute(
2503 surrounding_points[point_index]);
2504 chart_points[point_index] =
2505 pull_back(cell, surrounding_points[point_index], guess);
2508 if (chart_points[point_index][0] ==
2511 for (
unsigned int d = 0;
d < dim; ++
d)
2513 chart_points[point_index] =
2514 pull_back(cell, surrounding_points[point_index], guess);
2519 for (
unsigned int c = 0; c < nearby_cells.size(); ++c)
2523 bool inside_unit_cell =
true;
2524 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2526 compute_chart_point(cell, i);
2533 inside_unit_cell =
false;
2537 if (inside_unit_cell ==
true)
2544 if (c == nearby_cells.size() - 1 ||
2549 std::ostringstream message;
2550 for (
unsigned int b = 0;
b <= c; ++
b)
2554 message <<
"Looking at cell " << cell->id()
2555 <<
" with vertices: " << std::endl;
2557 message << cell->vertex(v) <<
" ";
2558 message << std::endl;
2559 message <<
"Transformation to chart coordinates: " << std::endl;
2560 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2562 compute_chart_point(cell, i);
2563 message << surrounding_points[i] <<
" -> " << chart_points[i]
2583 template <
int dim,
int spacedim>
2589 boost::container::small_vector<Point<dim>, 100> chart_points(
2590 surrounding_points.size());
2593 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2596 chart_manifold.get_new_point(chart_points_view, weights);
2603 template <
int dim,
int spacedim>
2613 boost::container::small_vector<Point<dim>, 100> chart_points(
2614 surrounding_points.size());
2617 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2619 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2621 chart_manifold.get_new_points(chart_points_view,
2624 new_points_on_chart.end()));
2626 for (
unsigned int row = 0; row < weights.size(0); ++row)
2627 new_points[row] =
push_forward(cell, new_points_on_chart[row]);
2633 #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 std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) 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
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
const std::vector< Point< dim > > & get_points() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() 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 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
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &vertices, const ArrayView< const double > &weights) const override
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
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
numbers::NumberTraits< Number >::real_type norm() 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)
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
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
std::array< unsigned int, 20 > get_possible_cells_around_points(const ArrayView< const Point< spacedim >> &surrounding_points) const
Triangulation< dim, spacedim >::cell_iterator compute_chart_points(const ArrayView< const Point< spacedim >> &surrounding_points, ArrayView< Point< dim >> chart_points) const
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()
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcEmptyObject()
#define AssertIndexRange(index, range)
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
Expression atan2(const Expression &y, const Expression &x)
Expression acos(const Expression &x)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Point< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number signed_angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b, const Tensor< 1, spacedim, Number > &axis)
long double gamma(const unsigned int n)
Point< spacedim > compute_normal(const Tensor< 1, spacedim > &, bool=false)
Tensor< 1, 3 > apply_exponential_map(const Tensor< 1, 3 > &u, const Tensor< 1, 3 > &dir)
static constexpr double invalid_pull_back_coordinate
Tensor< 1, 3 > projected_direction(const Tensor< 1, 3 > &u, const Tensor< 1, 3 > &v)
static constexpr double PI
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > sqrt(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 Point< dim, Number > project_to_unit_cell(const Point< dim, Number > &p)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
const ::Triangulation< dim, spacedim > & tria