Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
manifold_lib.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
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 const double tol = 1e-10;
450 (void)tol;
451
452 Assert(p1 != p2, ExcMessage("p1 and p2 should not concide."));
453
454 const Tensor<1, spacedim> v1 = p1 - p_center;
455 const Tensor<1, spacedim> v2 = p2 - p_center;
456 const double r1 = v1.norm();
457 const double r2 = v2.norm();
458
459 Assert(r1 > tol, ExcMessage("p1 cannot coincide with the center."));
460
461 Assert(r2 > tol, ExcMessage("p2 cannot coincide with the center."));
462
463 const Tensor<1, spacedim> e1 = v1 / r1;
464 const Tensor<1, spacedim> e2 = v2 / r2;
465
466 // Find the cosine of the angle gamma described by v1 and v2.
467 const double cosgamma = e1 * e2;
468
469 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
470 ExcMessage("p1 and p2 cannot lie on the same diameter and be opposite "
471 "respect to the center."));
472
473 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
474 return v2 - v1;
475
476 // Normal to v1 in the plane described by v1,v2,and the origin.
477 // Since p1 and p2 do not coincide n is not zero and well defined.
478 Tensor<1, spacedim> n = v2 - (v2 * e1) * e1;
479 const double n_norm = n.norm();
480 Assert(n_norm > 0,
481 ExcInternalError("n should be different from the null vector. "
482 "Probably, this means v1==v2 or v2==0."));
483
484 n /= n_norm;
485
486 // this is the derivative of the geodesic in get_intermediate_point
487 // derived with respect to w and inserting w=0.
488 const double gamma = std::acos(cosgamma);
489 return (r2 - r1) * e1 + r1 * gamma * n;
490}
491
492
493
494template <int dim, int spacedim>
498 const Point<spacedim> &p) const
499{
500 // Let us first test whether we are on a "horizontal" face
501 // (tangential to the sphere). In this case, the normal vector is
502 // easy to compute since it is proportional to the vector from the
503 // center to the point 'p'.
504 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
505 {
506 // So, if this is a "horizontal" face, then just compute the normal
507 // vector as the one from the center to the point 'p', adequately
508 // scaled.
509 const Tensor<1, spacedim> unnormalized_spherical_normal = p - p_center;
510 const Tensor<1, spacedim> normalized_spherical_normal =
511 unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
512 return normalized_spherical_normal;
513 }
514 else
515 // If it is not a horizontal face, just use the machinery of the
516 // base class.
518
519 return Tensor<1, spacedim>();
520}
521
522
523
524template <>
525void
532
533
534
535template <>
536void
543
544
545
546template <int dim, int spacedim>
547void
550 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
551 const
552{
553 // Let us first test whether we are on a "horizontal" face
554 // (tangential to the sphere). In this case, the normal vector is
555 // easy to compute since it is proportional to the vector from the
556 // center to the point 'p'.
557 if (spherical_face_is_horizontal<dim, spacedim>(face, p_center))
558 {
559 // So, if this is a "horizontal" face, then just compute the normal
560 // vector as the one from the center to the point 'p', adequately
561 // scaled.
562 for (unsigned int vertex = 0;
563 vertex < GeometryInfo<spacedim>::vertices_per_face;
564 ++vertex)
565 face_vertex_normals[vertex] = face->vertex(vertex) - p_center;
566 }
567 else
568 Manifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
569}
570
571
572
573template <int dim, int spacedim>
574void
576 const ArrayView<const Point<spacedim>> &surrounding_points,
577 const Table<2, double> &weights,
578 ArrayView<Point<spacedim>> new_points) const
579{
580 AssertDimension(new_points.size(), weights.size(0));
581 AssertDimension(surrounding_points.size(), weights.size(1));
582
583 do_get_new_points(surrounding_points, make_array_view(weights), new_points);
584
585 return;
586}
587
588
589
590template <int dim, int spacedim>
593 const ArrayView<const Point<spacedim>> &vertices,
594 const ArrayView<const double> &weights) const
595{
596 // To avoid duplicating all of the logic in get_new_points, simply call it
597 // for one position.
598 Point<spacedim> new_point;
599 do_get_new_points(vertices,
600 weights,
601 make_array_view(&new_point, &new_point + 1));
602
603 return new_point;
604}
605
606
607
608namespace internal
609{
611 {
612 namespace
613 {
614 template <int spacedim>
616 do_get_new_point(
617 const ArrayView<const Tensor<1, spacedim>> & /*directions*/,
618 const ArrayView<const double> & /*distances*/,
619 const ArrayView<const double> & /*weights*/,
620 const Point<spacedim> & /*candidate_point*/)
621 {
623 return {};
624 }
625
626 template <>
628 do_get_new_point(const ArrayView<const Tensor<1, 3>> &directions,
629 const ArrayView<const double> &distances,
630 const ArrayView<const double> &weights,
631 const Point<3> &candidate_point)
632 {
633 (void)distances;
634
635 AssertDimension(directions.size(), distances.size());
636 AssertDimension(directions.size(), weights.size());
637
638 Point<3> candidate = candidate_point;
639 const unsigned int n_merged_points = directions.size();
640 const double tolerance = 1e-10;
641 const int max_iterations = 10;
642
643 {
644 // If the candidate happens to coincide with a normalized
645 // direction, we return it. Otherwise, the Hessian would be singular.
646 for (unsigned int i = 0; i < n_merged_points; ++i)
647 {
648 const double squared_distance =
649 (candidate - directions[i]).norm_square();
650 if (squared_distance < tolerance * tolerance)
651 return candidate;
652 }
653
654 // check if we only have two points now, in which case we can use the
655 // get_intermediate_point function
656 if (n_merged_points == 2)
657 {
658 static const ::SphericalManifold<3, 3> unit_manifold;
659 Assert(std::abs(weights[0] + weights[1] - 1.0) < 1e-13,
660 ExcMessage("Weights do not sum up to 1"));
661 const Point<3> intermediate =
662 unit_manifold.get_intermediate_point(Point<3>(directions[0]),
663 Point<3>(directions[1]),
664 weights[1]);
665 return intermediate;
666 }
667
668 Tensor<1, 3> vPerp;
669 Tensor<2, 2> Hessian;
671 Tensor<1, 2> gradlocal;
672
673 // On success we exit the loop early.
674 // Otherwise, we just take the result after max_iterations steps.
675 for (unsigned int i = 0; i < max_iterations; ++i)
676 {
677 // Step 2a: Find new descent direction
678
679 // Get local basis for the estimate candidate
680 const Tensor<1, 3> Clocalx = internal::compute_normal(candidate);
681 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
682
683 // For each vertices vector, compute the tangent vector from
684 // candidate towards the vertices vector -- its length is the
685 // spherical length from candidate to the vertices vector. Then
686 // compute its contribution to the Hessian.
687 gradient = 0.;
688 Hessian = 0.;
689 for (unsigned int i = 0; i < n_merged_points; ++i)
690 if (std::abs(weights[i]) > 1.e-15)
691 {
692 vPerp =
693 internal::projected_direction(directions[i], candidate);
694 const double sinthetaSq = vPerp.norm_square();
695 const double sintheta = std::sqrt(sinthetaSq);
696 if (sintheta < tolerance)
697 {
698 Hessian[0][0] += weights[i];
699 Hessian[1][1] += weights[i];
700 }
701 else
702 {
703 const double costheta = (directions[i]) * candidate;
704 const double theta = std::atan2(sintheta, costheta);
705 const double sincthetaInv = theta / sintheta;
706
707 const double cosphi = vPerp * Clocalx;
708 const double sinphi = vPerp * Clocaly;
709
710 gradlocal[0] = cosphi;
711 gradlocal[1] = sinphi;
712 gradient += (weights[i] * sincthetaInv) * gradlocal;
713
714 const double wt = weights[i] / sinthetaSq;
715 const double sinphiSq = sinphi * sinphi;
716 const double cosphiSq = cosphi * cosphi;
717 const double tt = sincthetaInv * costheta;
718 const double offdiag =
719 cosphi * sinphi * wt * (1.0 - tt);
720 Hessian[0][0] += wt * (cosphiSq + tt * sinphiSq);
721 Hessian[0][1] += offdiag;
722 Hessian[1][0] += offdiag;
723 Hessian[1][1] += wt * (sinphiSq + tt * cosphiSq);
724 }
725 }
726
727 Assert(determinant(Hessian) > tolerance, ExcInternalError());
728
729 const Tensor<2, 2> inverse_Hessian = invert(Hessian);
730
731 const Tensor<1, 2> xDisplocal = inverse_Hessian * gradient;
732 const Tensor<1, 3> xDisp =
733 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
734
735 // Step 2b: rotate candidate in direction xDisp for a new
736 // candidate.
737 const Point<3> candidateOld = candidate;
738 candidate =
739 Point<3>(internal::apply_exponential_map(candidate, xDisp));
740
741 // Step 2c: return the new candidate if we didn't move
742 if ((candidate - candidateOld).norm_square() <
743 tolerance * tolerance)
744 break;
745 }
746 }
747 return candidate;
748 }
749 } // namespace
750 } // namespace SphericalManifold
751} // namespace internal
752
753
754
755template <int dim, int spacedim>
756void
758 const ArrayView<const Point<spacedim>> &surrounding_points,
759 const ArrayView<const double> &weights,
760 ArrayView<Point<spacedim>> new_points) const
761{
762 AssertDimension(weights.size(),
763 new_points.size() * surrounding_points.size());
764 const unsigned int weight_rows = new_points.size();
765 const unsigned int weight_columns = surrounding_points.size();
766
767 if (surrounding_points.size() == 2)
768 {
769 for (unsigned int row = 0; row < weight_rows; ++row)
770 new_points[row] =
772 surrounding_points[0],
773 surrounding_points[1],
774 weights[row * weight_columns + 1]);
775 return;
776 }
777
778 boost::container::small_vector<std::pair<double, Tensor<1, spacedim>>, 100>
779 new_candidates(new_points.size());
780 boost::container::small_vector<Tensor<1, spacedim>, 100> directions(
781 surrounding_points.size(), Point<spacedim>());
782 boost::container::small_vector<double, 100> distances(
783 surrounding_points.size(), 0.0);
784 double max_distance = 0.;
785 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
786 {
787 directions[i] = surrounding_points[i] - p_center;
788 distances[i] = directions[i].norm();
789
790 if (distances[i] != 0.)
791 directions[i] /= distances[i];
792 else
793 Assert(false,
794 ExcMessage("One of the vertices coincides with the center. "
795 "This is not allowed!"));
796
797 // Check if an estimate is good enough,
798 // this is often the case for sufficiently refined meshes.
799 for (unsigned int k = 0; k < i; ++k)
800 {
801 const double squared_distance =
802 (directions[i] - directions[k]).norm_square();
803 max_distance = std::max(max_distance, squared_distance);
804 }
805 }
806
807 // Step 1: Check for some special cases, create simple linear guesses
808 // otherwise.
809 const double tolerance = 1e-10;
810 boost::container::small_vector<bool, 100> accurate_point_was_found(
811 new_points.size(), false);
812 const ArrayView<const Tensor<1, spacedim>> array_directions =
813 make_array_view(directions.begin(), directions.end());
814 const ArrayView<const double> array_distances =
815 make_array_view(distances.begin(), distances.end());
816 for (unsigned int row = 0; row < weight_rows; ++row)
817 {
818 new_candidates[row] =
819 guess_new_point(array_directions,
820 array_distances,
821 ArrayView<const double>(&weights[row * weight_columns],
822 weight_columns));
823
824 // If the candidate is the center, mark it as found to avoid entering
825 // the Newton iteration in step 2, which would crash.
826 if (new_candidates[row].first == 0.0)
827 {
828 new_points[row] = p_center;
829 accurate_point_was_found[row] = true;
830 continue;
831 }
832
833 // If not in 3d, just use the implementation from PolarManifold
834 // after we verified that the candidate is not the center.
835 if (spacedim < 3)
836 new_points[row] = polar_manifold.get_new_point(
837 surrounding_points,
838 ArrayView<const double>(&weights[row * weight_columns],
839 weight_columns));
840 }
841
842 // In this case, we treated the case that the candidate is the center and
843 // obtained the new locations from the PolarManifold object otherwise.
844 if constexpr (spacedim < 3)
845 return;
846 else
847 {
848 // If all the points are close to each other, we expect the estimate to
849 // be good enough. This tolerance was chosen such that the first iteration
850 // for a at least three time refined HyperShell mesh with radii .5 and 1.
851 // doesn't already succeed.
852 if (max_distance < 2e-2)
853 {
854 for (unsigned int row = 0; row < weight_rows; ++row)
855 new_points[row] =
856 p_center + new_candidates[row].first * new_candidates[row].second;
857
858 return;
859 }
860
861 // Step 2:
862 // Do more expensive Newton-style iterations to improve the estimate.
863
864 // Search for duplicate directions and merge them to minimize the cost of
865 // the get_new_point function call below.
866 boost::container::small_vector<double, 1000> merged_weights(
867 weights.size());
868 boost::container::small_vector<Tensor<1, spacedim>, 100>
869 merged_directions(surrounding_points.size(), Point<spacedim>());
870 boost::container::small_vector<double, 100> merged_distances(
871 surrounding_points.size(), 0.0);
872
873 unsigned int n_unique_directions = 0;
874 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
875 {
876 bool found_duplicate = false;
877
878 // This inner loop is of @f$O(N^2)@f$ complexity, but
879 // surrounding_points.size() is usually at most 8 points large.
880 for (unsigned int j = 0; j < n_unique_directions; ++j)
881 {
882 const double squared_distance =
883 (directions[i] - directions[j]).norm_square();
884 if (!found_duplicate && squared_distance < 1e-28)
885 {
886 found_duplicate = true;
887 for (unsigned int row = 0; row < weight_rows; ++row)
888 merged_weights[row * weight_columns + j] +=
889 weights[row * weight_columns + i];
890 }
891 }
892
893 if (found_duplicate == false)
894 {
895 merged_directions[n_unique_directions] = directions[i];
896 merged_distances[n_unique_directions] = distances[i];
897 for (unsigned int row = 0; row < weight_rows; ++row)
898 merged_weights[row * weight_columns + n_unique_directions] =
899 weights[row * weight_columns + i];
900
901 ++n_unique_directions;
902 }
903 }
904
905 // Search for duplicate weight rows and merge them to minimize the cost of
906 // the get_new_point function call below.
907 boost::container::small_vector<unsigned int, 100> merged_weights_index(
908 new_points.size(), numbers::invalid_unsigned_int);
909 for (unsigned int row = 0; row < weight_rows; ++row)
910 {
911 for (unsigned int existing_row = 0; existing_row < row;
912 ++existing_row)
913 {
914 bool identical_weights = true;
915
916 for (unsigned int weight_index = 0;
917 weight_index < n_unique_directions;
918 ++weight_index)
919 if (std::abs(
920 merged_weights[row * weight_columns + weight_index] -
921 merged_weights[existing_row * weight_columns +
922 weight_index]) > tolerance)
923 {
924 identical_weights = false;
925 break;
926 }
927
928 if (identical_weights)
929 {
930 merged_weights_index[row] = existing_row;
931 break;
932 }
933 }
934 }
935
936 // Note that we only use the n_unique_directions first entries in the
937 // ArrayView
938 const ArrayView<const Tensor<1, spacedim>> array_merged_directions =
939 make_array_view(merged_directions.begin(),
940 merged_directions.begin() + n_unique_directions);
941 const ArrayView<const double> array_merged_distances =
942 make_array_view(merged_distances.begin(),
943 merged_distances.begin() + n_unique_directions);
944
945 for (unsigned int row = 0; row < weight_rows; ++row)
946 if (!accurate_point_was_found[row])
947 {
948 if (merged_weights_index[row] == numbers::invalid_unsigned_int)
949 {
950 const ArrayView<const double> array_merged_weights(
951 &merged_weights[row * weight_columns], n_unique_directions);
952 new_candidates[row].second =
953 internal::SphericalManifold::do_get_new_point(
954 array_merged_directions,
955 array_merged_distances,
956 array_merged_weights,
957 Point<spacedim>(new_candidates[row].second));
958 }
959 else
960 new_candidates[row].second =
961 new_candidates[merged_weights_index[row]].second;
962
963 new_points[row] =
964 p_center + new_candidates[row].first * new_candidates[row].second;
965 }
966 }
967}
968
969
970
971template <int dim, int spacedim>
972std::pair<double, Tensor<1, spacedim>>
974 const ArrayView<const Tensor<1, spacedim>> &directions,
975 const ArrayView<const double> &distances,
976 const ArrayView<const double> &weights) const
977{
978 const double tolerance = 1e-10;
979 double rho = 0.;
980 Tensor<1, spacedim> candidate;
981
982 // Perform a simple average ...
983 double total_weights = 0.;
984 for (unsigned int i = 0; i < directions.size(); ++i)
985 {
986 // if one weight is one, return its direction
987 if (std::abs(1 - weights[i]) < tolerance)
988 return std::make_pair(distances[i], directions[i]);
989
990 rho += distances[i] * weights[i];
991 candidate += directions[i] * weights[i];
992 total_weights += weights[i];
993 }
994
995 // ... and normalize if the candidate is different from the origin.
996 const double norm = candidate.norm();
997 if (norm == 0.)
998 return std::make_pair(0.0, Point<spacedim>());
999 candidate /= norm;
1000 rho /= total_weights;
1001
1002 return std::make_pair(rho, candidate);
1003}
1004
1005
1006
1007// ============================================================
1008// CylindricalManifold
1009// ============================================================
1010template <int dim, int spacedim>
1012 const double tolerance)
1013 : CylindricalManifold<dim, spacedim>(Point<spacedim>::unit_vector(axis),
1014 Point<spacedim>(),
1015 tolerance)
1016{
1017 // do not use static_assert to make dimension-independent programming
1018 // easier.
1019 Assert(spacedim == 3,
1020 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1021}
1022
1023
1024
1025template <int dim, int spacedim>
1027 const Tensor<1, spacedim> &direction,
1028 const Point<spacedim> &point_on_axis,
1029 const double tolerance)
1030 : ChartManifold<dim, spacedim, 3>(Tensor<1, 3>({0, 2. * numbers::PI, 0}))
1031 , normal_direction(internal::compute_normal(direction, true))
1032 , direction(direction / direction.norm())
1033 , point_on_axis(point_on_axis)
1034 , tolerance(tolerance)
1035 , dxn(cross_product_3d(this->direction, normal_direction))
1036{
1037 // do not use static_assert to make dimension-independent programming
1038 // easier.
1039 Assert(spacedim == 3,
1040 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1041}
1042
1043
1044
1045template <int dim, int spacedim>
1046std::unique_ptr<Manifold<dim, spacedim>>
1048{
1049 return std::make_unique<CylindricalManifold<dim, spacedim>>(*this);
1050}
1051
1052
1053
1054template <int dim, int spacedim>
1057 const ArrayView<const Point<spacedim>> &surrounding_points,
1058 const ArrayView<const double> &weights) const
1059{
1060 Assert(spacedim == 3,
1061 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1062
1063 // First check if the average in space lies on the axis.
1064 Point<spacedim> middle;
1065 double average_length = 0.;
1066 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
1067 {
1068 middle += surrounding_points[i] * weights[i];
1069 average_length += surrounding_points[i].square() * weights[i];
1070 }
1071 middle -= point_on_axis;
1072 const double lambda = middle * direction;
1073
1074 if ((middle - direction * lambda).square() < tolerance * average_length)
1075 return point_on_axis + direction * lambda;
1076 else // If not, using the ChartManifold should yield valid results.
1077 return ChartManifold<dim, spacedim, 3>::get_new_point(surrounding_points,
1078 weights);
1079}
1080
1081
1082
1083template <int dim, int spacedim>
1086 const Point<spacedim> &space_point) const
1087{
1088 Assert(spacedim == 3,
1089 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1090
1091 // First find the projection of the given point to the axis.
1092 const Tensor<1, spacedim> normalized_point = space_point - point_on_axis;
1093 const double lambda = normalized_point * direction;
1094 const Point<spacedim> projection = point_on_axis + direction * lambda;
1095 const Tensor<1, spacedim> p_diff = space_point - projection;
1096
1097 // Then compute the angle between the projection direction and
1098 // another vector orthogonal to the direction vector.
1099 const double phi = Physics::VectorRelations::signed_angle(normal_direction,
1100 p_diff,
1101 /*axis=*/direction);
1102
1103 // Return distance from the axis, angle and signed distance on the axis.
1104 return Point<3>(p_diff.norm(), phi, lambda);
1105}
1106
1107
1108
1109template <int dim, int spacedim>
1112 const Point<3> &chart_point) const
1113{
1114 Assert(spacedim == 3,
1115 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1116
1117 // Rotate the orthogonal direction by the given angle
1118 const double sine_r = std::sin(chart_point[1]) * chart_point[0];
1119 const double cosine_r = std::cos(chart_point[1]) * chart_point[0];
1120 const Tensor<1, spacedim> intermediate =
1121 normal_direction * cosine_r + dxn * sine_r;
1122
1123 // Finally, put everything together.
1124 return point_on_axis + direction * chart_point[2] + intermediate;
1125}
1126
1127
1128
1129template <int dim, int spacedim>
1132 const Point<3> &chart_point) const
1133{
1134 Assert(spacedim == 3,
1135 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1136
1137 Tensor<2, 3> derivatives;
1138
1139 // Rotate the orthogonal direction by the given angle
1140 const double sine = std::sin(chart_point[1]);
1141 const double cosine = std::cos(chart_point[1]);
1142 const Tensor<1, spacedim> intermediate =
1143 normal_direction * cosine + dxn * sine;
1144
1145 // avoid compiler warnings
1146 constexpr int s0 = 0 % spacedim;
1147 constexpr int s1 = 1 % spacedim;
1148 constexpr int s2 = 2 % spacedim;
1149
1150 // derivative w.r.t the radius
1151 derivatives[s0][s0] = intermediate[s0];
1152 derivatives[s1][s0] = intermediate[s1];
1153 derivatives[s2][s0] = intermediate[s2];
1154
1155 // derivatives w.r.t the angle
1156 derivatives[s0][s1] = -normal_direction[s0] * sine + dxn[s0] * cosine;
1157 derivatives[s1][s1] = -normal_direction[s1] * sine + dxn[s1] * cosine;
1158 derivatives[s2][s1] = -normal_direction[s2] * sine + dxn[s2] * cosine;
1159
1160 // derivatives w.r.t the direction of the axis
1161 derivatives[s0][s2] = direction[s0];
1162 derivatives[s1][s2] = direction[s1];
1163 derivatives[s2][s2] = direction[s2];
1164
1165 return derivatives;
1166}
1167
1168
1169
1170namespace
1171{
1172 template <int dim>
1174 check_and_normalize(const Tensor<1, dim> &t)
1175 {
1176 const double norm = t.norm();
1177 Assert(norm > 0.0, ExcMessage("The major axis must have a positive norm."));
1178 return t / norm;
1179 }
1180} // namespace
1181
1182
1183
1184// ============================================================
1185// EllipticalManifold
1186// ============================================================
1187template <int dim, int spacedim>
1189 const Point<spacedim> &center,
1190 const Tensor<1, spacedim> &major_axis_direction,
1191 const double eccentricity)
1192 : ChartManifold<dim, spacedim, spacedim>(
1193 EllipticalManifold<dim, spacedim>::get_periodicity())
1194 , direction(check_and_normalize(major_axis_direction))
1195 , center(center)
1196 , eccentricity(eccentricity)
1197 , cosh_u(1.0 / eccentricity)
1198 , sinh_u(std::sqrt(cosh_u * cosh_u - 1.0))
1199{
1200 // Throws an exception if dim!=2 || spacedim!=2.
1201 Assert(dim == 2 && spacedim == 2, ExcNotImplemented());
1202 // Throws an exception if eccentricity is not in range.
1203 Assert(std::signbit(cosh_u * cosh_u - 1.0) == false,
1204 ExcMessage(
1205 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1206}
1207
1208
1209
1210template <int dim, int spacedim>
1211std::unique_ptr<Manifold<dim, spacedim>>
1213{
1214 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1215}
1216
1217
1218
1219template <int dim, int spacedim>
1222{
1223 Tensor<1, spacedim> periodicity;
1224 // The second elliptical coordinate is periodic, while the first is not.
1225 // Enforce periodicity on the last variable.
1226 periodicity[spacedim - 1] = 2.0 * numbers::PI;
1227 return periodicity;
1228}
1229
1230
1231
1232template <int dim, int spacedim>
1239
1240
1241
1242template <>
1245{
1246 const double cs = std::cos(chart_point[1]);
1247 const double sn = std::sin(chart_point[1]);
1248 // Coordinates in the reference frame (i.e. major axis direction is
1249 // x-axis)
1250 const double x = chart_point[0] * cosh_u * cs;
1251 const double y = chart_point[0] * sinh_u * sn;
1252 // Rotate them according to the major axis direction
1253 const Point<2> p(direction[0] * x - direction[1] * y,
1254 direction[1] * x + direction[0] * y);
1255 return p + center;
1256}
1257
1258
1259
1260template <int dim, int spacedim>
1267
1268
1269
1270template <>
1273{
1274 // Moving space_point in the reference coordinate system.
1275 const double x0 = space_point[0] - center[0];
1276 const double y0 = space_point[1] - center[1];
1277 const double x = direction[0] * x0 + direction[1] * y0;
1278 const double y = -direction[1] * x0 + direction[0] * y0;
1279 const double pt0 =
1280 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1281 // If the radius is exactly zero, the point coincides with the origin.
1282 if (pt0 == 0.0)
1283 {
1284 return center;
1285 }
1286 double cos_eta = x / (pt0 * cosh_u);
1287 if (cos_eta < -1.0)
1288 {
1289 cos_eta = -1.0;
1290 }
1291 if (cos_eta > 1.0)
1292 {
1293 cos_eta = 1.0;
1294 }
1295 const double eta = std::acos(cos_eta);
1296 const double pt1 = (std::signbit(y) ? 2.0 * numbers::PI - eta : eta);
1297 return {pt0, pt1};
1298}
1299
1300
1301
1302template <int dim, int spacedim>
1310
1311
1312
1313template <>
1316 const Point<2> &chart_point) const
1317{
1318 const double cs = std::cos(chart_point[1]);
1319 const double sn = std::sin(chart_point[1]);
1320 Tensor<2, 2> dX;
1321 dX[0][0] = cosh_u * cs;
1322 dX[0][1] = -chart_point[0] * cosh_u * sn;
1323 dX[1][0] = sinh_u * sn;
1324 dX[1][1] = chart_point[0] * sinh_u * cs;
1325
1326 // rotate according to the major axis direction
1328 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1329
1330 return rot * dX;
1331}
1332
1333
1334
1335// ============================================================
1336// FunctionManifold
1337// ============================================================
1338template <int dim, int spacedim, int chartdim>
1340 const Function<chartdim> &push_forward_function,
1341 const Function<spacedim> &pull_back_function,
1342 const Tensor<1, chartdim> &periodicity,
1343 const double tolerance)
1344 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1345 , const_map()
1346 , push_forward_function(&push_forward_function)
1347 , pull_back_function(&pull_back_function)
1348 , tolerance(tolerance)
1349 , owns_pointers(false)
1350 , finite_difference_step(0)
1351{
1352 AssertDimension(push_forward_function.n_components, spacedim);
1353 AssertDimension(pull_back_function.n_components, chartdim);
1354}
1355
1356
1357
1358template <int dim, int spacedim, int chartdim>
1360 std::unique_ptr<Function<chartdim>> push_forward,
1361 std::unique_ptr<Function<spacedim>> pull_back,
1362 const Tensor<1, chartdim> &periodicity,
1363 const double tolerance)
1364 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1365 , const_map()
1366 , push_forward_function(push_forward.release())
1367 , pull_back_function(pull_back.release())
1368 , tolerance(tolerance)
1369 , owns_pointers(true)
1370 , finite_difference_step(0)
1371{
1372 AssertDimension(push_forward_function->n_components, spacedim);
1373 AssertDimension(pull_back_function->n_components, chartdim);
1374}
1375
1376
1377
1378template <int dim, int spacedim, int chartdim>
1380 const std::string push_forward_expression,
1381 const std::string pull_back_expression,
1382 const Tensor<1, chartdim> &periodicity,
1383 const typename FunctionParser<spacedim>::ConstMap const_map,
1384 const std::string chart_vars,
1385 const std::string space_vars,
1386 const double tolerance,
1387 const double h)
1388 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1389 , const_map(const_map)
1390 , tolerance(tolerance)
1391 , owns_pointers(true)
1392 , push_forward_expression(push_forward_expression)
1393 , pull_back_expression(pull_back_expression)
1394 , chart_vars(chart_vars)
1395 , space_vars(space_vars)
1396 , finite_difference_step(h)
1397{
1398 FunctionParser<chartdim> *pf = new FunctionParser<chartdim>(spacedim, 0.0, h);
1399 FunctionParser<spacedim> *pb = new FunctionParser<spacedim>(chartdim, 0.0, h);
1403 pull_back_function = pb;
1404}
1405
1406
1407
1408template <int dim, int spacedim, int chartdim>
1410{
1411 if (owns_pointers == true)
1412 {
1413 const Function<chartdim> *pf = push_forward_function;
1414 push_forward_function = nullptr;
1415 delete pf;
1416
1417 const Function<spacedim> *pb = pull_back_function;
1418 pull_back_function = nullptr;
1419 delete pb;
1420 }
1421}
1422
1423
1424
1425template <int dim, int spacedim, int chartdim>
1426std::unique_ptr<Manifold<dim, spacedim>>
1428{
1429 // This manifold can be constructed either by providing an expression for the
1430 // push forward and the pull back charts, or by providing two Function
1431 // objects. In the first case, the push_forward and pull_back functions are
1432 // created internally in FunctionManifold, and destroyed when this object is
1433 // deleted. In the second case, the function objects are destroyed if they
1434 // are passed as pointers upon construction.
1435 // We need to make sure that our cloned object is constructed in the
1436 // same way this class was constructed, and that its internal Function
1437 // pointers point either to the same Function objects used to construct this
1438 // function or that the newly generated manifold creates internally the
1439 // push_forward and pull_back functions using the same expressions that were
1440 // used to construct this class.
1441 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1442 {
1443 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1444 push_forward_expression,
1445 pull_back_expression,
1446 this->get_periodicity(),
1447 const_map,
1448 chart_vars,
1449 space_vars,
1450 tolerance,
1451 finite_difference_step);
1452 }
1453 else
1454 {
1455 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1456 *push_forward_function,
1457 *pull_back_function,
1458 this->get_periodicity(),
1459 tolerance);
1460 }
1461}
1462
1463
1464
1465template <int dim, int spacedim, int chartdim>
1468 const Point<chartdim> &chart_point) const
1469{
1470 Vector<double> pf(spacedim);
1471 Point<spacedim> result;
1472 push_forward_function->vector_value(chart_point, pf);
1473 for (unsigned int i = 0; i < spacedim; ++i)
1474 result[i] = pf[i];
1475
1476#ifdef DEBUG
1477 Vector<double> pb(chartdim);
1478 pull_back_function->vector_value(result, pb);
1479 for (unsigned int i = 0; i < chartdim; ++i)
1480 Assert(
1481 (chart_point.norm() > tolerance &&
1482 (std::abs(pb[i] - chart_point[i]) < tolerance * chart_point.norm())) ||
1483 (std::abs(pb[i] - chart_point[i]) < tolerance),
1484 ExcMessage(
1485 "The push forward is not the inverse of the pull back! Bailing out."));
1486#endif
1487
1488 return result;
1489}
1490
1491
1492
1493template <int dim, int spacedim, int chartdim>
1496 const Point<chartdim> &chart_point) const
1497{
1499 for (unsigned int i = 0; i < spacedim; ++i)
1500 {
1501 const auto gradient = push_forward_function->gradient(chart_point, i);
1502 for (unsigned int j = 0; j < chartdim; ++j)
1503 DF[i][j] = gradient[j];
1504 }
1505 return DF;
1506}
1507
1508
1509
1510template <int dim, int spacedim, int chartdim>
1513 const Point<spacedim> &space_point) const
1514{
1515 Vector<double> pb(chartdim);
1516 Point<chartdim> result;
1517 pull_back_function->vector_value(space_point, pb);
1518 for (unsigned int i = 0; i < chartdim; ++i)
1519 result[i] = pb[i];
1520 return result;
1521}
1522
1523
1524
1525// ============================================================
1526// TorusManifold
1527// ============================================================
1528template <int dim>
1531{
1532 double x = p[0];
1533 double z = p[1];
1534 double y = p[2];
1535 double phi = std::atan2(y, x);
1536 double theta = std::atan2(z, std::sqrt(x * x + y * y) - centerline_radius);
1537 double w =
1538 std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * centerline_radius) +
1539 Utilities::fixed_power<2>(x - std::cos(phi) * centerline_radius) +
1540 z * z) /
1541 inner_radius;
1542 return {phi, theta, w};
1543}
1544
1545
1546
1547template <int dim>
1550{
1551 double phi = chart_point[0];
1552 double theta = chart_point[1];
1553 double w = chart_point[2];
1554
1555 return {std::cos(phi) * centerline_radius +
1556 inner_radius * w * std::cos(theta) * std::cos(phi),
1557 inner_radius * w * std::sin(theta),
1558 std::sin(phi) * centerline_radius +
1559 inner_radius * w * std::cos(theta) * std::sin(phi)};
1560}
1561
1562
1563
1564template <int dim>
1565TorusManifold<dim>::TorusManifold(const double centerline_radius,
1566 const double inner_radius)
1567 : ChartManifold<dim, 3, 3>(Point<3>(2 * numbers::PI, 2 * numbers::PI, 0.0))
1568 , centerline_radius(centerline_radius)
1569 , inner_radius(inner_radius)
1570{
1572 ExcMessage("The centerline radius must be greater than the "
1573 "inner radius."));
1574 Assert(inner_radius > 0.0, ExcMessage("The inner radius must be positive."));
1575}
1576
1577
1578
1579template <int dim>
1580std::unique_ptr<Manifold<dim, 3>>
1582{
1583 return std::make_unique<TorusManifold<dim>>(centerline_radius, inner_radius);
1584}
1585
1586
1587
1588template <int dim>
1591{
1593
1594 double phi = chart_point[0];
1595 double theta = chart_point[1];
1596 double w = chart_point[2];
1597
1598 DX[0][0] = -std::sin(phi) * centerline_radius -
1599 inner_radius * w * std::cos(theta) * std::sin(phi);
1600 DX[0][1] = -inner_radius * w * std::sin(theta) * std::cos(phi);
1601 DX[0][2] = inner_radius * std::cos(theta) * std::cos(phi);
1602
1603 DX[1][0] = 0;
1604 DX[1][1] = inner_radius * w * std::cos(theta);
1605 DX[1][2] = inner_radius * std::sin(theta);
1606
1607 DX[2][0] = std::cos(phi) * centerline_radius +
1608 inner_radius * w * std::cos(theta) * std::cos(phi);
1609 DX[2][1] = -inner_radius * w * std::sin(theta) * std::sin(phi);
1610 DX[2][2] = inner_radius * std::cos(theta) * std::sin(phi);
1611
1612 return DX;
1613}
1614
1615
1616
1617// ============================================================
1618// TransfiniteInterpolationManifold
1619// ============================================================
1620template <int dim, int spacedim>
1623 : triangulation(nullptr)
1624 , level_coarse(-1)
1625{
1626 AssertThrow(dim > 1, ExcNotImplemented());
1627}
1628
1629
1630
1631template <int dim, int spacedim>
1633 spacedim>::~TransfiniteInterpolationManifold()
1634{
1635 if (clear_signal.connected())
1636 clear_signal.disconnect();
1637}
1638
1639
1640
1641template <int dim, int spacedim>
1642std::unique_ptr<Manifold<dim, spacedim>>
1644{
1646 if (triangulation)
1647 ptr->initialize(*triangulation);
1648 return std::unique_ptr<Manifold<dim, spacedim>>(ptr);
1649}
1650
1651
1652
1653template <int dim, int spacedim>
1654void
1657{
1658 this->triangulation = &triangulation;
1659 // In case the triangulation is cleared, remove the pointers by a signal:
1660 clear_signal.disconnect();
1661 clear_signal = triangulation.signals.clear.connect([&]() -> void {
1662 this->triangulation = nullptr;
1663 this->level_coarse = -1;
1664 });
1665 level_coarse = triangulation.last()->level();
1666 coarse_cell_is_flat.resize(triangulation.n_cells(level_coarse), false);
1667 quadratic_approximation.clear();
1668
1669 // In case of dim == spacedim we perform a quadratic approximation in
1670 // InverseQuadraticApproximation(), thus initialize the unit_points
1671 // vector with one subdivision to get 3^dim unit_points.
1672 //
1673 // In the co-dimension one case (meaning dim < spacedim) we have to fall
1674 // back to a simple GridTools::affine_cell_approximation<dim>() which
1675 // requires 2^dim points, instead. Thus, initialize the QIterated
1676 // quadrature with no subdivisions.
1677 const std::vector<Point<dim>> unit_points =
1678 QIterated<dim>(QTrapezoid<1>(), (dim == spacedim ? 2 : 1)).get_points();
1679 std::vector<Point<spacedim>> real_points(unit_points.size());
1680
1681 for (const auto &cell : triangulation.active_cell_iterators())
1682 {
1683 bool cell_is_flat = true;
1684 for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
1685 if (cell->line(l)->manifold_id() != cell->manifold_id() &&
1686 cell->line(l)->manifold_id() != numbers::flat_manifold_id)
1687 cell_is_flat = false;
1688 if (dim > 2)
1689 for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
1690 if (cell->quad(q)->manifold_id() != cell->manifold_id() &&
1691 cell->quad(q)->manifold_id() != numbers::flat_manifold_id)
1692 cell_is_flat = false;
1693 AssertIndexRange(static_cast<unsigned int>(cell->index()),
1694 coarse_cell_is_flat.size());
1695 coarse_cell_is_flat[cell->index()] = cell_is_flat;
1696
1697 // build quadratic interpolation
1698 for (unsigned int i = 0; i < unit_points.size(); ++i)
1699 real_points[i] = push_forward(cell, unit_points[i]);
1700 quadratic_approximation.emplace_back(real_points, unit_points);
1701 }
1702}
1703
1704
1705
1706namespace
1707{
1708 // version for 1d
1709 template <typename AccessorType>
1711 compute_transfinite_interpolation(const AccessorType &cell,
1712 const Point<1> &chart_point,
1713 const bool /*cell_is_flat*/)
1714 {
1715 return cell.vertex(0) * (1. - chart_point[0]) +
1716 cell.vertex(1) * chart_point[0];
1717 }
1718
1719 // version for 2d
1720 template <typename AccessorType>
1722 compute_transfinite_interpolation(const AccessorType &cell,
1723 const Point<2> &chart_point,
1724 const bool cell_is_flat)
1725 {
1726 const unsigned int dim = AccessorType::dimension;
1727 const unsigned int spacedim = AccessorType::space_dimension;
1728 const types::manifold_id my_manifold_id = cell.manifold_id();
1729 const Triangulation<dim, spacedim> &tria = cell.get_triangulation();
1730
1731 // formula see wikipedia
1732 // https://en.wikipedia.org/wiki/Transfinite_interpolation
1733 // S(u,v) = (1-v)c_1(u)+v c_3(u) + (1-u)c_2(v) + u c_4(v) -
1734 // [(1-u)(1-v)P_0 + u(1-v) P_1 + (1-u)v P_2 + uv P_3]
1735 const std::array<Point<spacedim>, 4> vertices{
1736 {cell.vertex(0), cell.vertex(1), cell.vertex(2), cell.vertex(3)}};
1737
1738 // this evaluates all bilinear shape functions because we need them
1739 // repeatedly. we will update this values in the complicated case with
1740 // curved lines below
1741 std::array<double, 4> weights_vertices{
1742 {(1. - chart_point[0]) * (1. - chart_point[1]),
1743 chart_point[0] * (1. - chart_point[1]),
1744 (1. - chart_point[0]) * chart_point[1],
1745 chart_point[0] * chart_point[1]}};
1746
1747 Point<spacedim> new_point;
1748 if (cell_is_flat)
1749 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1750 new_point += weights_vertices[v] * vertices[v];
1751 else
1752 {
1753 // The second line in the formula tells us to subtract the
1754 // contribution of the vertices. If a line employs the same manifold
1755 // as the cell, we can merge the weights of the line with the weights
1756 // of the vertex with a negative sign while going through the faces
1757 // (this is a bit artificial in 2d but it becomes clear in 3d where we
1758 // avoid looking at the faces' orientation and other complications).
1759
1760 // add the contribution from the lines around the cell (first line in
1761 // formula)
1762 std::array<double, GeometryInfo<2>::vertices_per_face> weights;
1763 std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_face> points;
1764 // note that the views are immutable, but the arrays are not
1765 const auto weights_view =
1766 make_array_view(weights.begin(), weights.end());
1767 const auto points_view = make_array_view(points.begin(), points.end());
1768
1769 for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell;
1770 ++line)
1771 {
1772 const double my_weight =
1773 (line % 2) ? chart_point[line / 2] : 1 - chart_point[line / 2];
1774 const double line_point = chart_point[1 - line / 2];
1775
1776 // Same manifold or invalid id which will go back to the same
1777 // class -> contribution should be added for the final point,
1778 // which means that we subtract the current weight from the
1779 // negative weight applied to the vertex
1780 const types::manifold_id line_manifold_id =
1781 cell.line(line)->manifold_id();
1782 if (line_manifold_id == my_manifold_id ||
1783 line_manifold_id == numbers::flat_manifold_id)
1784 {
1785 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,
1786 0)] -=
1787 my_weight * (1. - line_point);
1788 weights_vertices[GeometryInfo<2>::line_to_cell_vertices(line,
1789 1)] -=
1790 my_weight * line_point;
1791 }
1792 else
1793 {
1794 points[0] =
1796 points[1] =
1798 weights[0] = 1. - line_point;
1799 weights[1] = line_point;
1800 new_point +=
1801 my_weight * tria.get_manifold(line_manifold_id)
1802 .get_new_point(points_view, weights_view);
1803 }
1804 }
1805
1806 // subtract contribution from the vertices (second line in formula)
1807 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1808 new_point -= weights_vertices[v] * vertices[v];
1809 }
1810
1811 return new_point;
1812 }
1813
1814 // this is replicated from GeometryInfo::face_to_cell_vertices since we need
1815 // it very often in compute_transfinite_interpolation and the function is
1816 // performance critical
1817 static constexpr unsigned int face_to_cell_vertices_3d[6][4] = {{0, 2, 4, 6},
1818 {1, 3, 5, 7},
1819 {0, 4, 1, 5},
1820 {2, 6, 3, 7},
1821 {0, 1, 2, 3},
1822 {4, 5, 6, 7}};
1823
1824 // this is replicated from GeometryInfo::face_to_cell_lines since we need it
1825 // very often in compute_transfinite_interpolation and the function is
1826 // performance critical
1827 static constexpr unsigned int face_to_cell_lines_3d[6][4] = {{8, 10, 0, 4},
1828 {9, 11, 1, 5},
1829 {2, 6, 8, 9},
1830 {3, 7, 10, 11},
1831 {0, 1, 2, 3},
1832 {4, 5, 6, 7}};
1833
1834 // version for 3d
1835 template <typename AccessorType>
1837 compute_transfinite_interpolation(const AccessorType &cell,
1838 const Point<3> &chart_point,
1839 const bool cell_is_flat)
1840 {
1841 const unsigned int dim = AccessorType::dimension;
1842 const unsigned int spacedim = AccessorType::space_dimension;
1843 const types::manifold_id my_manifold_id = cell.manifold_id();
1844 const Triangulation<dim, spacedim> &tria = cell.get_triangulation();
1845
1846 // Same approach as in 2d, but adding the faces, subtracting the edges, and
1847 // adding the vertices
1848 const std::array<Point<spacedim>, 8> vertices{{cell.vertex(0),
1849 cell.vertex(1),
1850 cell.vertex(2),
1851 cell.vertex(3),
1852 cell.vertex(4),
1853 cell.vertex(5),
1854 cell.vertex(6),
1855 cell.vertex(7)}};
1856
1857 // store the components of the linear shape functions because we need them
1858 // repeatedly. we allow for 10 such shape functions to wrap around the
1859 // first four once again for easier face access.
1860 double linear_shapes[10];
1861 for (unsigned int d = 0; d < 3; ++d)
1862 {
1863 linear_shapes[2 * d] = 1. - chart_point[d];
1864 linear_shapes[2 * d + 1] = chart_point[d];
1865 }
1866
1867 // wrap linear shape functions around for access in face loop
1868 for (unsigned int d = 6; d < 10; ++d)
1869 linear_shapes[d] = linear_shapes[d - 6];
1870
1871 std::array<double, 8> weights_vertices;
1872 for (unsigned int i2 = 0, v = 0; i2 < 2; ++i2)
1873 for (unsigned int i1 = 0; i1 < 2; ++i1)
1874 for (unsigned int i0 = 0; i0 < 2; ++i0, ++v)
1875 weights_vertices[v] =
1876 (linear_shapes[4 + i2] * linear_shapes[2 + i1]) * linear_shapes[i0];
1877
1878 Point<spacedim> new_point;
1879 if (cell_is_flat)
1880 for (unsigned int v = 0; v < 8; ++v)
1881 new_point += weights_vertices[v] * vertices[v];
1882 else
1883 {
1884 // identify the weights for the lines to be accumulated (vertex
1885 // weights are set outside and coincide with the flat manifold case)
1886
1887 std::array<double, GeometryInfo<3>::lines_per_cell> weights_lines;
1888 std::fill(weights_lines.begin(), weights_lines.end(), 0.0);
1889
1890 // start with the contributions of the faces
1891 std::array<double, GeometryInfo<2>::vertices_per_cell> weights;
1892 std::array<Point<spacedim>, GeometryInfo<2>::vertices_per_cell> points;
1893 // note that the views are immutable, but the arrays are not
1894 const auto weights_view =
1895 make_array_view(weights.begin(), weights.end());
1896 const auto points_view = make_array_view(points.begin(), points.end());
1897
1898 for (const unsigned int face : GeometryInfo<3>::face_indices())
1899 {
1900 const double my_weight = linear_shapes[face];
1901 const unsigned int face_even = face - face % 2;
1902
1903 if (std::abs(my_weight) < 1e-13)
1904 continue;
1905
1906 // same manifold or invalid id which will go back to the same class
1907 // -> face will interpolate from the surrounding lines and vertices
1908 const types::manifold_id face_manifold_id =
1909 cell.face(face)->manifold_id();
1910 if (face_manifold_id == my_manifold_id ||
1911 face_manifold_id == numbers::flat_manifold_id)
1912 {
1913 for (unsigned int line = 0;
1914 line < GeometryInfo<2>::lines_per_cell;
1915 ++line)
1916 {
1917 const double line_weight =
1918 linear_shapes[face_even + 2 + line];
1919 weights_lines[face_to_cell_lines_3d[face][line]] +=
1920 my_weight * line_weight;
1921 }
1922 // as to the indices inside linear_shapes: we use the index
1923 // wrapped around at 2*d, ensuring the correct orientation of
1924 // the face's coordinate system with respect to the
1925 // lexicographic indices
1926 weights_vertices[face_to_cell_vertices_3d[face][0]] -=
1927 linear_shapes[face_even + 2] *
1928 (linear_shapes[face_even + 4] * my_weight);
1929 weights_vertices[face_to_cell_vertices_3d[face][1]] -=
1930 linear_shapes[face_even + 3] *
1931 (linear_shapes[face_even + 4] * my_weight);
1932 weights_vertices[face_to_cell_vertices_3d[face][2]] -=
1933 linear_shapes[face_even + 2] *
1934 (linear_shapes[face_even + 5] * my_weight);
1935 weights_vertices[face_to_cell_vertices_3d[face][3]] -=
1936 linear_shapes[face_even + 3] *
1937 (linear_shapes[face_even + 5] * my_weight);
1938 }
1939 else
1940 {
1941 for (const unsigned int v : GeometryInfo<2>::vertex_indices())
1942 points[v] = vertices[face_to_cell_vertices_3d[face][v]];
1943 weights[0] =
1944 linear_shapes[face_even + 2] * linear_shapes[face_even + 4];
1945 weights[1] =
1946 linear_shapes[face_even + 3] * linear_shapes[face_even + 4];
1947 weights[2] =
1948 linear_shapes[face_even + 2] * linear_shapes[face_even + 5];
1949 weights[3] =
1950 linear_shapes[face_even + 3] * linear_shapes[face_even + 5];
1951 new_point +=
1952 my_weight * tria.get_manifold(face_manifold_id)
1953 .get_new_point(points_view, weights_view);
1954 }
1955 }
1956
1957 // next subtract the contributions of the lines
1958 const auto weights_view_line =
1959 make_array_view(weights.begin(), weights.begin() + 2);
1960 const auto points_view_line =
1961 make_array_view(points.begin(), points.begin() + 2);
1962 for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell;
1963 ++line)
1964 {
1965 const double line_point =
1966 (line < 8 ? chart_point[1 - (line % 4) / 2] : chart_point[2]);
1967 double my_weight = 0.;
1968 if (line < 8)
1969 my_weight = linear_shapes[line % 4] * linear_shapes[4 + line / 4];
1970 else
1971 {
1972 const unsigned int subline = line - 8;
1973 my_weight =
1974 linear_shapes[subline % 2] * linear_shapes[2 + subline / 2];
1975 }
1976 my_weight -= weights_lines[line];
1977
1978 if (std::abs(my_weight) < 1e-13)
1979 continue;
1980
1981 const types::manifold_id line_manifold_id =
1982 cell.line(line)->manifold_id();
1983 if (line_manifold_id == my_manifold_id ||
1984 line_manifold_id == numbers::flat_manifold_id)
1985 {
1986 weights_vertices[GeometryInfo<3>::line_to_cell_vertices(line,
1987 0)] -=
1988 my_weight * (1. - line_point);
1989 weights_vertices[GeometryInfo<3>::line_to_cell_vertices(line,
1990 1)] -=
1991 my_weight * (line_point);
1992 }
1993 else
1994 {
1995 points[0] =
1997 points[1] =
1999 weights[0] = 1. - line_point;
2000 weights[1] = line_point;
2001 new_point -= my_weight * tria.get_manifold(line_manifold_id)
2002 .get_new_point(points_view_line,
2003 weights_view_line);
2004 }
2005 }
2006
2007 // finally add the contribution of the
2008 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2009 new_point += weights_vertices[v] * vertices[v];
2010 }
2011 return new_point;
2012 }
2013} // namespace
2014
2015
2016
2017template <int dim, int spacedim>
2021 const Point<dim> &chart_point) const
2022{
2023 AssertDimension(cell->level(), level_coarse);
2024
2025 // check that the point is in the unit cell which is the current chart
2026 // Tolerance 5e-4 chosen that the method also works with manifolds
2027 // that have some discretization error like SphericalManifold
2029 ExcMessage("chart_point is not in unit interval"));
2030
2031 return compute_transfinite_interpolation(*cell,
2032 chart_point,
2033 coarse_cell_is_flat[cell->index()]);
2034}
2035
2036
2037
2038template <int dim, int spacedim>
2042 const Point<dim> &chart_point,
2043 const Point<spacedim> &pushed_forward_chart_point) const
2044{
2045 // compute the derivative with the help of finite differences
2047 for (unsigned int d = 0; d < dim; ++d)
2048 {
2049 Point<dim> modified = chart_point;
2050 const double step = chart_point[d] > 0.5 ? -1e-8 : 1e-8;
2051
2052 // avoid checking outside of the unit interval
2053 modified[d] += step;
2054 Tensor<1, spacedim> difference =
2055 compute_transfinite_interpolation(*cell,
2056 modified,
2057 coarse_cell_is_flat[cell->index()]) -
2058 pushed_forward_chart_point;
2059 for (unsigned int e = 0; e < spacedim; ++e)
2060 grad[e][d] = difference[e] / step;
2061 }
2062 return grad;
2063}
2064
2065
2066
2067template <int dim, int spacedim>
2071 const Point<spacedim> &point,
2072 const Point<dim> &initial_guess) const
2073{
2074 Point<dim> outside;
2075 for (unsigned int d = 0; d < dim; ++d)
2077
2078 // project the user-given input to unit cell
2079 Point<dim> chart_point = cell->reference_cell().closest_point(initial_guess);
2080
2081 // run quasi-Newton iteration with a combination of finite differences for
2082 // the exact Jacobian and "Broyden's good method". As opposed to the various
2083 // mapping implementations, this class does not throw exception upon failure
2084 // as those are relatively expensive and failure occurs quite regularly in
2085 // the implementation of the compute_chart_points method.
2086 Tensor<1, spacedim> residual =
2087 point -
2088 compute_transfinite_interpolation(*cell,
2089 chart_point,
2090 coarse_cell_is_flat[cell->index()]);
2091 const double tolerance = 1e-21 * Utilities::fixed_power<2>(cell->diameter());
2092 double residual_norm_square = residual.norm_square();
2094 bool must_recompute_jacobian = true;
2095 for (unsigned int i = 0; i < 100; ++i)
2096 {
2097 if (residual_norm_square < tolerance)
2098 {
2099 // do a final update of the point with the last available Jacobian
2100 // information. The residual is close to zero due to the check
2101 // above, but me might improve some of the last digits by a final
2102 // Newton-like step with step length 1
2103 Tensor<1, dim> update;
2104 for (unsigned int d = 0; d < spacedim; ++d)
2105 for (unsigned int e = 0; e < dim; ++e)
2106 update[e] += inv_grad[d][e] * residual[d];
2107 return chart_point + update;
2108 }
2109
2110 // every 9 iterations, including the first time around, we create an
2111 // approximation of the Jacobian with finite differences. Broyden's
2112 // method usually does not need more than 5-8 iterations, but sometimes
2113 // we might have had a bad initial guess and then we can accelerate
2114 // convergence considerably with getting the actual Jacobian rather than
2115 // using secant-like methods (one gradient calculation in 3d costs as
2116 // much as 3 more iterations). this usually happens close to convergence
2117 // and one more step with the finite-differenced Jacobian leads to
2118 // convergence
2119 if (must_recompute_jacobian || i % 9 == 0)
2120 {
2121 // if the determinant is zero or negative, the mapping is either not
2122 // invertible or already has inverted and we are outside the valid
2123 // chart region. Note that the Jacobian here represents the
2124 // derivative of the forward map and should have a positive
2125 // determinant since we use properly oriented meshes.
2127 push_forward_gradient(cell,
2128 chart_point,
2129 Point<spacedim>(point - residual));
2130 if (grad.determinant() <= 0.0)
2131 return outside;
2132 inv_grad = grad.covariant_form();
2133 must_recompute_jacobian = false;
2134 }
2135 Tensor<1, dim> update;
2136 for (unsigned int d = 0; d < spacedim; ++d)
2137 for (unsigned int e = 0; e < dim; ++e)
2138 update[e] += inv_grad[d][e] * residual[d];
2139
2140 // Line search, accept step if the residual has decreased
2141 double alpha = 1.;
2142
2143 // check if point is inside 1.2 times the unit cell to avoid
2144 // hitting points very far away from valid ones in the manifolds
2145 while (
2146 !GeometryInfo<dim>::is_inside_unit_cell(chart_point + alpha * update,
2147 0.2) &&
2148 alpha > 1e-7)
2149 alpha *= 0.5;
2150
2151 const Tensor<1, spacedim> old_residual = residual;
2152 while (alpha > 1e-4)
2153 {
2154 Point<dim> guess = chart_point + alpha * update;
2155 const Tensor<1, spacedim> residual_guess =
2156 point - compute_transfinite_interpolation(
2157 *cell, guess, coarse_cell_is_flat[cell->index()]);
2158 const double residual_norm_new = residual_guess.norm_square();
2159 if (residual_norm_new < residual_norm_square)
2160 {
2161 residual = residual_guess;
2162 residual_norm_square = residual_norm_new;
2163 chart_point += alpha * update;
2164 break;
2165 }
2166 else
2167 alpha *= 0.5;
2168 }
2169 // If alpha got very small, it is likely due to a bad Jacobian
2170 // approximation with Broyden's method (relatively far away from the
2171 // zero), which can be corrected by the outer loop when a Newton update
2172 // is recomputed. The second case is when the Jacobian is actually bad
2173 // and we should fail as early as possible. Since we cannot really
2174 // distinguish the two, we must continue here in any case.
2175 if (alpha <= 1e-4)
2176 must_recompute_jacobian = true;
2177
2178 // update the inverse Jacobian with "Broyden's good method" and
2179 // Sherman-Morrison formula for the update of the inverse, see
2180 // https://en.wikipedia.org/wiki/Broyden%27s_method
2181 // J^{-1}_n = J^{-1}_{n-1} + (delta x_n - J^{-1}_{n-1} delta f_n) /
2182 // (delta x_n^T J_{-1}_{n-1} delta f_n) delta x_n^T J^{-1}_{n-1}
2183
2184 // switch sign in residual as compared to the formula above because we
2185 // use a negative definition of the residual with respect to the
2186 // Jacobian
2187 const Tensor<1, spacedim> delta_f = old_residual - residual;
2188
2189 Tensor<1, dim> Jinv_deltaf;
2190 for (unsigned int d = 0; d < spacedim; ++d)
2191 for (unsigned int e = 0; e < dim; ++e)
2192 Jinv_deltaf[e] += inv_grad[d][e] * delta_f[d];
2193
2194 const Tensor<1, dim> delta_x = alpha * update;
2195
2196 // prevent division by zero. This number should be scale-invariant
2197 // because Jinv_deltaf carries no units and x is in reference
2198 // coordinates.
2199 if (std::abs(delta_x * Jinv_deltaf) > 1e-12 && !must_recompute_jacobian)
2200 {
2201 const Tensor<1, dim> factor =
2202 (delta_x - Jinv_deltaf) / (delta_x * Jinv_deltaf);
2203 Tensor<1, spacedim> jac_update;
2204 for (unsigned int d = 0; d < spacedim; ++d)
2205 for (unsigned int e = 0; e < dim; ++e)
2206 jac_update[d] += delta_x[e] * inv_grad[d][e];
2207 for (unsigned int d = 0; d < spacedim; ++d)
2208 for (unsigned int e = 0; e < dim; ++e)
2209 inv_grad[d][e] += factor[e] * jac_update[d];
2210 }
2211 }
2212 return outside;
2213}
2214
2215
2216
2217template <int dim, int spacedim>
2218std::array<unsigned int, 20>
2221 const ArrayView<const Point<spacedim>> &points) const
2222{
2223 // The methods to identify cells around points in GridTools are all written
2224 // for the active cells, but we are here looking at some cells at the coarse
2225 // level.
2226 Assert(triangulation != nullptr, ExcNotInitialized());
2227 Assert(triangulation->begin_active()->level() >= level_coarse,
2228 ExcMessage("The manifold was initialized with level " +
2229 std::to_string(level_coarse) + " but there are now" +
2230 "active cells on a lower level. Coarsening the mesh is " +
2231 "currently not supported"));
2232
2233 // This computes the distance of the surrounding points transformed to the
2234 // unit cell from the unit cell.
2237 level_coarse),
2238 endc =
2240 level_coarse);
2241 boost::container::small_vector<std::pair<double, unsigned int>, 200>
2242 distances_and_cells;
2243 for (; cell != endc; ++cell)
2244 {
2245 // only consider cells where the current manifold is attached
2246 if (&cell->get_manifold() != this)
2247 continue;
2248
2249 std::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
2250 vertices;
2251 for (const unsigned int vertex_n : GeometryInfo<dim>::vertex_indices())
2252 {
2253 vertices[vertex_n] = cell->vertex(vertex_n);
2254 }
2255
2256 // cheap check: if any of the points is not inside a circle around the
2257 // center of the loop, we can skip the expensive part below (this assumes
2258 // that the manifold does not deform the grid too much)
2260 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2261 center += vertices[v];
2263 double radius_square = 0.;
2264 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2265 radius_square =
2266 std::max(radius_square, (center - vertices[v]).norm_square());
2267 bool inside_circle = true;
2268 for (unsigned int i = 0; i < points.size(); ++i)
2269 if ((center - points[i]).norm_square() > radius_square * 1.5)
2270 {
2271 inside_circle = false;
2272 break;
2273 }
2274 if (inside_circle == false)
2275 continue;
2276
2277 // slightly more expensive search
2278 double current_distance = 0;
2279 for (unsigned int i = 0; i < points.size(); ++i)
2280 {
2281 Point<dim> point =
2282 quadratic_approximation[cell->index()].compute(points[i]);
2283 current_distance += GeometryInfo<dim>::distance_to_unit_cell(point);
2284 }
2285 distances_and_cells.push_back(
2286 std::make_pair(current_distance, cell->index()));
2287 }
2288 // no coarse cell could be found -> transformation failed
2289 AssertThrow(distances_and_cells.size() > 0,
2291 std::sort(distances_and_cells.begin(), distances_and_cells.end());
2292 std::array<unsigned int, 20> cells;
2294 for (unsigned int i = 0; i < distances_and_cells.size() && i < cells.size();
2295 ++i)
2296 cells[i] = distances_and_cells[i].second;
2297
2298 return cells;
2299}
2300
2301
2302
2303template <int dim, int spacedim>
2306 const ArrayView<const Point<spacedim>> &surrounding_points,
2307 ArrayView<Point<dim>> chart_points) const
2308{
2309 Assert(surrounding_points.size() == chart_points.size(),
2310 ExcMessage("The chart points array view must be as large as the "
2311 "surrounding points array view."));
2312
2313 std::array<unsigned int, 20> nearby_cells =
2314 get_possible_cells_around_points(surrounding_points);
2315
2316 // This function is nearly always called to place new points on a cell or
2317 // cell face. In this case, the general structure of the surrounding points
2318 // is known (i.e., if there are eight surrounding points, then they will
2319 // almost surely be either eight points around a quadrilateral or the eight
2320 // vertices of a cube). Hence, making this assumption, we use two
2321 // optimizations (one for structdim == 2 and one for structdim == 3) that
2322 // guess the locations of some of the chart points more efficiently than the
2323 // affine map approximation. The affine map approximation is used whenever
2324 // we don't have a cheaper guess available.
2325
2326 // Function that can guess the location of a chart point by assuming that
2327 // the eight surrounding points are points on a two-dimensional object
2328 // (either a cell in 2d or the face of a hexahedron in 3d), arranged like
2329 //
2330 // 2 - 7 - 3
2331 // | |
2332 // 4 5
2333 // | |
2334 // 0 - 6 - 1
2335 //
2336 // This function assumes that the first three chart points have been
2337 // computed since there is no effective way to guess them.
2338 auto guess_chart_point_structdim_2 = [&](const unsigned int i) -> Point<dim> {
2339 Assert(surrounding_points.size() == 8 && 2 < i && i < 8,
2340 ExcMessage("This function assumes that there are eight surrounding "
2341 "points around a two-dimensional object. It also assumes "
2342 "that the first three chart points have already been "
2343 "computed."));
2344 switch (i)
2345 {
2346 case 0:
2347 case 1:
2348 case 2:
2350 break;
2351 case 3:
2352 return chart_points[1] + (chart_points[2] - chart_points[0]);
2353 case 4:
2354 return 0.5 * (chart_points[0] + chart_points[2]);
2355 case 5:
2356 return 0.5 * (chart_points[1] + chart_points[3]);
2357 case 6:
2358 return 0.5 * (chart_points[0] + chart_points[1]);
2359 case 7:
2360 return 0.5 * (chart_points[2] + chart_points[3]);
2361 default:
2363 }
2364
2365 return {};
2366 };
2367
2368 // Function that can guess the location of a chart point by assuming that
2369 // the eight surrounding points form the vertices of a hexahedron, arranged
2370 // like
2371 //
2372 // 6-------7
2373 // /| /|
2374 // / / |
2375 // / | / |
2376 // 4-------5 |
2377 // | 2- -|- -3
2378 // | / | /
2379 // | | /
2380 // |/ |/
2381 // 0-------1
2382 //
2383 // (where vertex 2 is the back left vertex) we can estimate where chart
2384 // points 5 - 7 are by computing the height (in chart coordinates) as c4 -
2385 // c0 and then adding that onto the appropriate bottom vertex.
2386 //
2387 // This function assumes that the first five chart points have been computed
2388 // since there is no effective way to guess them.
2389 auto guess_chart_point_structdim_3 = [&](const unsigned int i) -> Point<dim> {
2390 Assert(surrounding_points.size() == 8 && 4 < i && i < 8,
2391 ExcMessage("This function assumes that there are eight surrounding "
2392 "points around a three-dimensional object. It also "
2393 "assumes that the first five chart points have already "
2394 "been computed."));
2395 return chart_points[i - 4] + (chart_points[4] - chart_points[0]);
2396 };
2397
2398 // Check if we can use the two chart point shortcuts above before we start:
2399 bool use_structdim_2_guesses = false;
2400 bool use_structdim_3_guesses = false;
2401 // note that in the structdim 2 case: 0 - 6 and 2 - 7 should be roughly
2402 // parallel, while in the structdim 3 case, 0 - 6 and 2 - 7 should be roughly
2403 // orthogonal. Use the angle between these two vectors to figure out if we
2404 // should turn on either structdim optimization.
2405 if (surrounding_points.size() == 8)
2406 {
2407 const Tensor<1, spacedim> v06 =
2408 surrounding_points[6] - surrounding_points[0];
2409 const Tensor<1, spacedim> v27 =
2410 surrounding_points[7] - surrounding_points[2];
2411
2412 // note that we can save a call to sqrt() by rearranging
2413 const double cosine = scalar_product(v06, v27) /
2414 std::sqrt(v06.norm_square() * v27.norm_square());
2415 if (0.707 < cosine)
2416 // the angle is less than pi/4, so these vectors are roughly parallel:
2417 // enable the structdim 2 optimization
2418 use_structdim_2_guesses = true;
2419 else if (spacedim == 3)
2420 // otherwise these vectors are roughly orthogonal: enable the
2421 // structdim 3 optimization if we are in 3d
2422 use_structdim_3_guesses = true;
2423 }
2424 // we should enable at most one of the optimizations
2425 Assert((!use_structdim_2_guesses && !use_structdim_3_guesses) ||
2426 (use_structdim_2_guesses ^ use_structdim_3_guesses),
2428
2429
2430
2431 auto compute_chart_point =
2432 [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2433 const unsigned int point_index) {
2434 Point<dim> guess;
2435 // an optimization: keep track of whether or not we used the quadratic
2436 // approximation so that we don't call pull_back with the same
2437 // initial guess twice (i.e., if pull_back fails the first time,
2438 // don't try again with the same function arguments).
2439 bool used_quadratic_approximation = false;
2440 // if we have already computed three points, we can guess the fourth
2441 // to be the missing corner point of a rectangle
2442 if (point_index == 3 && surrounding_points.size() >= 8)
2443 guess = chart_points[1] + (chart_points[2] - chart_points[0]);
2444 else if (use_structdim_2_guesses && 3 < point_index)
2445 guess = guess_chart_point_structdim_2(point_index);
2446 else if (use_structdim_3_guesses && 4 < point_index)
2447 guess = guess_chart_point_structdim_3(point_index);
2448 else if (dim == 3 && point_index > 7 && surrounding_points.size() == 26)
2449 {
2450 if (point_index < 20)
2451 guess =
2452 0.5 * (chart_points[GeometryInfo<dim>::line_to_cell_vertices(
2453 point_index - 8, 0)] +
2455 point_index - 8, 1)]);
2456 else
2457 guess =
2458 0.25 * (chart_points[GeometryInfo<dim>::face_to_cell_vertices(
2459 point_index - 20, 0)] +
2461 point_index - 20, 1)] +
2463 point_index - 20, 2)] +
2465 point_index - 20, 3)]);
2466 }
2467 else
2468 {
2469 guess = quadratic_approximation[cell->index()].compute(
2470 surrounding_points[point_index]);
2471 used_quadratic_approximation = true;
2472 }
2473 chart_points[point_index] =
2474 pull_back(cell, surrounding_points[point_index], guess);
2475
2476 // the initial guess may not have been good enough: if applicable,
2477 // try again with the affine approximation (which is more accurate
2478 // than the cheap methods used above)
2479 if (chart_points[point_index][0] ==
2481 !used_quadratic_approximation)
2482 {
2483 guess = quadratic_approximation[cell->index()].compute(
2484 surrounding_points[point_index]);
2485 chart_points[point_index] =
2486 pull_back(cell, surrounding_points[point_index], guess);
2487 }
2488
2489 if (chart_points[point_index][0] ==
2491 {
2492 for (unsigned int d = 0; d < dim; ++d)
2493 guess[d] = 0.5;
2494 chart_points[point_index] =
2495 pull_back(cell, surrounding_points[point_index], guess);
2496 }
2497 };
2498
2499 // check whether all points are inside the unit cell of the current chart
2500 for (unsigned int c = 0; c < nearby_cells.size(); ++c)
2501 {
2503 triangulation, level_coarse, nearby_cells[c]);
2504 bool inside_unit_cell = true;
2505 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
2506 {
2507 compute_chart_point(cell, i);
2508
2509 // Tolerance 5e-4 chosen that the method also works with manifolds
2510 // that have some discretization error like SphericalManifold
2511 if (GeometryInfo<dim>::is_inside_unit_cell(chart_points[i], 5e-4) ==
2512 false)
2513 {
2514 inside_unit_cell = false;
2515 break;
2516 }
2517 }
2518 if (inside_unit_cell == true)
2519 {
2520 return cell;
2521 }
2522
2523 // if we did not find a point and this was the last valid cell (the next
2524 // iterate being the end of the array or an invalid tag), we must stop
2525 if (c == nearby_cells.size() - 1 ||
2526 nearby_cells[c + 1] == numbers::invalid_unsigned_int)
2527 {
2528 // generate additional information to help debugging why we did not
2529 // get a point
2530 std::ostringstream message;
2531 for (unsigned int b = 0; b <= c; ++b)
2532 {
2534 triangulation, level_coarse, nearby_cells[b]);
2535 message << "Looking at cell " << cell->id()
2536 << " with vertices: " << std::endl;
2537 for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
2538 message << cell->vertex(v) << " ";
2539 message << std::endl;
2540 message << "Transformation to chart coordinates: " << std::endl;
2541 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
2542 {
2543 compute_chart_point(cell, i);
2544 message << surrounding_points[i] << " -> " << chart_points[i]
2545 << std::endl;
2546 }
2547 }
2548
2549 AssertThrow(false,
2551 message.str())));
2552 }
2553 }
2554
2555 // a valid inversion should have returned a point above. an invalid
2556 // inversion should have triggered the assertion, so we should never end up
2557 // here
2560}
2561
2562
2563
2564template <int dim, int spacedim>
2567 const ArrayView<const Point<spacedim>> &surrounding_points,
2568 const ArrayView<const double> &weights) const
2569{
2570 boost::container::small_vector<Point<dim>, 100> chart_points(
2571 surrounding_points.size());
2572 ArrayView<Point<dim>> chart_points_view =
2573 make_array_view(chart_points.begin(), chart_points.end());
2574 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2575
2576 const Point<dim> p_chart =
2577 chart_manifold.get_new_point(chart_points_view, weights);
2578
2579 return push_forward(cell, p_chart);
2580}
2581
2582
2583
2584template <int dim, int spacedim>
2585void
2587 const ArrayView<const Point<spacedim>> &surrounding_points,
2588 const Table<2, double> &weights,
2589 ArrayView<Point<spacedim>> new_points) const
2590{
2591 Assert(weights.size(0) > 0, ExcEmptyObject());
2592 AssertDimension(surrounding_points.size(), weights.size(1));
2593
2594 boost::container::small_vector<Point<dim>, 100> chart_points(
2595 surrounding_points.size());
2596 ArrayView<Point<dim>> chart_points_view =
2597 make_array_view(chart_points.begin(), chart_points.end());
2598 const auto cell = compute_chart_points(surrounding_points, chart_points_view);
2599
2600 boost::container::small_vector<Point<dim>, 100> new_points_on_chart(
2601 weights.size(0));
2602 chart_manifold.get_new_points(chart_points_view,
2603 weights,
2604 make_array_view(new_points_on_chart.begin(),
2605 new_points_on_chart.end()));
2606
2607 for (unsigned int row = 0; row < weights.size(0); ++row)
2608 new_points[row] = push_forward(cell, new_points_on_chart[row]);
2609}
2610
2611
2612
2613// explicit instantiations
2614#include "manifold_lib.inst"
2615
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
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition manifold.h:305
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
double inner_radius
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
double centerline_radius
TorusManifold(const double centerline_radius, const double inner_radius)
virtual std::unique_ptr< Manifold< dim, 3 > > clone() const override
DerivativeForm< 1, dim, spacedim > push_forward_gradient(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &chart_point, const Point< spacedim > &pushed_forward_chart_point) const
void initialize(const Triangulation< dim, spacedim > &triangulation)
Triangulation< dim, spacedim >::cell_iterator compute_chart_points(const ArrayView< const Point< spacedim > > &surrounding_points, ArrayView< Point< dim > > chart_points) const
std::array< unsigned int, 20 > get_possible_cells_around_points(const ArrayView< const Point< spacedim > > &surrounding_points) const
Point< dim > pull_back(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< spacedim > &p, const Point< dim > &initial_guess) const
virtual void get_new_points(const ArrayView< const Point< spacedim > > &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim > > new_points) const override
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Point< spacedim > push_forward(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const Point< dim > &chart_point) const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
cell_iterator begin(const unsigned int level=0) const
cell_iterator end() const
cell_iterator last() const
unsigned int n_cells() const
Signals signals
Definition tria.h:2513
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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
IteratorRange< active_cell_iterator > active_cell_iterators() const
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)
constexpr T fixed_power(const T t)
Definition utilities.h:942
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:259
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()
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 > &)