16 #include <deal.II/base/std_cxx14/memory.h> 17 #include <deal.II/base/table.h> 18 #include <deal.II/base/tensor.h> 20 #include <deal.II/fe/mapping.h> 22 #include <deal.II/grid/manifold_lib.h> 23 #include <deal.II/grid/tria.h> 24 #include <deal.II/grid/tria_accessor.h> 25 #include <deal.II/grid/tria_iterator.h> 27 #include <deal.II/lac/vector.h> 31 DEAL_II_NAMESPACE_OPEN
40 static constexpr
double invalid_pull_back_coordinate = 20.0;
48 const double theta = dir.
norm();
56 const Tensor<1, 3> tmp = cos(theta) * u + sin(theta) * dirUnit;
57 return tmp / tmp.
norm();
74 template <
int spacedim>
83 compute_normal(
const Tensor<1, 3> &vector,
bool normalize =
false)
86 ExcMessage(
"The direction parameter must not be zero!"));
88 if (std::abs(vector[0]) >= std::abs(vector[1]) &&
89 std::abs(vector[0]) >= std::abs(vector[2]))
93 normal[0] = (vector[1] + vector[2]) / vector[0];
95 else if (std::abs(vector[1]) >= std::abs(vector[0]) &&
96 std::abs(vector[1]) >= std::abs(vector[2]))
100 normal[1] = (vector[0] + vector[2]) / vector[1];
106 normal[2] = (vector[0] + vector[1]) / vector[2];
109 normal /= normal.
norm();
120 template <
int dim,
int spacedim>
129 template <
int dim,
int spacedim>
130 std::unique_ptr<Manifold<dim, spacedim>>
133 return std_cxx14::make_unique<PolarManifold<dim, spacedim>>(
center);
138 template <
int dim,
int spacedim>
155 template <
int dim,
int spacedim>
160 Assert(spherical_point[0] >= 0.0,
161 ExcMessage(
"Negative radius for given point."));
162 const double rho = spherical_point[0];
163 const double theta = spherical_point[1];
170 p[0] = rho * cos(theta);
171 p[1] = rho * sin(theta);
175 const double phi = spherical_point[2];
176 p[0] = rho * sin(theta) * cos(phi);
177 p[1] = rho * sin(theta) * sin(phi);
178 p[2] = rho * cos(theta);
189 template <
int dim,
int spacedim>
195 const double rho = R.
norm();
204 p[1] = atan2(R[1], R[0]);
212 const double z = R[2];
213 p[2] = atan2(R[1], R[0]);
216 p[1] = atan2(sqrt(R[0] * R[0] + R[1] * R[1]), z);
228 template <
int dim,
int spacedim>
233 Assert(spherical_point[0] >= 0.0,
234 ExcMessage(
"Negative radius for given point."));
235 const double rho = spherical_point[0];
236 const double theta = spherical_point[1];
244 DX[0][0] = cos(theta);
245 DX[0][1] = -rho * sin(theta);
246 DX[1][0] = sin(theta);
247 DX[1][1] = rho * cos(theta);
253 const double phi = spherical_point[2];
254 DX[0][0] = sin(theta) * cos(phi);
255 DX[0][1] = rho * cos(theta) * cos(phi);
256 DX[0][2] = -rho * sin(theta) * sin(phi);
258 DX[1][0] = sin(theta) * sin(phi);
259 DX[1][1] = rho * cos(theta) * sin(phi);
260 DX[1][2] = rho * sin(theta) * cos(phi);
262 DX[2][0] = cos(theta);
263 DX[2][1] = -rho * sin(theta);
280 template <
int dim,
int spacedim>
284 , polar_manifold(center)
289 template <
int dim,
int spacedim>
290 std::unique_ptr<Manifold<dim, spacedim>>
293 return std_cxx14::make_unique<SphericalManifold<dim, spacedim>>(
center);
298 template <
int dim,
int spacedim>
303 const double w)
const 305 const double tol = 1e-10;
307 if ((p1 - p2).norm_square() < tol * tol || std::abs(w) < tol)
309 else if (std::abs(w - 1.0) < tol)
319 const double r1 = v1.
norm();
320 const double r2 = v2.
norm();
322 Assert(r1 > tol && r2 > tol,
323 ExcMessage(
"p1 and p2 cannot coincide with the center."));
329 const double cosgamma = e1 * e2;
332 if (cosgamma < -1 + 8. * std::numeric_limits<double>::epsilon())
336 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
342 const double sigma = w * std::acos(cosgamma);
347 const double n_norm = n.
norm();
350 "Probably, this means v1==v2 or v2==0."));
364 template <
int dim,
int spacedim>
370 const double tol = 1e-10;
377 const double r1 = v1.
norm();
378 const double r2 = v2.
norm();
388 const double cosgamma = e1 * e2;
390 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
391 ExcMessage(
"p1 and p2 cannot lie on the same diameter and be opposite " 392 "respect to the center."));
394 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
400 const double n_norm = n.
norm();
403 "Probably, this means v1==v2 or v2==0."));
409 const double gamma = std::acos(cosgamma);
410 return (r2 - r1) * e1 + r1 * gamma * n;
415 template <
int dim,
int spacedim>
426 std::array<double, n_vertices> distances_to_center;
427 std::array<double, n_vertices - 1> distances_to_first_vertex;
428 distances_to_center[0] = (face->vertex(0) -
center).norm_square();
429 for (
unsigned int i = 1; i < n_vertices; ++i)
431 distances_to_center[i] = (face->vertex(i) -
center).norm_square();
432 distances_to_first_vertex[i - 1] =
433 (face->vertex(i) - face->vertex(0)).norm_square();
435 const auto minmax_distance =
436 std::minmax_element(distances_to_center.begin(), distances_to_center.end());
437 const auto min_distance_to_first_vertex =
438 std::min_element(distances_to_first_vertex.begin(),
439 distances_to_first_vertex.end());
441 if (*minmax_distance.second - *minmax_distance.first <
442 1.e-10 * *min_distance_to_first_vertex)
446 unnormalized_spherical_normal / unnormalized_spherical_normal.
norm();
447 return normalized_spherical_normal;
476 template <
int dim,
int spacedim>
488 std::array<double, n_vertices> distances_to_center;
489 std::array<double, n_vertices - 1> distances_to_first_vertex;
490 distances_to_center[0] = (face->vertex(0) -
center).norm_square();
491 for (
unsigned int i = 1; i < n_vertices; ++i)
493 distances_to_center[i] = (face->vertex(i) -
center).norm_square();
494 distances_to_first_vertex[i - 1] =
495 (face->vertex(i) - face->vertex(0)).norm_square();
497 const auto minmax_distance =
498 std::minmax_element(distances_to_center.begin(), distances_to_center.end());
499 const auto min_distance_to_first_vertex =
500 std::min_element(distances_to_first_vertex.begin(),
501 distances_to_first_vertex.end());
503 if (*minmax_distance.second - *minmax_distance.first <
504 1.e-10 * *min_distance_to_first_vertex)
506 for (
unsigned int vertex = 0; vertex < n_vertices; ++vertex)
507 face_vertex_normals[vertex] = face->vertex(vertex) -
center;
515 template <
int dim,
int spacedim>
525 get_new_points(surrounding_points, make_array_view(weights), new_points);
532 template <
int dim,
int spacedim>
543 make_array_view(&new_point, &new_point + 1));
550 template <
int dim,
int spacedim>
558 new_points.size() * surrounding_points.size());
559 const unsigned int weight_rows = new_points.size();
560 const unsigned int weight_columns = surrounding_points.size();
562 if (surrounding_points.size() == 2)
564 for (
unsigned int row = 0; row < weight_rows; ++row)
567 surrounding_points[1],
568 weights[row * weight_columns + 1]);
572 boost::container::small_vector<std::pair<double, Tensor<1, spacedim>>, 100>
573 new_candidates(new_points.size());
574 boost::container::small_vector<Tensor<1, spacedim>, 100> directions(
576 boost::container::small_vector<double, 100> distances(
577 surrounding_points.size(), 0.0);
578 double max_distance = 0.;
579 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
581 directions[i] = surrounding_points[i] -
center;
582 distances[i] = directions[i].
norm();
584 if (distances[i] != 0.)
585 directions[i] /= distances[i];
588 ExcMessage(
"One of the vertices coincides with the center. " 589 "This is not allowed!"));
593 for (
unsigned int k = 0; k < i; ++k)
595 const double squared_distance =
596 (directions[i] - directions[k]).norm_square();
597 max_distance = std::max(max_distance, squared_distance);
603 const double tolerance = 1e-10;
604 boost::container::small_vector<bool, 100> accurate_point_was_found(
605 new_points.size(),
false);
607 make_array_view(directions.begin(), directions.end());
609 make_array_view(distances.begin(), distances.end());
610 for (
unsigned int row = 0; row < weight_rows; ++row)
612 new_candidates[row] =
620 if (new_candidates[row].first == 0.0)
623 accurate_point_was_found[row] =
true;
645 if (max_distance < 2e-2)
647 for (
unsigned int row = 0; row < weight_rows; ++row)
649 center + new_candidates[row].first * new_candidates[row].second;
659 boost::container::small_vector<double, 1000> merged_weights(weights.
size());
660 boost::container::small_vector<Tensor<1, spacedim>, 100> merged_directions(
662 boost::container::small_vector<double, 100> merged_distances(
663 surrounding_points.size(), 0.0);
665 unsigned int n_unique_directions = 0;
666 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
668 bool found_duplicate =
false;
672 for (
unsigned int j = 0; j < n_unique_directions; ++j)
674 const double squared_distance =
675 (directions[i] - directions[j]).norm_square();
676 if (!found_duplicate && squared_distance < 1e-28)
678 found_duplicate =
true;
679 for (
unsigned int row = 0; row < weight_rows; ++row)
680 merged_weights[row * weight_columns + j] +=
681 weights[row * weight_columns + i];
685 if (found_duplicate ==
false)
687 merged_directions[n_unique_directions] = directions[i];
688 merged_distances[n_unique_directions] = distances[i];
689 for (
unsigned int row = 0; row < weight_rows; ++row)
690 merged_weights[row * weight_columns + n_unique_directions] =
691 weights[row * weight_columns + i];
693 ++n_unique_directions;
699 boost::container::small_vector<unsigned int, 100> merged_weights_index(
701 for (
unsigned int row = 0; row < weight_rows; ++row)
703 for (
unsigned int existing_row = 0; existing_row < row; ++existing_row)
705 bool identical_weights =
true;
707 for (
unsigned int weight_index = 0;
708 weight_index < n_unique_directions;
710 if (std::abs(merged_weights[row * weight_columns + weight_index] -
711 merged_weights[existing_row * weight_columns +
712 weight_index]) > tolerance)
714 identical_weights =
false;
718 if (identical_weights)
720 merged_weights_index[row] = existing_row;
729 make_array_view(merged_directions.begin(),
730 merged_directions.begin() + n_unique_directions);
732 make_array_view(merged_distances.begin(),
733 merged_distances.begin() + n_unique_directions);
735 for (
unsigned int row = 0; row < weight_rows; ++row)
736 if (!accurate_point_was_found[row])
741 &merged_weights[row * weight_columns], n_unique_directions);
742 new_candidates[row].second =
744 array_merged_distances,
745 array_merged_weights,
749 new_candidates[row].second =
750 new_candidates[merged_weights_index[row]].second;
753 center + new_candidates[row].first * new_candidates[row].second;
759 template <
int dim,
int spacedim>
760 std::pair<double, Tensor<1, spacedim>>
766 const double tolerance = 1e-10;
771 double total_weights = 0.;
772 for (
unsigned int i = 0; i < directions.size(); i++)
775 if (std::abs(1 - weights[i]) < tolerance)
776 return std::make_pair(distances[i], directions[i]);
778 rho += distances[i] * weights[i];
779 candidate += directions[i] * weights[i];
780 total_weights += weights[i];
784 const double norm = candidate.
norm();
788 rho /= total_weights;
790 return std::make_pair(rho, candidate);
796 template <
int spacedim>
819 Point<3> candidate = candidate_point;
820 const unsigned int n_merged_points = directions.size();
821 const double tolerance = 1e-10;
822 const int max_iterations = 10;
827 for (
unsigned int i = 0; i < n_merged_points; ++i)
829 const double squared_distance =
830 (candidate - directions[i]).norm_square();
831 if (squared_distance < tolerance * tolerance)
837 if (n_merged_points == 2)
840 Assert(std::abs(weights[0] + weights[1] - 1.0) < 1e-13,
856 for (
unsigned int i = 0; i < max_iterations; ++i)
861 const Tensor<1, 3> Clocalx = internal::compute_normal(candidate);
862 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
870 for (
unsigned int i = 0; i < n_merged_points; ++i)
871 if (std::abs(weights[i]) > 1.e-15)
873 vPerp = internal::projected_direction(directions[i], candidate);
875 const double sintheta = std::sqrt(sinthetaSq);
876 if (sintheta < tolerance)
878 Hessian[0][0] += weights[i];
879 Hessian[1][1] += weights[i];
883 const double costheta = (directions[i]) * candidate;
884 const double theta = atan2(sintheta, costheta);
885 const double sincthetaInv = theta / sintheta;
887 const double cosphi = vPerp * Clocalx;
888 const double sinphi = vPerp * Clocaly;
890 gradlocal[0] = cosphi;
891 gradlocal[1] = sinphi;
892 gradient += (weights[i] * sincthetaInv) * gradlocal;
894 const double wt = weights[i] / sinthetaSq;
895 const double sinphiSq = sinphi * sinphi;
896 const double cosphiSq = cosphi * cosphi;
897 const double tt = sincthetaInv * costheta;
898 const double offdiag = cosphi * sinphi * wt * (1.0 - tt);
899 Hessian[0][0] += wt * (cosphiSq + tt * sinphiSq);
900 Hessian[0][1] += offdiag;
901 Hessian[1][0] += offdiag;
902 Hessian[1][1] += wt * (sinphiSq + tt * cosphiSq);
910 const Tensor<1, 2> xDisplocal = inverse_Hessian * gradient;
912 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
915 const Point<3> candidateOld = candidate;
917 Point<3>(internal::apply_exponential_map(candidate, xDisp));
920 if ((candidate - candidateOld).norm_square() < tolerance * tolerance)
930 template <
int dim,
int spacedim>
950 const Point<3> & candidate_point)
const 952 return do_get_new_point(directions, distances, weights, candidate_point);
963 const Point<3> & candidate_point)
const 965 return do_get_new_point(directions, distances, weights, candidate_point);
976 const Point<3> & candidate_point)
const 978 return do_get_new_point(directions, distances, weights, candidate_point);
986 template <
int dim,
int spacedim>
988 const double tolerance)
996 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1001 template <
int dim,
int spacedim>
1015 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1020 template <
int dim,
int spacedim>
1021 std::unique_ptr<Manifold<dim, spacedim>>
1024 return std_cxx14::make_unique<CylindricalManifold<dim, spacedim>>(
1030 template <
int dim,
int spacedim>
1037 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1041 double average_length = 0.;
1042 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
1044 middle += surrounding_points[i] * weights[i];
1045 average_length += surrounding_points[i].square() * weights[i];
1048 const double lambda = middle *
direction;
1050 if ((middle - direction * lambda).square() <
tolerance * average_length)
1059 template <
int dim,
int spacedim>
1065 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1069 const double lambda = normalized_point *
direction;
1070 const Point<spacedim> projection = point_on_axis + direction * lambda;
1077 const double phi = std::atan2(det, dot);
1080 return Point<3>(p_diff.norm(), phi, lambda);
1085 template <
int dim,
int spacedim>
1091 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1098 const double sine_r = std::sin(chart_point(1)) * chart_point(0);
1099 const double cosine_r = std::cos(chart_point(1)) * chart_point(0);
1110 template <
int dim,
int spacedim>
1116 ExcMessage(
"CylindricalManifold can only be used for spacedim==3!"));
1125 const double sine = std::sin(chart_point(1));
1126 const double cosine = std::cos(chart_point(1));
1132 derivatives[0][0] = intermediate[0];
1133 derivatives[1][0] = intermediate[1];
1134 derivatives[2][0] = intermediate[2];
1154 template <
int dim,
int spacedim,
int chartdim>
1161 , push_forward_function(&push_forward_function)
1162 , pull_back_function(&pull_back_function)
1163 , tolerance(tolerance)
1164 , owns_pointers(false)
1165 , finite_difference_step(0)
1173 template <
int dim,
int spacedim,
int chartdim>
1184 , const_map(const_map)
1185 , tolerance(tolerance)
1187 , push_forward_expression(push_forward_expression)
1188 , pull_back_expression(pull_back_expression)
1189 , chart_vars(chart_vars)
1190 , space_vars(space_vars)
1195 pf->
initialize(chart_vars, push_forward_expression, const_map);
1196 pb->
initialize(space_vars, pull_back_expression, const_map);
1203 template <
int dim,
int spacedim,
int chartdim>
1220 template <
int dim,
int spacedim,
int chartdim>
1221 std::unique_ptr<Manifold<dim, spacedim>>
1236 return std_cxx14::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1247 return std_cxx14::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1256 template <
int dim,
int spacedim,
int chartdim>
1264 for (
unsigned int i = 0; i < spacedim; ++i)
1270 for (
unsigned int i = 0; i < chartdim; ++i)
1273 (std::abs(pb[i] - chart_point[i]) <
tolerance * chart_point.
norm())) ||
1274 (std::abs(pb[i] - chart_point[i]) <
tolerance),
1276 "The push forward is not the inverse of the pull back! Bailing out."));
1284 template <
int dim,
int spacedim,
int chartdim>
1290 for (
unsigned int i = 0; i < spacedim; ++i)
1293 for (
unsigned int j = 0; j < chartdim; ++j)
1294 DF[i][j] = gradient[j];
1301 template <
int dim,
int spacedim,
int chartdim>
1309 for (
unsigned int i = 0; i < chartdim; ++i)
1326 double phi = atan2(y, x);
1327 double theta = atan2(z, std::sqrt(x * x + y * y) - R);
1329 std::sqrt(pow(y - sin(phi) * R, 2.0) + pow(x - cos(phi) * R, 2.0) + z * z) /
1340 double phi = chart_point(0);
1341 double theta = chart_point(1);
1342 double w = chart_point(2);
1344 return Point<3>(cos(phi) * R + r * w * cos(theta) * cos(phi),
1346 sin(phi) * R + r * w * cos(theta) * sin(phi));
1358 ExcMessage(
"Outer radius R must be greater than the inner " 1366 std::unique_ptr<Manifold<dim, 3>>
1369 return std_cxx14::make_unique<TorusManifold<dim>>(R, r);
1380 double phi = chart_point(0);
1381 double theta = chart_point(1);
1382 double w = chart_point(2);
1384 DX[0][0] = -sin(phi) * R - r * w * cos(theta) * sin(phi);
1385 DX[0][1] = -r * w * sin(theta) * cos(phi);
1386 DX[0][2] = r * cos(theta) * cos(phi);
1389 DX[1][1] = r * w * cos(theta);
1390 DX[1][2] = r * sin(theta);
1392 DX[2][0] = cos(phi) * R + r * w * cos(theta) * cos(phi);
1393 DX[2][1] = -r * w * sin(theta) * sin(phi);
1394 DX[2][2] = r * cos(theta) * sin(phi);
1404 template <
int dim,
int spacedim>
1407 : triangulation(nullptr)
1415 template <
int dim,
int spacedim>
1425 template <
int dim,
int spacedim>
1426 std::unique_ptr<Manifold<dim, spacedim>>
1432 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1437 template <
int dim,
int spacedim>
1445 this->triangulation =
nullptr;
1453 for (; cell != endc; ++cell)
1455 bool cell_is_flat =
true;
1456 for (
unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1457 if (cell->line(l)->manifold_id() != cell->manifold_id() &&
1459 cell_is_flat =
false;
1461 for (
unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1462 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1464 cell_is_flat =
false;
1476 template <
typename AccessorType>
1478 compute_transfinite_interpolation(
const AccessorType &cell,
1482 return cell.vertex(0) * (1. - chart_point[0]) +
1483 cell.vertex(1) * chart_point[0];
1487 template <
typename AccessorType>
1489 compute_transfinite_interpolation(
const AccessorType &cell,
1491 const bool cell_is_flat)
1493 const unsigned int dim = AccessorType::dimension;
1494 const unsigned int spacedim = AccessorType::space_dimension;
1502 const std::array<Point<spacedim>, 4> vertices{
1503 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1508 std::array<double, 4> weights_vertices{
1509 {(1. - chart_point[0]) * (1. - chart_point[1]),
1510 chart_point[0] * (1. - chart_point[1]),
1511 (1. - chart_point[0]) * chart_point[1],
1512 chart_point[0] * chart_point[1]}};
1516 for (
unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
1517 new_point += weights_vertices[v] * vertices[v];
1529 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1532 const auto weights_view =
1533 make_array_view(weights.begin(), weights.end());
1534 const auto points_view = make_array_view(points.begin(), points.end());
1536 for (
unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1539 const double my_weight =
1540 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1541 const double line_point = chart_point[1 - line / 2];
1548 cell.line(line)->manifold_id();
1549 if (line_manifold_id == my_manifold_id ||
1554 my_weight * (1. - line_point);
1557 my_weight * line_point;
1565 weights[0] = 1. - line_point;
1566 weights[1] = line_point;
1569 .get_new_point(points_view, weights_view);
1574 for (
unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
1575 new_point -= weights_vertices[v] * vertices[v];
1584 static constexpr
unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1594 static constexpr
unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1602 template <
typename AccessorType>
1604 compute_transfinite_interpolation(
const AccessorType &cell,
1606 const bool cell_is_flat)
1608 const unsigned int dim = AccessorType::dimension;
1609 const unsigned int spacedim = AccessorType::space_dimension;
1615 const std::array<Point<spacedim>, 8> vertices{{cell.vertex(0),
1627 double linear_shapes[10];
1628 for (
unsigned int d = 0; d < 3; ++d)
1630 linear_shapes[2 * d] = 1. - chart_point[d];
1631 linear_shapes[2 * d + 1] = chart_point[d];
1635 for (
unsigned int d = 6; d < 10; ++d)
1636 linear_shapes[d] = linear_shapes[d - 6];
1638 std::array<double, 8> weights_vertices;
1639 for (
unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1640 for (
unsigned int i1 = 0; i1 < 2; ++i1)
1641 for (
unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1642 weights_vertices[v] =
1643 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1647 for (
unsigned int v = 0; v < 8; ++v)
1648 new_point += weights_vertices[v] * vertices[v];
1655 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1657 weights_lines[line] = 0;
1660 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1663 const auto weights_view =
1664 make_array_view(weights.begin(), weights.end());
1665 const auto points_view = make_array_view(points.begin(), points.end());
1667 for (
unsigned int face = 0; face < GeometryInfo<3>::faces_per_cell;
1670 const double my_weight = linear_shapes[face];
1671 const unsigned int face_even = face - face % 2;
1673 if (std::abs(my_weight) < 1e-13)
1679 cell.face(face)->manifold_id();
1680 if (face_manifold_id == my_manifold_id ||
1683 for (
unsigned int line = 0;
1684 line < GeometryInfo<2>::lines_per_cell;
1687 const double line_weight =
1688 linear_shapes[face_even + 2 + line];
1689 weights_lines[face_to_cell_lines_3d[face][line]] +=
1690 my_weight * line_weight;
1696 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1697 linear_shapes[face_even + 2] *
1698 (linear_shapes[face_even + 4] * my_weight);
1699 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1700 linear_shapes[face_even + 3] *
1701 (linear_shapes[face_even + 4] * my_weight);
1702 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1703 linear_shapes[face_even + 2] *
1704 (linear_shapes[face_even + 5] * my_weight);
1705 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1706 linear_shapes[face_even + 3] *
1707 (linear_shapes[face_even + 5] * my_weight);
1711 for (
unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell;
1713 points[v] = vertices[face_to_cell_vertices_3d[face][v]];
1715 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1717 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1719 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1721 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1724 .get_new_point(points_view, weights_view);
1729 const auto weights_view_line =
1730 make_array_view(weights.begin(), weights.begin() + 2);
1731 const auto points_view_line =
1732 make_array_view(points.begin(), points.begin() + 2);
1733 for (
unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1736 const double line_point =
1737 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1738 double my_weight = 0.;
1740 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1743 const unsigned int subline = line - 8;
1745 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1747 my_weight -= weights_lines[line];
1749 if (std::abs(my_weight) < 1e-13)
1753 cell.line(line)->manifold_id();
1754 if (line_manifold_id == my_manifold_id ||
1759 my_weight * (1. - line_point);
1762 my_weight * (line_point);
1770 weights[0] = 1. - line_point;
1771 weights[1] = line_point;
1772 new_point -= my_weight * tria.
get_manifold(line_manifold_id)
1773 .get_new_point(points_view_line,
1779 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
1780 new_point += weights_vertices[v] * vertices[v];
1788 template <
int dim,
int spacedim>
1800 ExcMessage(
"chart_point is not in unit interval"));
1802 return compute_transfinite_interpolation(*cell,
1809 template <
int dim,
int spacedim>
1818 for (
unsigned int d = 0; d < dim; ++d)
1821 const double step = chart_point[d] > 0.5 ? -1e-8 : 1e-8;
1824 modified[d] += step;
1826 compute_transfinite_interpolation(*cell,
1829 pushed_forward_chart_point;
1830 for (
unsigned int e = 0; e < spacedim; ++e)
1831 grad[e][d] = difference[e] / step;
1838 template <
int dim,
int spacedim>
1846 for (
unsigned int d = 0; d < dim; ++d)
1847 outside[d] = internal::invalid_pull_back_coordinate;
1860 compute_transfinite_interpolation(*cell,
1863 const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter());
1864 double residual_norm_square = residual.
norm_square();
1866 for (
unsigned int i = 0; i < 100; ++i)
1868 if (residual_norm_square < tolerance)
1875 for (
unsigned int d = 0; d < spacedim; ++d)
1876 for (
unsigned int e = 0; e < dim; ++e)
1877 update[e] += inv_grad[d][e] * residual[d];
1878 return chart_point + update;
1906 for (
unsigned int d = 0; d < spacedim; ++d)
1907 for (
unsigned int e = 0; e < dim; ++e)
1908 update[e] += inv_grad[d][e] * residual[d];
1922 while (alpha > 1e-7)
1924 Point<dim> guess = chart_point + alpha * update;
1926 point - compute_transfinite_interpolation(
1928 const double residual_norm_new = residual.
norm_square();
1929 if (residual_norm_new < residual_norm_square)
1931 residual_norm_square = residual_norm_new;
1932 chart_point += alpha * update;
1952 for (
unsigned int d = 0; d < spacedim; ++d)
1953 for (
unsigned int e = 0; e < dim; ++e)
1954 Jinv_deltaf[e] += inv_grad[d][e] * delta_f[d];
1956 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
1958 for (
unsigned int d = 0; d < spacedim; ++d)
1959 for (
unsigned int e = 0; e < dim; ++e)
1960 jac_update[d] += delta_x[e] * inv_grad[d][e];
1961 for (
unsigned int d = 0; d < spacedim; ++d)
1962 for (
unsigned int e = 0; e < dim; ++e)
1963 inv_grad[d][e] += factor[e] * jac_update[d];
1970 template <
int dim,
int spacedim>
1971 std::array<unsigned int, 20>
1981 ExcMessage(
"The manifold was initialized with level " +
1983 "active cells on a lower level. Coarsening the mesh is " +
1984 "currently not supported"));
1994 boost::container::small_vector<std::pair<double, unsigned int>, 200>
1995 distances_and_cells;
1996 for (; cell != endc; ++cell)
1999 if (&cell->get_manifold() !=
this)
2004 for (
unsigned int vertex_n = 0;
2005 vertex_n < GeometryInfo<dim>::vertices_per_cell;
2008 vertices[vertex_n] = cell->vertex(vertex_n);
2015 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2016 center += vertices[v];
2018 double radius_square = 0.;
2019 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
2021 std::max(radius_square, (center - vertices[v]).norm_square());
2022 bool inside_circle =
true;
2023 for (
unsigned int i = 0; i < points.size(); ++i)
2024 if ((center - points[i]).norm_square() > radius_square * 1.5)
2026 inside_circle =
false;
2029 if (inside_circle ==
false)
2033 double current_distance = 0;
2034 for (
unsigned int i = 0; i < points.size(); ++i)
2037 cell->real_to_unit_cell_affine_approximation(points[i]);
2040 distances_and_cells.emplace_back(current_distance, cell->index());
2045 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2046 std::array<unsigned int, 20> cells;
2048 for (
unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2050 cells[i] = distances_and_cells[i].second;
2057 template <
int dim,
int spacedim>
2063 Assert(surrounding_points.size() == chart_points.size(),
2064 ExcMessage(
"The chart points array view must be as large as the " 2065 "surrounding points array view."));
2067 std::array<unsigned int, 20> nearby_cells =
2092 auto guess_chart_point_structdim_2 = [&](
const unsigned int i) ->
Point<dim> {
2093 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2094 ExcMessage(
"This function assumes that there are eight surrounding " 2095 "points around a two-dimensional object. It also assumes " 2096 "that the first three chart points have already been " 2106 return chart_points[1] + (chart_points[2] - chart_points[0]);
2108 return 0.5 * (chart_points[0] + chart_points[2]);
2110 return 0.5 * (chart_points[1] + chart_points[3]);
2112 return 0.5 * (chart_points[0] + chart_points[1]);
2114 return 0.5 * (chart_points[2] + chart_points[3]);
2143 auto guess_chart_point_structdim_3 = [&](
const unsigned int i) ->
Point<dim> {
2144 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2145 ExcMessage(
"This function assumes that there are eight surrounding " 2146 "points around a three-dimensional object. It also " 2147 "assumes that the first five chart points have already " 2149 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2153 bool use_structdim_2_guesses =
false;
2154 bool use_structdim_3_guesses =
false;
2159 if (surrounding_points.size() == 8)
2162 surrounding_points[6] - surrounding_points[0];
2164 surrounding_points[7] - surrounding_points[2];
2172 use_structdim_2_guesses =
true;
2173 else if (spacedim == 3)
2176 use_structdim_3_guesses =
true;
2179 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2180 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2184 for (
unsigned int c = 0; c < nearby_cells.size(); ++c)
2188 bool inside_unit_cell =
true;
2189 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2196 bool used_affine_approximation =
false;
2199 if (i == 3 && surrounding_points.size() == 8)
2200 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2201 else if (use_structdim_2_guesses && 3 < i)
2202 guess = guess_chart_point_structdim_2(i);
2203 else if (use_structdim_3_guesses && 4 < i)
2204 guess = guess_chart_point_structdim_3(i);
2207 guess = cell->real_to_unit_cell_affine_approximation(
2208 surrounding_points[i]);
2209 used_affine_approximation =
true;
2211 chart_points[i] =
pull_back(cell, surrounding_points[i], guess);
2216 if (chart_points[i][0] == internal::invalid_pull_back_coordinate &&
2217 !used_affine_approximation)
2219 guess = cell->real_to_unit_cell_affine_approximation(
2220 surrounding_points[i]);
2221 chart_points[i] =
pull_back(cell, surrounding_points[i], guess);
2229 inside_unit_cell =
false;
2233 if (inside_unit_cell ==
true)
2240 if (c == nearby_cells.size() - 1 ||
2245 std::ostringstream message;
2246 for (
unsigned int b = 0; b <= c; ++b)
2250 message <<
"Looking at cell " << cell->id()
2251 <<
" with vertices: " << std::endl;
2252 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
2254 message << cell->vertex(v) <<
" ";
2255 message << std::endl;
2256 message <<
"Transformation to chart coordinates: " << std::endl;
2257 for (
unsigned int i = 0; i < surrounding_points.size(); ++i)
2260 << surrounding_points[i] <<
" -> " 2262 surrounding_points[i],
2263 cell->real_to_unit_cell_affine_approximation(
2264 surrounding_points[i]))
2284 template <
int dim,
int spacedim>
2290 boost::container::small_vector<Point<dim>, 100> chart_points(
2291 surrounding_points.size());
2293 make_array_view(chart_points.begin(), chart_points.end());
2304 template <
int dim,
int spacedim>
2314 boost::container::small_vector<Point<dim>, 100> chart_points(
2315 surrounding_points.size());
2317 make_array_view(chart_points.begin(), chart_points.end());
2320 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2324 make_array_view(new_points_on_chart.begin(),
2325 new_points_on_chart.end()));
2327 for (
unsigned int row = 0; row < weights.
size(0); ++row)
2328 new_points[row] =
push_forward(cell, new_points_on_chart[row]);
2334 #include "manifold_lib.inst" 2336 DEAL_II_NAMESPACE_CLOSE
Point< dim > pull_back(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_guess) const
FlatManifold< dim > chart_manifold
std::array< unsigned int, 20 > get_possible_cells_around_points(const ArrayView< const Point< spacedim >> &surrounding_points) const
static::ExceptionBase & ExcTransformationFailed()
static const unsigned int invalid_unsigned_int
virtual Point< spacedim > push_forward(const Point< 3 > &chart_point) const override
#define AssertDimension(dim1, dim2)
const unsigned int n_components
Number determinant(const SymmetricTensor< 2, dim, Number > &)
cell_iterator end() const
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)
static Tensor< 1, spacedim > get_periodicity()
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
const Triangulation< dim, spacedim > * triangulation
virtual Tensor< 1, dim, RangeNumberType > gradient(const Point< dim > &p, const unsigned int component=0) const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
const Point< spacedim > point_on_axis
std::vector< bool > coarse_cell_is_flat
virtual std::unique_ptr< Manifold< dim, 3 > > clone() const override
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
#define AssertIndexRange(index, range)
const double finite_difference_step
const std::string pull_back_expression
Point< spacedim > push_forward(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &chart_point) const
virtual DerivativeForm< 1, 3, 3 > push_forward_gradient(const Point< 3 > &chart_point) const override
static::ExceptionBase & ExcNotInitialized()
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &vertices, const ArrayView< const double > &weights) const override
cell_iterator last() const
virtual Point< 3 > pull_back(const Point< 3 > &p) const override
Triangulation< dim, spacedim >::cell_iterator compute_chart_points(const ArrayView< const Point< spacedim >> &surrounding_points, ArrayView< Point< dim >> chart_points) const
numbers::NumberTraits< Number >::real_type norm() const
#define AssertThrow(cond, exc)
TorusManifold(const double R, const double r)
cell_iterator begin(const unsigned int level=0) const
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
void initialize(const Triangulation< dim, spacedim > &triangulation)
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
const Point< spacedim > center
static double distance_to_unit_cell(const Point< dim > &p)
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
virtual Point< 3 > push_forward(const Point< 3 > &chart_point) const override
static::ExceptionBase & ExcMessage(std::string arg1)
static::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual DerivativeForm< 1, 3, spacedim > push_forward_gradient(const Point< 3 > &chart_point) const override
CylindricalManifold(const unsigned int axis=0, const double tolerance=1e-10)
unsigned int n_cells() const
SmartPointer< const Function< spacedim >, FunctionManifold< dim, spacedim, chartdim > > pull_back_function
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
PolarManifold(const Point< spacedim > center=Point< spacedim >())
size_type size(const unsigned int i) const
const PolarManifold< spacedim > polar_manifold
#define Assert(cond, exc)
virtual ~FunctionManifold() override
const std::string chart_vars
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const override
Abstract base class for mapping classes.
virtual Point< 3 > pull_back(const Point< spacedim > &space_point) const override
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual Point< spacedim > push_forward(const Point< spacedim > &chart_point) const override
SmartPointer< const Function< chartdim >, FunctionManifold< dim, spacedim, chartdim > > push_forward_function
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
const Point< spacedim > center
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 Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
static Point< dim > project_to_unit_cell(const Point< dim > &p)
const types::manifold_id invalid_manifold_id
virtual ~TransfiniteInterpolationManifold() override
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
void initialize(const std::string &vars, const std::vector< std::string > &expressions, const ConstMap &constants, const bool time_dependent=false)
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const override
const std::string space_vars
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
const Tensor< 1, chartdim > & get_periodicity() const
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
static::ExceptionBase & ExcEmptyObject()
const Tensor< 1, spacedim > direction
const std::string push_forward_expression
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
const Tensor< 1, spacedim > normal_direction
static::ExceptionBase & ExcNotImplemented()
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
const FunctionParser< spacedim >::ConstMap const_map
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
boost::signals2::connection clear_signal
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
numbers::NumberTraits< Number >::real_type norm_square() const
SphericalManifold(const Point< spacedim > center=Point< spacedim >())
ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
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, FaceVertexNormals &face_vertex_normals) const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
TriaIterator< CellAccessor< dim, spacedim >> cell_iterator
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
static::ExceptionBase & ExcInternalError()
Triangulation< dim, spacedim > & get_triangulation()