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