deal.II version GIT relicensing-1809-g201751ccb9 2024-09-08 15:20:00+00:00
Searching...
No Matches
manifold_lib.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
4// Copyright (C) 2014 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/table.h>
16#include <deal.II/base/tensor.h>
17
18#include <deal.II/fe/mapping.h>
20
22#include <deal.II/grid/tria.h>
25
26#include <deal.II/lac/vector.h>
27
29
30#include <boost/container/small_vector.hpp>
31
32#include <cmath>
33#include <limits>
34#include <memory>
35
37
38
39namespace internal
40{
41 // The pull_back function fails regularly in the compute_chart_points
42 // method, and, instead of throwing an exception, returns a point outside
43 // the unit cell. The individual coordinates of that point are given by the
44 // value below.
45 static constexpr double invalid_pull_back_coordinate = 20.0;
46
47 // Rotate a given unit vector u around the axis dir
48 // where the angle is given by the length of dir.
49 // This is the exponential map for a sphere.
52 {
53 const double theta = dir.norm();
54 if (theta < 1.e-10)
55 {
56 return u;
57 }
58 else
59 {
60 const Tensor<1, 3> dirUnit = dir / theta;
61 const Tensor<1, 3> tmp =
62 std::cos(theta) * u + std::sin(theta) * dirUnit;
63 return tmp / tmp.norm();
64 }
65 }
66
67 // Returns the direction to go from v to u
68 // projected to the plane perpendicular to the unit vector v.
69 // This one is more stable when u and v are nearly equal.
72 {
73 Tensor<1, 3> ans = u - v;
74 ans -= (ans * v) * v;
75 return ans; // ans = (u-v) - ((u-v)*v)*v
76 }
77
78 // helper function to compute a vector orthogonal to a given one.
79 // does nothing unless spacedim == 3.
80 template <int spacedim>
83 bool /*normalize*/ = false)
84 {
85 return {};
86 }
87
89 compute_normal(const Tensor<1, 3> &vector, bool normalize = false)
90 {
91 Assert(vector.norm_square() != 0.,
92 ExcMessage("The direction parameter must not be zero!"));
93 Point<3> normal;
94 if (std::abs(vector[0]) >= std::abs(vector[1]) &&
95 std::abs(vector[0]) >= std::abs(vector[2]))
96 {
97 normal[1] = -1.;
98 normal[2] = -1.;
99 normal[0] = (vector[1] + vector[2]) / vector[0];
100 }
101 else if (std::abs(vector[1]) >= std::abs(vector[0]) &&
102 std::abs(vector[1]) >= std::abs(vector[2]))
103 {
104 normal[0] = -1.;
105 normal[2] = -1.;
106 normal[1] = (vector[0] + vector[2]) / vector[1];
107 }
108 else
109 {
110 normal[0] = -1.;
111 normal[1] = -1.;
112 normal[2] = (vector[0] + vector[1]) / vector[2];
113 }
114 if (normalize)
115 normal /= normal.norm();
116 return normal;
117 }
118} // namespace internal
119
120
121
122// ============================================================
123// PolarManifold
124// ============================================================
125
126template <int dim, int spacedim>
128 : ChartManifold<dim, spacedim, spacedim>(
129 PolarManifold<dim, spacedim>::get_periodicity())
130 , center(center)
131 , p_center(center)
132{}
133
134
135
136template <int dim, int spacedim>
137std::unique_ptr<Manifold<dim, spacedim>>
139{
140 return std::make_unique<PolarManifold<dim, spacedim>>(*this);
141}
142
143
144
145template <int dim, int spacedim>
148{
149 Tensor<1, spacedim> periodicity;
150 // In two dimensions, theta is periodic.
151 // In three dimensions things are a little more complicated, since the only
152 // variable that is truly periodic is phi, while theta should be bounded
153 // between 0 and pi. There is currently no way to enforce this, so here we
154 // only fix periodicity for the last variable, corresponding to theta in 2d
155 // and phi in 3d.
156 periodicity[spacedim - 1] = 2 * numbers::PI;
157 return periodicity;
158}
159
160
161
162template <int dim, int spacedim>
165 const Point<spacedim> &spherical_point) const
166{
167 Assert(spherical_point[0] >= 0.0,
168 ExcMessage("Negative radius for given point."));
169 const double rho = spherical_point[0];
170 const double theta = spherical_point[1];
171
173 if (rho > 1e-10)
174 switch (spacedim)
175 {
176 case 2:
177 p[0] = rho * std::cos(theta);
178 p[1] = rho * std::sin(theta);
179 break;
180 case 3:
181 {
182 const double phi = spherical_point[2];
183 p[0] = rho * std::sin(theta) * std::cos(phi);
184 p[1] = rho * std::sin(theta) * std::sin(phi);
185 p[2] = rho * std::cos(theta);
186 break;
187 }
188 default:
190 }
191 return p + p_center;
192}
193
194
195
196template <int dim, int spacedim>
199 const Point<spacedim> &space_point) const
200{
201 const Tensor<1, spacedim> R = space_point - p_center;
202 const double rho = R.norm();
203
205 p[0] = rho;
206
207 switch (spacedim)
208 {
209 case 2:
210 {
211 p[1] = std::atan2(R[1], R[0]);
212 if (p[1] < 0)
213 p[1] += 2 * numbers::PI;
214 break;
215 }
216
217 case 3:
218 {
219 const double z = R[2];
220 p[2] = std::atan2(R[1], R[0]); // phi
221 if (p[2] < 0)
222 p[2] += 2 * numbers::PI; // phi is periodic
223 p[1] = std::atan2(std::sqrt(R[0] * R[0] + R[1] * R[1]), z); // theta
224 break;
225 }
226
227 default:
229 }
230 return p;
231}
232
233
234
235template <int dim, int spacedim>
238 const Point<spacedim> &spherical_point) const
239{
240 Assert(spherical_point[0] >= 0.0,
241 ExcMessage("Negative radius for given point."));
242 const double rho = spherical_point[0];
243 const double theta = spherical_point[1];
244
246 if (rho > 1e-10)
247 switch (spacedim)
248 {
249 case 2:
250 {
251 DX[0][0] = std::cos(theta);
252 DX[0][1] = -rho * std::sin(theta);
253 DX[1][0] = std::sin(theta);
254 DX[1][1] = rho * std::cos(theta);
255 break;
256 }
257
258 case 3:
259 {
260 const double phi = spherical_point[2];
261 DX[0][0] = std::sin(theta) * std::cos(phi);
262 DX[0][1] = rho * std::cos(theta) * std::cos(phi);
263 DX[0][2] = -rho * std::sin(theta) * std::sin(phi);
264
265 DX[1][0] = std::sin(theta) * std::sin(phi);
266 DX[1][1] = rho * std::cos(theta) * std::sin(phi);
267 DX[1][2] = rho * std::sin(theta) * std::cos(phi);
268
269 DX[2][0] = std::cos(theta);
270 DX[2][1] = -rho * std::sin(theta);
271 DX[2][2] = 0;
272 break;
273 }
274
275 default:
277 }
278 return DX;
279}
280
281
282
283namespace
284{
285 template <int dim, int spacedim>
286 bool
287 spherical_face_is_horizontal(
289 const Point<spacedim> &manifold_center)
290 {
291 // We test whether a face is horizontal by checking that the vertices
292 // all have roughly the same distance from the center: If the
293 // maximum deviation for the distances from the vertices to the
294 // center is less than 1.e-5 of the distance between vertices (as
295 // measured by the minimum distance from any of the other vertices
296 // to the first vertex), then we call this a horizontal face.
297 constexpr unsigned int n_vertices =
299 std::array<double, n_vertices> sqr_distances_to_center;
300 std::array<double, n_vertices - 1> sqr_distances_to_first_vertex;
301 sqr_distances_to_center[0] =
302 (face->vertex(0) - manifold_center).norm_square();
303 for (unsigned int i = 1; i < n_vertices; ++i)
304 {
305 sqr_distances_to_center[i] =
306 (face->vertex(i) - manifold_center).norm_square();
307 sqr_distances_to_first_vertex[i - 1] =
308 (face->vertex(i) - face->vertex(0)).norm_square();
309 }
310 const auto minmax_sqr_distance =
311 std::minmax_element(sqr_distances_to_center.begin(),
312 sqr_distances_to_center.end());
313 const auto min_sqr_distance_to_first_vertex =
314 std::min_element(sqr_distances_to_first_vertex.begin(),
315 sqr_distances_to_first_vertex.end());
316
317 return (*minmax_sqr_distance.second - *minmax_sqr_distance.first <
318 1.e-10 * *min_sqr_distance_to_first_vertex);
319 }
320} // namespace
321
322
323
324template <int dim, int spacedim>
328 const Point<spacedim> &p) const
329{
330 // Let us first test whether we are on a "horizontal" face
331 // (tangential to the sphere). In this case, the normal vector is
332 // easy to compute since it is proportional to the vector from the
333 // center to the point 'p'.
334 if (spherical_face_is_horizontal<dim, spacedim>(face, center))
335 {
336 // So, if this is a "horizontal" face, then just compute the normal
337 // vector as the one from the center to the point 'p', adequately
338 // scaled.
339 const Tensor<1, spacedim> unnormalized_spherical_normal = p - center;
340 const Tensor<1, spacedim> normalized_spherical_normal =
341 unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
342 return normalized_spherical_normal;
343 }
344 else
345 // If it is not a horizontal face, just use the machinery of the
346 // base class.
348
349 return Tensor<1, spacedim>();
350}
351
352
353
354// ============================================================
355// SphericalManifold
356// ============================================================
357
358template <int dim, int spacedim>
361 : center(center)
362 , p_center(center)
363 , polar_manifold(center)
364{}
365
366
367
368template <int dim, int spacedim>
369std::unique_ptr<Manifold<dim, spacedim>>
371{
372 return std::make_unique<SphericalManifold<dim, spacedim>>(*this);
373}
374
375
376
377template <int dim, int spacedim>
380 const Point<spacedim> &p1,
381 const Point<spacedim> &p2,
382 const double w) const
383{
384 const double tol = 1e-10;
385
386 if ((p1 - p2).norm_square() < tol * tol || std::abs(w) < tol)
387 return p1;
388 else if (std::abs(w - 1.0) < tol)
389 return p2;
390
391 // If the points are one dimensional then there is no need for anything but
392 // a linear combination.
393 if (spacedim == 1)
394 return Point<spacedim>(w * p2 + (1 - w) * p1);
395
396 const Tensor<1, spacedim> v1 = p1 - p_center;
397 const Tensor<1, spacedim> v2 = p2 - p_center;
398 const double r1 = v1.norm();
399 const double r2 = v2.norm();
400
401 Assert(r1 > tol && r2 > tol,
402 ExcMessage("p1 and p2 cannot coincide with the center."));
403
404 const Tensor<1, spacedim> e1 = v1 / r1;
405 const Tensor<1, spacedim> e2 = v2 / r2;
406
407 // Find the cosine of the angle gamma described by v1 and v2.
408 const double cosgamma = e1 * e2;
409
410 // Points are collinear with the center (allow for 8*eps as a tolerance)
411 if (cosgamma < -1 + 8. * std::numeric_limits<double>::epsilon())
412 return p_center;
413
414 // Points are along a line, in which case e1 and e2 are essentially the same.
415 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
416 return Point<spacedim>(p_center + w * v2 + (1 - w) * v1);
417
418 // Find the angle sigma that corresponds to arclength equal to w. acos
419 // should never be undefined because we have ruled out the two special cases
420 // above.
421 const double sigma = w * std::acos(cosgamma);
422
423 // Normal to v1 in the plane described by v1,v2,and the origin.
424 // Since p1 and p2 do not coincide n is not zero and well defined.
425 Tensor<1, spacedim> n = v2 - (v2 * e1) * e1;
426 const double n_norm = n.norm();
427 Assert(n_norm > 0,
428 ExcInternalError("n should be different from the null vector. "
429 "Probably, this means v1==v2 or v2==0."));
430
431 n /= n_norm;
432
433 // Find the point Q along O,v1 such that
434 // P1,V,P2 has measure sigma.
435 const Tensor<1, spacedim> P = std::cos(sigma) * e1 + std::sin(sigma) * n;
436
437 // Project this point on the manifold.
438 return Point<spacedim>(p_center + (w * r2 + (1.0 - w) * r1) * P);
439}
440
441
442
443template <int dim, int spacedim>
446 const Point<spacedim> &p1,
447 const Point<spacedim> &p2) const
448{
449 [[maybe_unused]] const double tol = 1e-10;
450
451 Assert(p1 != p2, ExcMessage("p1 and p2 should not concide."));
452
453 const Tensor<1, spacedim> v1 = p1 - p_center;
454 const Tensor<1, spacedim> v2 = p2 - p_center;
455 const double r1 = v1.norm();
456 const double r2 = v2.norm();
457
458 Assert(r1 > tol, ExcMessage("p1 cannot coincide with the center."));
459
460 Assert(r2 > tol, ExcMessage("p2 cannot coincide with the center."));
461
462 const Tensor<1, spacedim> e1 = v1 / r1;
463 const Tensor<1, spacedim> e2 = v2 / r2;
464
465 // Find the cosine of the angle gamma described by v1 and v2.
466 const double cosgamma = e1 * e2;
467
468 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
469 ExcMessage("p1 and p2 cannot lie on the same diameter and be opposite "
470 "respect to the center."));
471
472 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
473 return v2 - v1;
474
475 // Normal to v1 in the plane described by v1,v2,and the origin.
476 // Since p1 and p2 do not coincide n is not zero and well defined.
477 Tensor<1, spacedim> n = v2 - (v2 * e1) * e1;
478 const double n_norm = n.norm();
479 Assert(n_norm > 0,
480 ExcInternalError("n should be different from the null vector. "
481 "Probably, this means v1==v2 or v2==0."));
482
483 n /= n_norm;
484
485 // this is the derivative of the geodesic in get_intermediate_point
486 // derived with respect to w and inserting w=0.
487 const double gamma = std::acos(cosgamma);
488 return (r2 - r1) * e1 + r1 * gamma * n;
489}
490
491
492
493template <int dim, int spacedim>
497 const Point<spacedim> &p) const
498{
499 // Let us first test whether we are on a "horizontal" face
500 // (tangential to the sphere). In this case, the normal vector is
501 // easy to compute since it is proportional to the vector from the
502 // center to the point 'p'.
503 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
504 {
505 // So, if this is a "horizontal" face, then just compute the normal
506 // vector as the one from the center to the point 'p', adequately
507 // scaled.
508 const Tensor<1, spacedim> unnormalized_spherical_normal = p - p_center;
509 const Tensor<1, spacedim> normalized_spherical_normal =
510 unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
511 return normalized_spherical_normal;
512 }
513 else
514 // If it is not a horizontal face, just use the machinery of the
515 // base class.
517
518 return Tensor<1, spacedim>();
519}
520
521
522
523template <>
524void
531
532
533
534template <>
535void
542
543
544
545template <int dim, int spacedim>
546void
549 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
550 const
551{
552 // Let us first test whether we are on a "horizontal" face
553 // (tangential to the sphere). In this case, the normal vector is
554 // easy to compute since it is proportional to the vector from the
555 // center to the point 'p'.
556 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
557 {
558 // So, if this is a "horizontal" face, then just compute the normal
559 // vector as the one from the center to the point 'p', adequately
560 // scaled.
561 for (unsigned int vertex = 0;
562 vertex < GeometryInfo<spacedim>::vertices_per_face;
563 ++vertex)
564 face_vertex_normals[vertex] = face->vertex(vertex) - p_center;
565 }
566 else
567 Manifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
568}
569
570
571
572template <int dim, int spacedim>
573void
575 const ArrayView<const Point<spacedim>> &surrounding_points,
576 const Table<2, double> &weights,
577 ArrayView<Point<spacedim>> new_points) const
578{
579 AssertDimension(new_points.size(), weights.size(0));
580 AssertDimension(surrounding_points.size(), weights.size(1));
581
582 do_get_new_points(surrounding_points, make_array_view(weights), new_points);
583
584 return;
585}
586
587
588
589template <int dim, int spacedim>
592 const ArrayView<const Point<spacedim>> &vertices,
593 const ArrayView<const double> &weights) const
594{
595 // To avoid duplicating all of the logic in get_new_points, simply call it
596 // for one position.
597 Point<spacedim> new_point;
598 do_get_new_points(vertices,
599 weights,
600 make_array_view(&new_point, &new_point + 1));
601
602 return new_point;
603}
604
605
606
607namespace internal
608{
610 {
611 namespace
612 {
613 template <int spacedim>
615 do_get_new_point(
616 const ArrayView<const Tensor<1, spacedim>> & /*directions*/,
617 const ArrayView<const double> & /*distances*/,
618 const ArrayView<const double> & /*weights*/,
619 const Point<spacedim> & /*candidate_point*/)
620 {
622 return {};
623 }
624
625 template <>
627 do_get_new_point(
628 const ArrayView<const Tensor<1, 3>> &directions,
629 [[maybe_unused]] const ArrayView<const double> &distances,
630 const ArrayView<const double> &weights,
631 const Point<3> &candidate_point)
632 {
633 AssertDimension(directions.size(), distances.size());
634 AssertDimension(directions.size(), weights.size());
635
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;
640
641 {
642 // If the candidate happens to coincide with a normalized
643 // direction, we return it. Otherwise, the Hessian would be singular.
644 for (unsigned int i = 0; i < n_merged_points; ++i)
645 {
646 const double squared_distance =
647 (candidate - directions[i]).norm_square();
648 if (squared_distance < tolerance * tolerance)
649 return candidate;
650 }
651
652 // check if we only have two points now, in which case we can use the
653 // get_intermediate_point function
654 if (n_merged_points == 2)
655 {
656 static const ::SphericalManifold<3, 3> unit_manifold;
657 Assert(std::abs(weights[0] + weights[1] - 1.0) < 1e-13,
658 ExcMessage("Weights do not sum up to 1"));
659 const Point<3> intermediate =
660 unit_manifold.get_intermediate_point(Point<3>(directions[0]),
661 Point<3>(directions[1]),
662 weights[1]);
663 return intermediate;
664 }
665
666 Tensor<1, 3> vPerp;
667 Tensor<2, 2> Hessian;
669 Tensor<1, 2> gradlocal;
670
671 // On success we exit the loop early.
672 // Otherwise, we just take the result after max_iterations steps.
673 for (unsigned int i = 0; i < max_iterations; ++i)
674 {
675 // Step 2a: Find new descent direction
676
677 // Get local basis for the estimate candidate
678 const Tensor<1, 3> Clocalx = internal::compute_normal(candidate);
679 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
680
681 // For each vertices vector, compute the tangent vector from
682 // candidate towards the vertices vector -- its length is the
683 // spherical length from candidate to the vertices vector. Then
684 // compute its contribution to the Hessian.
685 gradient = 0.;
686 Hessian = 0.;
687 for (unsigned int i = 0; i < n_merged_points; ++i)
688 if (std::abs(weights[i]) > 1.e-15)
689 {
690 vPerp =
691 internal::projected_direction(directions[i], candidate);
692 const double sinthetaSq = vPerp.norm_square();
693 const double sintheta = std::sqrt(sinthetaSq);
694 if (sintheta < tolerance)
695 {
696 Hessian[0][0] += weights[i];
697 Hessian[1][1] += weights[i];
698 }
699 else
700 {
701 const double costheta = (directions[i]) * candidate;
702 const double theta = std::atan2(sintheta, costheta);
703 const double sincthetaInv = theta / sintheta;
704
705 const double cosphi = vPerp * Clocalx;
706 const double sinphi = vPerp * Clocaly;
707
708 gradlocal[0] = cosphi;
709 gradlocal[1] = sinphi;
710 gradient += (weights[i] * sincthetaInv) * gradlocal;
711
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);
722 }
723 }
724
725 Assert(determinant(Hessian) > tolerance, ExcInternalError());
726
727 const Tensor<2, 2> inverse_Hessian = invert(Hessian);
728
729 const Tensor<1, 2> xDisplocal = inverse_Hessian * gradient;
730 const Tensor<1, 3> xDisp =
731 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
732
733 // Step 2b: rotate candidate in direction xDisp for a new
734 // candidate.
735 const Point<3> candidateOld = candidate;
736 candidate =
737 Point<3>(internal::apply_exponential_map(candidate, xDisp));
738
739 // Step 2c: return the new candidate if we didn't move
740 if ((candidate - candidateOld).norm_square() <
741 tolerance * tolerance)
742 break;
743 }
744 }
745 return candidate;
746 }
747 } // namespace
748 } // namespace SphericalManifold
749} // namespace internal
750
751
752
753template <int dim, int spacedim>
754void
756 const ArrayView<const Point<spacedim>> &surrounding_points,
757 const ArrayView<const double> &weights,
758 ArrayView<Point<spacedim>> new_points) const
759{
760 AssertDimension(weights.size(),
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();
764
765 if (surrounding_points.size() == 2)
766 {
767 for (unsigned int row = 0; row < weight_rows; ++row)
768 new_points[row] =
770 surrounding_points[0],
771 surrounding_points[1],
772 weights[row * weight_columns + 1]);
773 return;
774 }
775
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(
779 surrounding_points.size(), Point<spacedim>());
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)
784 {
785 directions[i] = surrounding_points[i] - p_center;
786 distances[i] = directions[i].norm();
787
788 if (distances[i] != 0.)
789 directions[i] /= distances[i];
790 else
791 Assert(false,
792 ExcMessage("One of the vertices coincides with the center. "
793 "This is not allowed!"));
794
795 // Check if an estimate is good enough,
796 // this is often the case for sufficiently refined meshes.
797 for (unsigned int k = 0; k < i; ++k)
798 {
799 const double squared_distance =
800 (directions[i] - directions[k]).norm_square();
801 max_distance = std::max(max_distance, squared_distance);
802 }
803 }
804
805 // Step 1: Check for some special cases, create simple linear guesses
806 // otherwise.
807 const double tolerance = 1e-10;
808 boost::container::small_vector<bool, 100> accurate_point_was_found(
809 new_points.size(), false);
810 const ArrayView<const Tensor<1, spacedim>> array_directions =
811 make_array_view(directions.begin(), directions.end());
812 const ArrayView<const double> array_distances =
813 make_array_view(distances.begin(), distances.end());
814 for (unsigned int row = 0; row < weight_rows; ++row)
815 {
816 new_candidates[row] =
817 guess_new_point(array_directions,
818 array_distances,
819 ArrayView<const double>(&weights[row * weight_columns],
820 weight_columns));
821
822 // If the candidate is the center, mark it as found to avoid entering
823 // the Newton iteration in step 2, which would crash.
824 if (new_candidates[row].first == 0.0)
825 {
826 new_points[row] = p_center;
827 accurate_point_was_found[row] = true;
828 continue;
829 }
830
831 // If not in 3d, just use the implementation from PolarManifold
832 // after we verified that the candidate is not the center.
833 if (spacedim < 3)
834 new_points[row] = polar_manifold.get_new_point(
835 surrounding_points,
836 ArrayView<const double>(&weights[row * weight_columns],
837 weight_columns));
838 }
839
840 // In this case, we treated the case that the candidate is the center and
841 // obtained the new locations from the PolarManifold object otherwise.
842 if constexpr (spacedim < 3)
843 return;
844 else
845 {
846 // If all the points are close to each other, we expect the estimate to
847 // be good enough. This tolerance was chosen such that the first iteration
848 // for a at least three time refined HyperShell mesh with radii .5 and 1.
849 // doesn't already succeed.
850 if (max_distance < 2e-2)
851 {
852 for (unsigned int row = 0; row < weight_rows; ++row)
853 new_points[row] =
854 p_center + new_candidates[row].first * new_candidates[row].second;
855
856 return;
857 }
858
859 // Step 2:
860 // Do more expensive Newton-style iterations to improve the estimate.
861
862 // Search for duplicate directions and merge them to minimize the cost of
863 // the get_new_point function call below.
864 boost::container::small_vector<double, 1000> merged_weights(
865 weights.size());
866 boost::container::small_vector<Tensor<1, spacedim>, 100>
867 merged_directions(surrounding_points.size(), Point<spacedim>());
868 boost::container::small_vector<double, 100> merged_distances(
869 surrounding_points.size(), 0.0);
870
871 unsigned int n_unique_directions = 0;
872 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
873 {
874 bool found_duplicate = false;
875
876 // This inner loop is of @f$O(N^2)@f$ complexity, but
877 // surrounding_points.size() is usually at most 8 points large.
878 for (unsigned int j = 0; j < n_unique_directions; ++j)
879 {
880 const double squared_distance =
881 (directions[i] - directions[j]).norm_square();
882 if (!found_duplicate && squared_distance < 1e-28)
883 {
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];
888 }
889 }
890
891 if (found_duplicate == false)
892 {
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];
898
899 ++n_unique_directions;
900 }
901 }
902
903 // Search for duplicate weight rows and merge them to minimize the cost of
904 // the get_new_point function call below.
905 boost::container::small_vector<unsigned int, 100> merged_weights_index(
906 new_points.size(), numbers::invalid_unsigned_int);
907 for (unsigned int row = 0; row < weight_rows; ++row)
908 {
909 for (unsigned int existing_row = 0; existing_row < row;
910 ++existing_row)
911 {
912 bool identical_weights = true;
913
914 for (unsigned int weight_index = 0;
915 weight_index < n_unique_directions;
916 ++weight_index)
917 if (std::abs(
918 merged_weights[row * weight_columns + weight_index] -
919 merged_weights[existing_row * weight_columns +
920 weight_index]) > tolerance)
921 {
922 identical_weights = false;
923 break;
924 }
925
926 if (identical_weights)
927 {
928 merged_weights_index[row] = existing_row;
929 break;
930 }
931 }
932 }
933
934 // Note that we only use the n_unique_directions first entries in the
935 // ArrayView
936 const ArrayView<const Tensor<1, spacedim>> array_merged_directions =
937 make_array_view(merged_directions.begin(),
938 merged_directions.begin() + n_unique_directions);
939 const ArrayView<const double> array_merged_distances =
940 make_array_view(merged_distances.begin(),
941 merged_distances.begin() + n_unique_directions);
942
943 for (unsigned int row = 0; row < weight_rows; ++row)
944 if (!accurate_point_was_found[row])
945 {
946 if (merged_weights_index[row] == numbers::invalid_unsigned_int)
947 {
948 const ArrayView<const double> array_merged_weights(
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,
955 Point<spacedim>(new_candidates[row].second));
956 }
957 else
958 new_candidates[row].second =
959 new_candidates[merged_weights_index[row]].second;
960
961 new_points[row] =
962 p_center + new_candidates[row].first * new_candidates[row].second;
963 }
964 }
965}
966
967
968
969template <int dim, int spacedim>
970std::pair<double, Tensor<1, spacedim>>
972 const ArrayView<const Tensor<1, spacedim>> &directions,
973 const ArrayView<const double> &distances,
974 const ArrayView<const double> &weights) const
975{
976 const double tolerance = 1e-10;
977 double rho = 0.;
978 Tensor<1, spacedim> candidate;
979
980 // Perform a simple average ...
981 double total_weights = 0.;
982 for (unsigned int i = 0; i < directions.size(); ++i)
983 {
984 // if one weight is one, return its direction
985 if (std::abs(1 - weights[i]) < tolerance)
986 return std::make_pair(distances[i], directions[i]);
987
988 rho += distances[i] * weights[i];
989 candidate += directions[i] * weights[i];
990 total_weights += weights[i];
991 }
992
993 // ... and normalize if the candidate is different from the origin.
994 const double norm = candidate.norm();
995 if (norm == 0.)
996 return std::make_pair(0.0, Point<spacedim>());
997 candidate /= norm;
998 rho /= total_weights;
999
1000 return std::make_pair(rho, candidate);
1001}
1002
1003
1004
1005// ============================================================
1006// CylindricalManifold
1007// ============================================================
1008template <int dim, int spacedim>
1010 const double tolerance)
1011 : CylindricalManifold<dim, spacedim>(Point<spacedim>::unit_vector(axis),
1012 Point<spacedim>(),
1013 tolerance)
1014{
1015 // do not use static_assert to make dimension-independent programming
1016 // easier.
1017 Assert(spacedim == 3,
1018 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1019}
1020
1021
1022
1023template <int dim, int spacedim>
1025 const Tensor<1, spacedim> &direction,
1026 const Point<spacedim> &point_on_axis,
1027 const double tolerance)
1028 : ChartManifold<dim, spacedim, 3>(Tensor<1, 3>({0, 2. * numbers::PI, 0}))
1029 , normal_direction(internal::compute_normal(direction, true))
1030 , direction(direction / direction.norm())
1031 , point_on_axis(point_on_axis)
1032 , tolerance(tolerance)
1033 , dxn(cross_product_3d(this->direction, normal_direction))
1034{
1035 // do not use static_assert to make dimension-independent programming
1036 // easier.
1037 Assert(spacedim == 3,
1038 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1039}
1040
1041
1042
1043template <int dim, int spacedim>
1044std::unique_ptr<Manifold<dim, spacedim>>
1046{
1047 return std::make_unique<CylindricalManifold<dim, spacedim>>(*this);
1048}
1049
1050
1051
1052template <int dim, int spacedim>
1055 const ArrayView<const Point<spacedim>> &surrounding_points,
1056 const ArrayView<const double> &weights) const
1057{
1058 Assert(spacedim == 3,
1059 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1060
1061 // First check if the average in space lies on the axis.
1062 Point<spacedim> middle;
1063 double average_length = 0.;
1064 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
1065 {
1066 middle += surrounding_points[i] * weights[i];
1067 average_length += surrounding_points[i].square() * weights[i];
1068 }
1069 middle -= point_on_axis;
1070 const double lambda = middle * direction;
1071
1072 if ((middle - direction * lambda).square() < tolerance * average_length)
1073 return point_on_axis + direction * lambda;
1074 else // If not, using the ChartManifold should yield valid results.
1075 return ChartManifold<dim, spacedim, 3>::get_new_point(surrounding_points,
1076 weights);
1077}
1078
1079
1080
1081template <int dim, int spacedim>
1084 const Point<spacedim> &space_point) const
1085{
1086 Assert(spacedim == 3,
1087 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1088
1089 // First find the projection of the given point to the axis.
1090 const Tensor<1, spacedim> normalized_point = space_point - point_on_axis;
1091 const double lambda = normalized_point * direction;
1092 const Point<spacedim> projection = point_on_axis + direction * lambda;
1093 const Tensor<1, spacedim> p_diff = space_point - projection;
1094
1095 // Then compute the angle between the projection direction and
1096 // another vector orthogonal to the direction vector.
1097 const double phi = Physics::VectorRelations::signed_angle(normal_direction,
1098 p_diff,
1099 /*axis=*/direction);
1100
1101 // Return distance from the axis, angle and signed distance on the axis.
1102 return Point<3>(p_diff.norm(), phi, lambda);
1103}
1104
1105
1106
1107template <int dim, int spacedim>
1110 const Point<3> &chart_point) const
1111{
1112 Assert(spacedim == 3,
1113 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1114
1115 // Rotate the orthogonal direction by the given angle
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];
1118 const Tensor<1, spacedim> intermediate =
1119 normal_direction * cosine_r + dxn * sine_r;
1120
1121 // Finally, put everything together.
1122 return point_on_axis + direction * chart_point[2] + intermediate;
1123}
1124
1125
1126
1127template <int dim, int spacedim>
1130 const Point<3> &chart_point) const
1131{
1132 Assert(spacedim == 3,
1133 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1134
1135 Tensor<2, 3> derivatives;
1136
1137 // Rotate the orthogonal direction by the given angle
1138 const double sine = std::sin(chart_point[1]);
1139 const double cosine = std::cos(chart_point[1]);
1140 const Tensor<1, spacedim> intermediate =
1141 normal_direction * cosine + dxn * sine;
1142
1143 // avoid compiler warnings
1144 constexpr int s0 = 0 % spacedim;
1145 constexpr int s1 = 1 % spacedim;
1146 constexpr int s2 = 2 % spacedim;
1147
1148 // derivative w.r.t the radius
1149 derivatives[s0][s0] = intermediate[s0];
1150 derivatives[s1][s0] = intermediate[s1];
1151 derivatives[s2][s0] = intermediate[s2];
1152
1153 // derivatives w.r.t the angle
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;
1157
1158 // derivatives w.r.t the direction of the axis
1159 derivatives[s0][s2] = direction[s0];
1160 derivatives[s1][s2] = direction[s1];
1161 derivatives[s2][s2] = direction[s2];
1162
1163 return derivatives;
1164}
1165
1166
1167
1168namespace
1169{
1170 template <int dim>
1172 check_and_normalize(const Tensor<1, dim> &t)
1173 {
1174 const double norm = t.norm();
1175 Assert(norm > 0.0, ExcMessage("The major axis must have a positive norm."));
1176 return t / norm;
1177 }
1178} // namespace
1179
1180
1181
1182// ============================================================
1183// EllipticalManifold
1184// ============================================================
1185template <int dim, int spacedim>
1187 const Point<spacedim> &center,
1188 const Tensor<1, spacedim> &major_axis_direction,
1189 const double eccentricity)
1190 : ChartManifold<dim, spacedim, spacedim>(
1191 EllipticalManifold<dim, spacedim>::get_periodicity())
1192 , direction(check_and_normalize(major_axis_direction))
1193 , center(center)
1194 , eccentricity(eccentricity)
1195 , cosh_u(1.0 / eccentricity)
1196 , sinh_u(std::sqrt(cosh_u * cosh_u - 1.0))
1197{
1198 // Throws an exception if dim!=2 || spacedim!=2.
1199 Assert(dim == 2 && spacedim == 2, ExcNotImplemented());
1200 // Throws an exception if eccentricity is not in range.
1201 Assert(std::signbit(cosh_u * cosh_u - 1.0) == false,
1202 ExcMessage(
1203 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1204}
1205
1206
1207
1208template <int dim, int spacedim>
1209std::unique_ptr<Manifold<dim, spacedim>>
1211{
1212 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1213}
1214
1215
1216
1217template <int dim, int spacedim>
1220{
1221 Tensor<1, spacedim> periodicity;
1222 // The second elliptical coordinate is periodic, while the first is not.
1223 // Enforce periodicity on the last variable.
1224 periodicity[spacedim - 1] = 2.0 * numbers::PI;
1225 return periodicity;
1226}
1227
1228
1229
1230template <int dim, int spacedim>
1233{
1235 return {};
1236}
1237
1238
1239
1240template <>
1243{
1244 const double cs = std::cos(chart_point[1]);
1245 const double sn = std::sin(chart_point[1]);
1246 // Coordinates in the reference frame (i.e. major axis direction is
1247 // x-axis)
1248 const double x = chart_point[0] * cosh_u * cs;
1249 const double y = chart_point[0] * sinh_u * sn;
1250 // Rotate them according to the major axis direction
1251 const Point<2> p(direction[0] * x - direction[1] * y,
1252 direction[1] * x + direction[0] * y);
1253 return p + center;
1254}
1255
1256
1257
1258template <int dim, int spacedim>
1261{
1263 return {};
1264}
1265
1266
1267
1268template <>
1271{
1272 // Moving space_point in the reference coordinate system.
1273 const double x0 = space_point[0] - center[0];
1274 const double y0 = space_point[1] - center[1];
1275 const double x = direction[0] * x0 + direction[1] * y0;
1276 const double y = -direction[1] * x0 + direction[0] * y0;
1277 const double pt0 =
1278 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1279 // If the radius is exactly zero, the point coincides with the origin.
1280 if (pt0 == 0.0)
1281 {
1282 return center;
1283 }
1284 double cos_eta = x / (pt0 * cosh_u);
1285 if (cos_eta < -1.0)
1286 {
1287 cos_eta = -1.0;
1288 }
1289 if (cos_eta > 1.0)
1290 {
1291 cos_eta = 1.0;
1292 }
1293 const double eta = std::acos(cos_eta);
1294 const double pt1 = (std::signbit(y) ? 2.0 * numbers::PI - eta : eta);
1295 return {pt0, pt1};
1296}
1297
1298
1299
1300template <int dim, int spacedim>
1303 const Point<spacedim> &) const
1304{
1306 return {};
1307}
1308
1309
1310
1311template <>
1314 const Point<2> &chart_point) const
1315{
1316 const double cs = std::cos(chart_point[1]);
1317 const double sn = std::sin(chart_point[1]);
1318 Tensor<2, 2> dX;
1319 dX[0][0] = cosh_u * cs;
1320 dX[0][1] = -chart_point[0] * cosh_u * sn;
1321 dX[1][0] = sinh_u * sn;
1322 dX[1][1] = chart_point[0] * sinh_u * cs;
1323
1324 // rotate according to the major axis direction
1326 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1327
1328 return rot * dX;
1329}
1330
1331
1332
1333// ============================================================
1334// FunctionManifold
1335// ============================================================
1336template <int dim, int spacedim, int chartdim>
1338 const Function<chartdim> &push_forward_function,
1339 const Function<spacedim> &pull_back_function,
1340 const Tensor<1, chartdim> &periodicity,
1341 const double tolerance)
1342 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1343 , const_map()
1344 , push_forward_function(&push_forward_function)
1345 , pull_back_function(&pull_back_function)
1346 , tolerance(tolerance)
1347 , owns_pointers(false)
1348 , finite_difference_step(0)
1349{
1350 AssertDimension(push_forward_function.n_components, spacedim);
1351 AssertDimension(pull_back_function.n_components, chartdim);
1352}
1353
1354
1355
1356template <int dim, int spacedim, int chartdim>
1358 std::unique_ptr<Function<chartdim>> push_forward,
1359 std::unique_ptr<Function<spacedim>> pull_back,
1360 const Tensor<1, chartdim> &periodicity,
1361 const double tolerance)
1362 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1363 , const_map()
1364 , push_forward_function(push_forward.release())
1365 , pull_back_function(pull_back.release())
1366 , tolerance(tolerance)
1367 , owns_pointers(true)
1368 , finite_difference_step(0)
1369{
1370 AssertDimension(push_forward_function->n_components, spacedim);
1371 AssertDimension(pull_back_function->n_components, chartdim);
1372}
1373
1374
1375
1376template <int dim, int spacedim, int chartdim>
1378 const std::string push_forward_expression,
1379 const std::string pull_back_expression,
1380 const Tensor<1, chartdim> &periodicity,
1381 const typename FunctionParser<spacedim>::ConstMap const_map,
1382 const std::string chart_vars,
1383 const std::string space_vars,
1384 const double tolerance,
1385 const double h)
1386 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1387 , const_map(const_map)
1388 , tolerance(tolerance)
1389 , owns_pointers(true)
1390 , push_forward_expression(push_forward_expression)
1391 , pull_back_expression(pull_back_expression)
1392 , chart_vars(chart_vars)
1393 , space_vars(space_vars)
1394 , finite_difference_step(h)
1395{
1396 FunctionParser<chartdim> *pf = new FunctionParser<chartdim>(spacedim, 0.0, h);
1397 FunctionParser<spacedim> *pb = new FunctionParser<spacedim>(chartdim, 0.0, h);
1401 pull_back_function = pb;
1402}
1403
1404
1405
1406template <int dim, int spacedim, int chartdim>
1408{
1409 if (owns_pointers == true)
1410 {
1411 const Function<chartdim> *pf = push_forward_function;
1412 push_forward_function = nullptr;
1413 delete pf;
1414
1415 const Function<spacedim> *pb = pull_back_function;
1416 pull_back_function = nullptr;
1417 delete pb;
1418 }
1419}
1420
1421
1422
1423template <int dim, int spacedim, int chartdim>
1424std::unique_ptr<Manifold<dim, spacedim>>
1426{
1427 // This manifold can be constructed either by providing an expression for the
1428 // push forward and the pull back charts, or by providing two Function
1429 // objects. In the first case, the push_forward and pull_back functions are
1430 // created internally in FunctionManifold, and destroyed when this object is
1431 // deleted. In the second case, the function objects are destroyed if they
1432 // are passed as pointers upon construction.
1433 // We need to make sure that our cloned object is constructed in the
1434 // same way this class was constructed, and that its internal Function
1435 // pointers point either to the same Function objects used to construct this
1436 // function or that the newly generated manifold creates internally the
1437 // push_forward and pull_back functions using the same expressions that were
1438 // used to construct this class.
1439 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1440 {
1441 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1442 push_forward_expression,
1443 pull_back_expression,
1444 this->get_periodicity(),
1445 const_map,
1446 chart_vars,
1447 space_vars,
1448 tolerance,
1449 finite_difference_step);
1450 }
1451 else
1452 {
1453 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1454 *push_forward_function,
1455 *pull_back_function,
1456 this->get_periodicity(),
1457 tolerance);
1458 }
1459}
1460
1461
1462
1463template <int dim, int spacedim, int chartdim>
1466 const Point<chartdim> &chart_point) const
1467{
1468 Vector<double> pf(spacedim);
1469 Point<spacedim> result;
1470 push_forward_function->vector_value(chart_point, pf);
1471 for (unsigned int i = 0; i < spacedim; ++i)
1472 result[i] = pf[i];
1473
1474#ifdef DEBUG
1475 Vector<double> pb(chartdim);
1476 pull_back_function->vector_value(result, pb);
1477 for (unsigned int i = 0; i < chartdim; ++i)
1478 Assert(
1479 (chart_point.norm() > tolerance &&
1480 (std::abs(pb[i] - chart_point[i]) < tolerance * chart_point.norm())) ||
1481 (std::abs(pb[i] - chart_point[i]) < tolerance),
1482 ExcMessage(
1483 "The push forward is not the inverse of the pull back! Bailing out."));
1484#endif
1485
1486 return result;
1487}
1488
1489
1490
1491template <int dim, int spacedim, int chartdim>
1494 const Point<chartdim> &chart_point) const
1495{
1497 for (unsigned int i = 0; i < spacedim; ++i)
1498 {
1499 const auto gradient = push_forward_function->gradient(chart_point, i);
1500 for (unsigned int j = 0; j < chartdim; ++j)
1501 DF[i][j] = gradient[j];
1502 }
1503 return DF;
1504}
1505
1506
1507
1508template <int dim, int spacedim, int chartdim>
1511 const Point<spacedim> &space_point) const
1512{
1513 Vector<double> pb(chartdim);
1514 Point<chartdim> result;
1515 pull_back_function->vector_value(space_point, pb);
1516 for (unsigned int i = 0; i < chartdim; ++i)
1517 result[i] = pb[i];
1518 return result;
1519}
1520
1521
1522
1523// ============================================================
1524// TorusManifold
1525// ============================================================
1526template <int dim>
1529{
1530 double x = p[0];
1531 double z = p[1];
1532 double y = p[2];
1533 double phi = std::atan2(y, x);
1534 double theta = std::atan2(z, std::sqrt(x * x + y * y) - centerline_radius);
1535 double w =
1536 std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * centerline_radius) +
1537 Utilities::fixed_power<2>(x - std::cos(phi) * centerline_radius) +
1538 z * z) /
1540 return {phi, theta, w};
1541}
1542
1543
1544
1545template <int dim>
1548{
1549 double phi = chart_point[0];
1550 double theta = chart_point[1];
1551 double w = chart_point[2];
1552
1553 return {std::cos(phi) * centerline_radius +
1554 inner_radius * w * std::cos(theta) * std::cos(phi),
1555 inner_radius * w * std::sin(theta),
1556 std::sin(phi) * centerline_radius +
1557 inner_radius * w * std::cos(theta) * std::sin(phi)};
1558}
1559
1560
1561
1562template <int dim>
1564 const double inner_radius)
1565 : ChartManifold<dim, 3, 3>(Point<3>(2 * numbers::PI, 2 * numbers::PI, 0.0))
1568{
1570 ExcMessage("The centerline radius must be greater than the "
1572 Assert(inner_radius > 0.0, ExcMessage("The inner radius must be positive."));
1573}
1574
1575
1576
1577template <int dim>
1578std::unique_ptr<Manifold<dim, 3>>
1580{
1582}
1583
1584
1585
1586template <int dim>
1589{
1591
1592 double phi = chart_point[0];
1593 double theta = chart_point[1];
1594 double w = chart_point[2];
1595
1596 DX[0][0] = -std::sin(phi) * centerline_radius -
1597 inner_radius * w * std::cos(theta) * std::sin(phi);
1598 DX[0][1] = -inner_radius * w * std::sin(theta) * std::cos(phi);
1599 DX[0][2] = inner_radius * std::cos(theta) * std::cos(phi);
1600
1601 DX[1][0] = 0;
1602 DX[1][1] = inner_radius * w * std::cos(theta);
1603 DX[1][2] = inner_radius * std::sin(theta);
1604
1605 DX[2][0] = std::cos(phi) * centerline_radius +
1606 inner_radius * w * std::cos(theta) * std::cos(phi);
1607 DX[2][1] = -inner_radius * w * std::sin(theta) * std::sin(phi);
1608 DX[2][2] = inner_radius * std::cos(theta) * std::sin(phi);
1609
1610 return DX;
1611}
1612
1613
1614
1615// ============================================================
1616// TransfiniteInterpolationManifold
1617// ============================================================
1618template <int dim, int spacedim>
1621 : triangulation(nullptr)
1622 , level_coarse(-1)
1623{
1624 AssertThrow(dim > 1, ExcNotImplemented());
1625}
1626
1627
1628
1629template <int dim, int spacedim>
1631 spacedim>::~TransfiniteInterpolationManifold()
1632{
1633 if (clear_signal.connected())
1634 clear_signal.disconnect();
1635}
1636
1637
1638
1639template <int dim, int spacedim>
1640std::unique_ptr<Manifold<dim, spacedim>>
1642{
1644 if (triangulation)
1645 ptr->initialize(*triangulation);
1646 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1647}
1648
1649
1650
1651template <int dim, int spacedim>
1652void
1655{
1656 this->triangulation = &triangulation;
1657 // In case the triangulation is cleared, remove the pointers by a signal:
1658 clear_signal.disconnect();
1659 clear_signal = triangulation.signals.clear.connect([&]() -> void {
1660 this->triangulation = nullptr;
1661 this->level_coarse = -1;
1662 });
1663 level_coarse = triangulation.last()->level();
1664 coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false);
1666
1667 // In case of dim == spacedim we perform a quadratic approximation in
1668 // InverseQuadraticApproximation(), thus initialize the unit_points
1669 // vector with one subdivision to get 3^dim unit_points.
1670 //
1671 // In the co-dimension one case (meaning dim < spacedim) we have to fall
1672 // back to a simple GridTools::affine_cell_approximation<dim>() which
1673 // requires 2^dim points, instead. Thus, initialize the QIterated
1674 // quadrature with no subdivisions.
1675 const std::vector<Point<dim>> unit_points =
1676 QIterated<dim>(QTrapezoid<1>(), (dim == spacedim ? 2 : 1)).get_points();
1677 std::vector<Point<spacedim>> real_points(unit_points.size());
1678
1679 for (const auto &cell : triangulation.active_cell_iterators())
1680 {
1681 bool cell_is_flat = true;
1682 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1683 if (cell->line(l)->manifold_id() != cell->manifold_id() &&
1684 cell->line(l)->manifold_id() != numbers::flat_manifold_id)
1685 cell_is_flat = false;
1686 if (dim > 2)
1687 for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1688 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1689 cell->quad(q)->manifold_id() != numbers::flat_manifold_id)
1690 cell_is_flat = false;
1691 AssertIndexRange(static_cast<unsigned int>(cell->index()),
1692 coarse_cell_is_flat.size());
1693 coarse_cell_is_flat[cell->index()] = cell_is_flat;
1694
1695 // build quadratic interpolation
1696 for (unsigned int i = 0; i < unit_points.size(); ++i)
1697 real_points[i] = push_forward(cell, unit_points[i]);
1699 }
1700}
1701
1702
1703
1704namespace
1705{
1706 // version for 1d
1707 template <typename AccessorType>
1709 compute_transfinite_interpolation(const AccessorType &cell,
1710 const Point<1> &chart_point,
1711 const bool /*cell_is_flat*/)
1712 {
1713 return cell.vertex(0) * (1. - chart_point[0]) +
1714 cell.vertex(1) * chart_point[0];
1715 }
1716
1717 // version for 2d
1718 template <typename AccessorType>
1720 compute_transfinite_interpolation(const AccessorType &cell,
1721 const Point<2> &chart_point,
1722 const bool cell_is_flat)
1723 {
1724 const unsigned int dim = AccessorType::dimension;
1725 const unsigned int spacedim = AccessorType::space_dimension;
1726 const types::manifold_id my_manifold_id = cell.manifold_id();
1728
1729 // formula see wikipedia
1730 // https://en.wikipedia.org/wiki/Transfinite_interpolation
1731 // S(u,v) = (1-v)c_1(u)+v c_3(u) + (1-u)c_2(v) + u c_4(v) -
1732 // [(1-u)(1-v)P_0 + u(1-v) P_1 + (1-u)v P_2 + uv P_3]
1733 const std::array<Point<spacedim>, 4> vertices{
1734 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1735
1736 // this evaluates all bilinear shape functions because we need them
1737 // repeatedly. we will update this values in the complicated case with
1738 // curved lines below
1739 std::array<double, 4> weights_vertices{
1740 {(1. - chart_point[0]) * (1. - chart_point[1]),
1741 chart_point[0] * (1. - chart_point[1]),
1742 (1. - chart_point[0]) * chart_point[1],
1743 chart_point[0] * chart_point[1]}};
1744
1745 Point<spacedim> new_point;
1746 if (cell_is_flat)
1747 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1748 new_point += weights_vertices[v] * vertices[v];
1749 else
1750 {
1751 // The second line in the formula tells us to subtract the
1752 // contribution of the vertices. If a line employs the same manifold
1753 // as the cell, we can merge the weights of the line with the weights
1754 // of the vertex with a negative sign while going through the faces
1755 // (this is a bit artificial in 2d but it becomes clear in 3d where we
1756 // avoid looking at the faces' orientation and other complications).
1757
1758 // add the contribution from the lines around the cell (first line in
1759 // formula)
1760 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1761 std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_face> points;
1762 // note that the views are immutable, but the arrays are not
1763 const auto weights_view =
1764 make_array_view(weights.begin(), weights.end());
1765 const auto points_view = make_array_view(points.begin(), points.end());
1766
1767 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1768 ++line)
1769 {
1770 const double my_weight =
1771 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1772 const double line_point = chart_point[1 - line / 2];
1773
1774 // Same manifold or invalid id which will go back to the same
1775 // class -> contribution should be added for the final point,
1776 // which means that we subtract the current weight from the
1777 // negative weight applied to the vertex
1778 const types::manifold_id line_manifold_id =
1779 cell.line(line)->manifold_id();
1780 if (line_manifold_id == my_manifold_id ||
1781 line_manifold_id == numbers::flat_manifold_id)
1782 {
1783 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,
1784 0)] -=
1785 my_weight * (1. - line_point);
1786 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,
1787 1)] -=
1788 my_weight * line_point;
1789 }
1790 else
1791 {
1792 points[0] =
1794 points[1] =
1796 weights[0] = 1. - line_point;
1797 weights[1] = line_point;
1798 new_point +=
1799 my_weight * tria.get_manifold(line_manifold_id)
1800 .get_new_point(points_view, weights_view);
1801 }
1802 }
1803
1804 // subtract contribution from the vertices (second line in formula)
1805 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1806 new_point -= weights_vertices[v] * vertices[v];
1807 }
1808
1809 return new_point;
1810 }
1811
1812 // this is replicated from GeometryInfo::face_to_cell_vertices since we need
1813 // it very often in compute_transfinite_interpolation and the function is
1814 // performance critical
1815 static constexpr unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1816 {1, 3, 5, 7},
1817 {0, 4, 1, 5},
1818 {2, 6, 3, 7},
1819 {0, 1, 2, 3},
1820 {4, 5, 6, 7}};
1821
1822 // this is replicated from GeometryInfo::face_to_cell_lines since we need it
1823 // very often in compute_transfinite_interpolation and the function is
1824 // performance critical
1825 static constexpr unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1826 {9, 11, 1, 5},
1827 {2, 6, 8, 9},
1828 {3, 7, 10, 11},
1829 {0, 1, 2, 3},
1830 {4, 5, 6, 7}};
1831
1832 // version for 3d
1833 template <typename AccessorType>
1835 compute_transfinite_interpolation(const AccessorType &cell,
1836 const Point<3> &chart_point,
1837 const bool cell_is_flat)
1838 {
1839 const unsigned int dim = AccessorType::dimension;
1840 const unsigned int spacedim = AccessorType::space_dimension;
1841 const types::manifold_id my_manifold_id = cell.manifold_id();
1843
1844 // Same approach as in 2d, but adding the faces, subtracting the edges, and
1845 // adding the vertices
1846 const std::array<Point<spacedim>, 8> vertices{{cell.vertex(0),
1847 cell.vertex(1),
1848 cell.vertex(2),
1849 cell.vertex(3),
1850 cell.vertex(4),
1851 cell.vertex(5),
1852 cell.vertex(6),
1853 cell.vertex(7)}};
1854
1855 // store the components of the linear shape functions because we need them
1856 // repeatedly. we allow for 10 such shape functions to wrap around the
1857 // first four once again for easier face access.
1858 double linear_shapes[10];
1859 for (unsigned int d = 0; d < 3; ++d)
1860 {
1861 linear_shapes[2 * d] = 1. - chart_point[d];
1862 linear_shapes[2 * d + 1] = chart_point[d];
1863 }
1864
1865 // wrap linear shape functions around for access in face loop
1866 for (unsigned int d = 6; d < 10; ++d)
1867 linear_shapes[d] = linear_shapes[d - 6];
1868
1869 std::array<double, 8> weights_vertices;
1870 for (unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1871 for (unsigned int i1 = 0; i1 < 2; ++i1)
1872 for (unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1873 weights_vertices[v] =
1874 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1875
1876 Point<spacedim> new_point;
1877 if (cell_is_flat)
1878 for (unsigned int v = 0; v < 8; ++v)
1879 new_point += weights_vertices[v] * vertices[v];
1880 else
1881 {
1882 // identify the weights for the lines to be accumulated (vertex
1883 // weights are set outside and coincide with the flat manifold case)
1884
1885 std::array<double, GeometryInfo<3>::lines_per_cell> weights_lines;
1886 std::fill(weights_lines.begin(), weights_lines.end(), 0.0);
1887
1888 // start with the contributions of the faces
1889 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1890 std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> points;
1891 // note that the views are immutable, but the arrays are not
1892 const auto weights_view =
1893 make_array_view(weights.begin(), weights.end());
1894 const auto points_view = make_array_view(points.begin(), points.end());
1895
1896 for (const unsigned int face : GeometryInfo<3>::face_indices())
1897 {
1898 const double my_weight = linear_shapes[face];
1899 const unsigned int face_even = face - face % 2;
1900
1901 if (std::abs(my_weight) < 1e-13)
1902 continue;
1903
1904 // same manifold or invalid id which will go back to the same class
1905 // -> face will interpolate from the surrounding lines and vertices
1906 const types::manifold_id face_manifold_id =
1907 cell.face(face)->manifold_id();
1908 if (face_manifold_id == my_manifold_id ||
1909 face_manifold_id == numbers::flat_manifold_id)
1910 {
1911 for (unsigned int line = 0;
1912 line < GeometryInfo<2>::lines_per_cell;
1913 ++line)
1914 {
1915 const double line_weight =
1916 linear_shapes[face_even + 2 + line];
1917 weights_lines[face_to_cell_lines_3d[face][line]] +=
1918 my_weight * line_weight;
1919 }
1920 // as to the indices inside linear_shapes: we use the index
1921 // wrapped around at 2*d, ensuring the correct orientation of
1922 // the face's coordinate system with respect to the
1923 // lexicographic indices
1924 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1925 linear_shapes[face_even + 2] *
1926 (linear_shapes[face_even + 4] * my_weight);
1927 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1928 linear_shapes[face_even + 3] *
1929 (linear_shapes[face_even + 4] * my_weight);
1930 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1931 linear_shapes[face_even + 2] *
1932 (linear_shapes[face_even + 5] * my_weight);
1933 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1934 linear_shapes[face_even + 3] *
1935 (linear_shapes[face_even + 5] * my_weight);
1936 }
1937 else
1938 {
1939 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1940 points[v] = vertices[face_to_cell_vertices_3d[face][v]];
1941 weights[0] =
1942 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1943 weights[1] =
1944 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1945 weights[2] =
1946 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1947 weights[3] =
1948 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1949 new_point +=
1950 my_weight * tria.get_manifold(face_manifold_id)
1951 .get_new_point(points_view, weights_view);
1952 }
1953 }
1954
1955 // next subtract the contributions of the lines
1956 const auto weights_view_line =
1957 make_array_view(weights.begin(), weights.begin() + 2);
1958 const auto points_view_line =
1959 make_array_view(points.begin(), points.begin() + 2);
1960 for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1961 ++line)
1962 {
1963 const double line_point =
1964 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1965 double my_weight = 0.;
1966 if (line < 8)
1967 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1968 else
1969 {
1970 const unsigned int subline = line - 8;
1971 my_weight =
1972 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1973 }
1974 my_weight -= weights_lines[line];
1975
1976 if (std::abs(my_weight) < 1e-13)
1977 continue;
1978
1979 const types::manifold_id line_manifold_id =
1980 cell.line(line)->manifold_id();
1981 if (line_manifold_id == my_manifold_id ||
1982 line_manifold_id == numbers::flat_manifold_id)
1983 {
1984 weights_vertices[GeometryInfo<3>::line_to_cell_vertices(line,
1985 0)] -=
1986 my_weight * (1. - line_point);
1987 weights_vertices[GeometryInfo<3>::line_to_cell_vertices(line,
1988 1)] -=
1989 my_weight * (line_point);
1990 }
1991 else
1992 {
1993 points[0] =
1995 points[1] =
1997 weights[0] = 1. - line_point;
1998 weights[1] = line_point;
1999 new_point -= my_weight * tria.get_manifold(line_manifold_id)
2000 .get_new_point(points_view_line,
2001 weights_view_line);
2002 }
2003 }
2004
2005 // finally add the contribution of the
2006 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2007 new_point += weights_vertices[v] * vertices[v];
2008 }
2009 return new_point;
2010 }
2011} // namespace
2012
2013
2014
2015template <int dim, int spacedim>
2019 const Point<dim> &chart_point) const
2020{
2021 AssertDimension(cell->level(), level_coarse);
2022
2023 // check that the point is in the unit cell which is the current chart
2024 // Tolerance 5e-4 chosen that the method also works with manifolds
2025 // that have some discretization error like SphericalManifold
2027 ExcMessage("chart_point is not in unit interval"));
2028
2029 return compute_transfinite_interpolation(*cell,
2030 chart_point,
2031 coarse_cell_is_flat[cell->index()]);
2032}
2033
2034
2035
2036template <int dim, int spacedim>
2040 const Point<dim> &chart_point,
2041 const Point<spacedim> &pushed_forward_chart_point) const
2042{
2043 // compute the derivative with the help of finite differences
2045 for (unsigned int d = 0; d < dim; ++d)
2046 {
2047 Point<dim> modified = chart_point;
2048 const double step = chart_point[d] > 0.5 ? -1e-8 : 1e-8;
2049
2050 // avoid checking outside of the unit interval
2051 modified[d] += step;
2052 Tensor<1, spacedim> difference =
2053 compute_transfinite_interpolation(*cell,
2054 modified,
2055 coarse_cell_is_flat[cell->index()]) -
2056 pushed_forward_chart_point;
2057 for (unsigned int e = 0; e < spacedim; ++e)
2058 grad[e][d] = difference[e] / step;
2059 }
2061}
2062
2063
2064
2065template <int dim, int spacedim>
2069 const Point<spacedim> &point,
2070 const Point<dim> &initial_guess) const
2071{
2072 Point<dim> outside;
2073 for (unsigned int d = 0; d < dim; ++d)
2075
2076 // project the user-given input to unit cell
2077 Point<dim> chart_point = cell->reference_cell().closest_point(initial_guess);
2078
2079 // run quasi-Newton iteration with a combination of finite differences for
2080 // the exact Jacobian and "Broyden's good method". As opposed to the various
2081 // mapping implementations, this class does not throw exception upon failure
2082 // as those are relatively expensive and failure occurs quite regularly in
2083 // the implementation of the compute_chart_points method.
2084 Tensor<1, spacedim> residual =
2085 point -
2086 compute_transfinite_interpolation(*cell,
2087 chart_point,
2088 coarse_cell_is_flat[cell->index()]);
2089 const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter());
2090 double residual_norm_square = residual.norm_square();
2092 bool must_recompute_jacobian = true;
2093 for (unsigned int i = 0; i < 100; ++i)
2094 {
2095 if (residual_norm_square < tolerance)
2096 {
2097 // do a final update of the point with the last available Jacobian
2098 // information. The residual is close to zero due to the check
2099 // above, but me might improve some of the last digits by a final
2100 // Newton-like step with step length 1
2101 Tensor<1, dim> update;
2102 for (unsigned int d = 0; d < spacedim; ++d)
2103 for (unsigned int e = 0; e < dim; ++e)
2104 update[e] += inv_grad[d][e] * residual[d];
2105 return chart_point + update;
2106 }
2107
2108 // every 9 iterations, including the first time around, we create an
2109 // approximation of the Jacobian with finite differences. Broyden's
2110 // method usually does not need more than 5-8 iterations, but sometimes
2111 // we might have had a bad initial guess and then we can accelerate
2112 // convergence considerably with getting the actual Jacobian rather than
2113 // using secant-like methods (one gradient calculation in 3d costs as
2114 // much as 3 more iterations). this usually happens close to convergence
2115 // and one more step with the finite-differenced Jacobian leads to
2116 // convergence
2117 if (must_recompute_jacobian || i % 9 == 0)
2118 {
2119 // if the determinant is zero or negative, the mapping is either not
2120 // invertible or already has inverted and we are outside the valid
2121 // chart region. Note that the Jacobian here represents the
2122 // derivative of the forward map and should have a positive
2123 // determinant since we use properly oriented meshes.
2126 chart_point,
2127 Point<spacedim>(point - residual));
2128 if (grad.determinant() <= 0.0)
2129 return outside;
2131 must_recompute_jacobian = false;
2132 }
2133 Tensor<1, dim> update;
2134 for (unsigned int d = 0; d < spacedim; ++d)
2135 for (unsigned int e = 0; e < dim; ++e)
2136 update[e] += inv_grad[d][e] * residual[d];
2137
2138 // Line search, accept step if the residual has decreased
2139 double alpha = 1.;
2140
2141 // check if point is inside 1.2 times the unit cell to avoid
2142 // hitting points very far away from valid ones in the manifolds
2143 while (
2144 !GeometryInfo<dim>::is_inside_unit_cell(chart_point + alpha * update,
2145 0.2) &&
2146 alpha > 1e-7)
2147 alpha *= 0.5;
2148
2149 const Tensor<1, spacedim> old_residual = residual;
2150 while (alpha > 1e-4)
2151 {
2152 Point<dim> guess = chart_point + alpha * update;
2153 const Tensor<1, spacedim> residual_guess =
2154 point - compute_transfinite_interpolation(
2155 *cell, guess, coarse_cell_is_flat[cell->index()]);
2156 const double residual_norm_new = residual_guess.norm_square();
2157 if (residual_norm_new < residual_norm_square)
2158 {
2159 residual = residual_guess;
2160 residual_norm_square = residual_norm_new;
2161 chart_point += alpha * update;
2162 break;
2163 }
2164 else
2165 alpha *= 0.5;
2166 }
2167 // If alpha got very small, it is likely due to a bad Jacobian
2168 // approximation with Broyden's method (relatively far away from the
2169 // zero), which can be corrected by the outer loop when a Newton update
2170 // is recomputed. The second case is when the Jacobian is actually bad
2171 // and we should fail as early as possible. Since we cannot really
2172 // distinguish the two, we must continue here in any case.
2173 if (alpha <= 1e-4)
2174 must_recompute_jacobian = true;
2175
2176 // update the inverse Jacobian with "Broyden's good method" and
2177 // Sherman-Morrison formula for the update of the inverse, see
2178 // https://en.wikipedia.org/wiki/Broyden%27s_method
2179 // J^{-1}_n = J^{-1}_{n-1} + (delta x_n - J^{-1}_{n-1} delta f_n) /
2180 // (delta x_n^T J_{-1}_{n-1} delta f_n) delta x_n^T J^{-1}_{n-1}
2181
2182 // switch sign in residual as compared to the formula above because we
2183 // use a negative definition of the residual with respect to the
2184 // Jacobian
2185 const Tensor<1, spacedim> delta_f = old_residual - residual;
2186
2187 Tensor<1, dim> Jinv_deltaf;
2188 for (unsigned int d = 0; d < spacedim; ++d)
2189 for (unsigned int e = 0; e < dim; ++e)
2190 Jinv_deltaf[e] += inv_grad[d][e] * delta_f[d];
2191
2192 const Tensor<1, dim> delta_x = alpha * update;
2193
2194 // prevent division by zero. This number should be scale-invariant
2195 // because Jinv_deltaf carries no units and x is in reference
2196 // coordinates.
2197 if (std::abs(delta_x * Jinv_deltaf) > 1e-12 && !must_recompute_jacobian)
2198 {
2199 const Tensor<1, dim> factor =
2200 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
2201 Tensor<1, spacedim> jac_update;
2202 for (unsigned int d = 0; d < spacedim; ++d)
2203 for (unsigned int e = 0; e < dim; ++e)
2204 jac_update[d] += delta_x[e] * inv_grad[d][e];
2205 for (unsigned int d = 0; d < spacedim; ++d)
2206 for (unsigned int e = 0; e < dim; ++e)
2207 inv_grad[d][e] += factor[e] * jac_update[d];
2208 }
2209 }
2210 return outside;
2211}
2212
2213
2214
2215template <int dim, int spacedim>
2216std::array<unsigned int, 20>
2219 const ArrayView<const Point<spacedim>> &points) const
2220{
2221 // The methods to identify cells around points in GridTools are all written
2222 // for the active cells, but we are here looking at some cells at the coarse
2223 // level.
2224 Assert(triangulation != nullptr, ExcNotInitialized());
2225 Assert(triangulation->begin_active()->level() >= level_coarse,
2226 ExcMessage("The manifold was initialized with level " +
2227 std::to_string(level_coarse) + " but there are now" +
2228 "active cells on a lower level. Coarsening the mesh is " +
2229 "currently not supported"));
2230
2231 // This computes the distance of the surrounding points transformed to the
2232 // unit cell from the unit cell.
2234 triangulation->begin(
2235 level_coarse),
2236 endc =
2237 triangulation->end(
2238 level_coarse);
2239 boost::container::small_vector<std::pair<double, unsigned int>, 200>
2240 distances_and_cells;
2241 for (; cell != endc; ++cell)
2242 {
2243 // only consider cells where the current manifold is attached
2244 if (&cell->get_manifold() != this)
2245 continue;
2246
2247 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
2248 vertices;
2249 for (const unsigned int vertex_n : GeometryInfo<dim>::vertex_indices())
2250 {
2251 vertices[vertex_n] = cell->vertex(vertex_n);
2252 }
2253
2254 // cheap check: if any of the points is not inside a circle around the
2255 // center of the loop, we can skip the expensive part below (this assumes
2256 // that the manifold does not deform the grid too much)
2258 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2259 center += vertices[v];
2261 double radius_square = 0.;
2262 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2264 std::max(radius_square, (center - vertices[v]).norm_square());
2265 bool inside_circle = true;
2266 for (unsigned int i = 0; i < points.size(); ++i)
2267 if ((center - points[i]).norm_square() > radius_square * 1.5)
2268 {
2269 inside_circle = false;
2270 break;
2271 }
2272 if (inside_circle == false)
2273 continue;
2274
2275 // slightly more expensive search
2276 double current_distance = 0;
2277 for (unsigned int i = 0; i < points.size(); ++i)
2278 {
2279 Point<dim> point =
2281 current_distance += GeometryInfo<dim>::distance_to_unit_cell(point);
2282 }
2283 distances_and_cells.push_back(
2284 std::make_pair(current_distance, cell->index()));
2285 }
2286 // no coarse cell could be found -> transformation failed
2287 AssertThrow(distances_and_cells.size() > 0,
2289 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2290 std::array<unsigned int, 20> cells;
2292 for (unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2293 ++i)
2294 cells[i] = distances_and_cells[i].second;
2295
2296 return cells;
2297}
2298
2299
2300
2301template <int dim, int spacedim>
2304 const ArrayView<const Point<spacedim>> &surrounding_points,
2305 ArrayView<Point<dim>> chart_points) const
2306{
2307 Assert(surrounding_points.size() == chart_points.size(),
2308 ExcMessage("The chart points array view must be as large as the "
2309 "surrounding points array view."));
2310
2311 std::array<unsigned int, 20> nearby_cells =
2312 get_possible_cells_around_points(surrounding_points);
2313
2314 // This function is nearly always called to place new points on a cell or
2315 // cell face. In this case, the general structure of the surrounding points
2316 // is known (i.e., if there are eight surrounding points, then they will
2317 // almost surely be either eight points around a quadrilateral or the eight
2318 // vertices of a cube). Hence, making this assumption, we use two
2319 // optimizations (one for structdim == 2 and one for structdim == 3) that
2320 // guess the locations of some of the chart points more efficiently than the
2321 // affine map approximation. The affine map approximation is used whenever
2322 // we don't have a cheaper guess available.
2323
2324 // Function that can guess the location of a chart point by assuming that
2325 // the eight surrounding points are points on a two-dimensional object
2326 // (either a cell in 2d or the face of a hexahedron in 3d), arranged like
2327 //
2328 // 2 - 7 - 3
2329 // | |
2330 // 4 5
2331 // | |
2332 // 0 - 6 - 1
2333 //
2334 // This function assumes that the first three chart points have been
2335 // computed since there is no effective way to guess them.
2336 auto guess_chart_point_structdim_2 = [&](const unsigned int i) -> Point<dim> {
2337 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2338 ExcMessage("This function assumes that there are eight surrounding "
2339 "points around a two-dimensional object. It also assumes "
2340 "that the first three chart points have already been "
2341 "computed."));
2342 switch (i)
2343 {
2344 case 0:
2345 case 1:
2346 case 2:
2348 break;
2349 case 3:
2350 return chart_points[1] + (chart_points[2] - chart_points[0]);
2351 case 4:
2352 return 0.5 * (chart_points[0] + chart_points[2]);
2353 case 5:
2354 return 0.5 * (chart_points[1] + chart_points[3]);
2355 case 6:
2356 return 0.5 * (chart_points[0] + chart_points[1]);
2357 case 7:
2358 return 0.5 * (chart_points[2] + chart_points[3]);
2359 default:
2361 }
2362
2363 return {};
2364 };
2365
2366 // Function that can guess the location of a chart point by assuming that
2367 // the eight surrounding points form the vertices of a hexahedron, arranged
2368 // like
2369 //
2370 // 6-------7
2371 // /| /|
2372 // / / |
2373 // / | / |
2374 // 4-------5 |
2375 // | 2- -|- -3
2376 // | / | /
2377 // | | /
2378 // |/ |/
2379 // 0-------1
2380 //
2381 // (where vertex 2 is the back left vertex) we can estimate where chart
2382 // points 5 - 7 are by computing the height (in chart coordinates) as c4 -
2383 // c0 and then adding that onto the appropriate bottom vertex.
2384 //
2385 // This function assumes that the first five chart points have been computed
2386 // since there is no effective way to guess them.
2387 auto guess_chart_point_structdim_3 = [&](const unsigned int i) -> Point<dim> {
2388 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2389 ExcMessage("This function assumes that there are eight surrounding "
2390 "points around a three-dimensional object. It also "
2391 "assumes that the first five chart points have already "
2392 "been computed."));
2393 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2394 };
2395
2396 // Check if we can use the two chart point shortcuts above before we start:
2397 bool use_structdim_2_guesses = false;
2398 bool use_structdim_3_guesses = false;
2399 // note that in the structdim 2 case: 0 - 6 and 2 - 7 should be roughly
2400 // parallel, while in the structdim 3 case, 0 - 6 and 2 - 7 should be roughly
2401 // orthogonal. Use the angle between these two vectors to figure out if we
2402 // should turn on either structdim optimization.
2403 if (surrounding_points.size() == 8)
2404 {
2405 const Tensor<1, spacedim> v06 =
2406 surrounding_points[6] - surrounding_points[0];
2407 const Tensor<1, spacedim> v27 =
2408 surrounding_points[7] - surrounding_points[2];
2409
2410 // note that we can save a call to sqrt() by rearranging
2411 const double cosine = scalar_product(v06, v27) /
2412 std::sqrt(v06.norm_square() * v27.norm_square());
2413 if (0.707 < cosine)
2414 // the angle is less than pi/4, so these vectors are roughly parallel:
2415 // enable the structdim 2 optimization
2416 use_structdim_2_guesses = true;
2417 else if (spacedim == 3)
2418 // otherwise these vectors are roughly orthogonal: enable the
2419 // structdim 3 optimization if we are in 3d
2420 use_structdim_3_guesses = true;
2421 }
2422 // we should enable at most one of the optimizations
2423 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2424 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2426
2427
2428
2429 auto compute_chart_point =
2430 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2431 const unsigned int point_index) {
2432 Point<dim> guess;
2433 // an optimization: keep track of whether or not we used the quadratic
2434 // approximation so that we don't call pull_back with the same
2435 // initial guess twice (i.e., if pull_back fails the first time,
2436 // don't try again with the same function arguments).
2437 bool used_quadratic_approximation = false;
2438 // if we have already computed three points, we can guess the fourth
2439 // to be the missing corner point of a rectangle
2440 if (point_index == 3 && surrounding_points.size() >= 8)
2441 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2442 else if (use_structdim_2_guesses && 3 < point_index)
2443 guess = guess_chart_point_structdim_2(point_index);
2444 else if (use_structdim_3_guesses && 4 < point_index)
2445 guess = guess_chart_point_structdim_3(point_index);
2446 else if (dim == 3 && point_index > 7 && surrounding_points.size() == 26)
2447 {
2448 if (point_index < 20)
2449 guess =
2450 0.5 * (chart_points[GeometryInfo<dim>::line_to_cell_vertices(
2451 point_index - 8, 0)] +
2453 point_index - 8, 1)]);
2454 else
2455 guess =
2456 0.25 * (chart_points[GeometryInfo<dim>::face_to_cell_vertices(
2457 point_index - 20, 0)] +
2459 point_index - 20, 1)] +
2461 point_index - 20, 2)] +
2463 point_index - 20, 3)]);
2464 }
2465 else
2466 {
2467 guess = quadratic_approximation[cell->index()].compute(
2468 surrounding_points[point_index]);
2469 used_quadratic_approximation = true;
2470 }
2471 chart_points[point_index] =
2472 pull_back(cell, surrounding_points[point_index], guess);
2473
2474 // the initial guess may not have been good enough: if applicable,
2475 // try again with the affine approximation (which is more accurate
2476 // than the cheap methods used above)
2477 if (chart_points[point_index][0] ==
2480 {
2481 guess = quadratic_approximation[cell->index()].compute(
2482 surrounding_points[point_index]);
2483 chart_points[point_index] =
2484 pull_back(cell, surrounding_points[point_index], guess);
2485 }
2486
2487 if (chart_points[point_index][0] ==
2489 {
2490 for (unsigned int d = 0; d < dim; ++d)
2491 guess[d] = 0.5;
2492 chart_points[point_index] =
2493 pull_back(cell, surrounding_points[point_index], guess);
2494 }
2495 };
2496
2497 // check whether all points are inside the unit cell of the current chart
2498 for (unsigned int c = 0; c < nearby_cells.size(); ++c)
2499 {
2501 triangulation, level_coarse, nearby_cells[c]);
2502 bool inside_unit_cell = true;
2503 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
2504 {
2505 compute_chart_point(cell, i);
2506
2507 // Tolerance 5e-4 chosen that the method also works with manifolds
2508 // that have some discretization error like SphericalManifold
2509 if (GeometryInfo<dim>::is_inside_unit_cell(chart_points[i], 5e-4) ==
2510 false)
2511 {
2512 inside_unit_cell = false;
2513 break;
2514 }
2515 }
2516 if (inside_unit_cell == true)
2517 {
2518 return cell;
2519 }
2520
2521 // if we did not find a point and this was the last valid cell (the next
2522 // iterate being the end of the array or an invalid tag), we must stop
2523 if (c == nearby_cells.size() - 1 ||
2524 nearby_cells[c + 1] == numbers::invalid_unsigned_int)
2525 {
2526 // generate additional information to help debugging why we did not
2527 // get a point
2528 std::ostringstream message;
2529 for (unsigned int b = 0; b <= c; ++b)
2530 {
2532 triangulation, level_coarse, nearby_cells[b]);
2533 message << "Looking at cell " << cell->id()
2534 << " with vertices: " << std::endl;
2535 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2536 message << cell->vertex(v) << " ";
2537 message << std::endl;
2538 message << "Transformation to chart coordinates: " << std::endl;
2539 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
2540 {
2541 compute_chart_point(cell, i);
2542 message << surrounding_points[i] << " -> " << chart_points[i]
2543 << std::endl;
2544 }
2545 }
2546
2547 AssertThrow(false,
2549 message.str())));
2550 }
2551 }
2552
2553 // a valid inversion should have returned a point above. an invalid
2554 // inversion should have triggered the assertion, so we should never end up
2555 // here
2558}
2559
2560
2561
2562template <int dim, int spacedim>
2565 const ArrayView<const Point<spacedim>> &surrounding_points,
2566 const ArrayView<const double> &weights) const
2567{
2568 boost::container::small_vector<Point<dim>, 100> chart_points(
2569 surrounding_points.size());
2570 ArrayView<Point<dim>> chart_points_view =
2571 make_array_view(chart_points.begin(), chart_points.end());
2572 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2573
2574 const Point<dim> p_chart =
2575 chart_manifold.get_new_point(chart_points_view, weights);
2576
2577 return push_forward(cell, p_chart);
2578}
2579
2580
2581
2582template <int dim, int spacedim>
2583void
2585 const ArrayView<const Point<spacedim>> &surrounding_points,
2586 const Table<2, double> &weights,
2587 ArrayView<Point<spacedim>> new_points) const
2588{
2589 Assert(weights.size(0) > 0, ExcEmptyObject());
2590 AssertDimension(surrounding_points.size(), weights.size(1));
2591
2592 boost::container::small_vector<Point<dim>, 100> chart_points(
2593 surrounding_points.size());
2594 ArrayView<Point<dim>> chart_points_view =
2595 make_array_view(chart_points.begin(), chart_points.end());
2596 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2597
2598 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2599 weights.size(0));
2600 chart_manifold.get_new_points(chart_points_view,
2601 weights,
2602 make_array_view(new_points_on_chart.begin(),
2603 new_points_on_chart.end()));
2604
2605 for (unsigned int row = 0; row < weights.size(0); ++row)
2606 new_points[row] = push_forward(cell, new_points_on_chart[row]);
2607}
2608
2609
2610
2611// explicit instantiations
2612#include "manifold_lib.inst"
2613
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
std::size_t size() const
Definition array_view.h:684
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition manifold.cc:1027
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
Number determinant() const
DerivativeForm< 1, dim, spacedim, Number > covariant_form() const
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 > &center, const Tensor< 1, spacedim > &major_axis_direction, const double eccentricity)
virtual DerivativeForm< 1, spacedim, spacedim > push_forward_gradient(const Point< spacedim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
static Tensor< 1, spacedim > get_periodicity()
const double cosh_u
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
Definition manifold.h:306
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.
Definition mapping.h:318
Definition point.h:111
PolarManifold(const Point< spacedim > center=Point< spacedim >())
static Tensor< 1, spacedim > get_periodicity()
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual DerivativeForm< 1, spacedim, spacedim > push_forward_gradient(const Point< spacedim > &chart_point) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
virtual Point< spacedim > push_forward(const Point< spacedim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &vertices, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
std::pair< double, Tensor< 1, spacedim > > guess_new_point(const ArrayView< const Tensor< 1, spacedim > > &directions, const ArrayView< const double > &distances, const ArrayView< const double > &weights) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
SphericalManifold(const Point< spacedim > center=Point< spacedim >())
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
void do_get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights, ArrayView< Point< spacedim > > new_points) const
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
virtual Point< 3 > pull_back(const Point< 3 > &p) const override
virtual DerivativeForm< 1, 3, 3 > push_forward_gradient(const Point< 3 > &chart_point) const override
virtual Point< 3 > push_forward(const Point< 3 > &chart_point) const override
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()
Signals signals
Definition tria.h:2524
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:500
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:501
Point< 3 > center
Point< 3 > vertices[4]
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int vertex_indices[2]
const unsigned int v1
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
Definition tria.h:1556
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
Definition numbers.h:254
static const unsigned int invalid_unsigned_int
Definition types.h:220
const types::manifold_id flat_manifold_id
Definition types.h:325
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > acos(const ::VectorizedArray< Number, width > &x)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
static double distance_to_unit_cell(const Point< dim > &p)
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
static unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()
boost::signals2::signal< void()> clear
Definition tria.h:2385
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 > &)