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