Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
manifold_lib.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2013 - 2022 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
32#include <boost/container/small_vector.hpp>
34
35#include <cmath>
36#include <memory>
37
39
40
41namespace internal
42{
43 // The pull_back function fails regularly in the compute_chart_points
44 // method, and, instead of throwing an exception, returns a point outside
45 // the unit cell. The individual coordinates of that point are given by the
46 // value below.
47 static constexpr double invalid_pull_back_coordinate = 20.0;
48
49 // Rotate a given unit vector u around the axis dir
50 // where the angle is given by the length of dir.
51 // This is the exponential map for a sphere.
54 {
55 const double theta = dir.norm();
56 if (theta < 1.e-10)
57 {
58 return u;
59 }
60 else
61 {
62 const Tensor<1, 3> dirUnit = dir / theta;
63 const Tensor<1, 3> tmp =
64 std::cos(theta) * u + std::sin(theta) * dirUnit;
65 return tmp / tmp.norm();
66 }
67 }
68
69 // Returns the direction to go from v to u
70 // projected to the plane perpendicular to the unit vector v.
71 // This one is more stable when u and v are nearly equal.
74 {
75 Tensor<1, 3> ans = u - v;
76 ans -= (ans * v) * v;
77 return ans; // ans = (u-v) - ((u-v)*v)*v
78 }
79
80 // helper function to compute a vector orthogonal to a given one.
81 // does nothing unless spacedim == 3.
82 template <int spacedim>
85 bool /*normalize*/ = false)
86 {
87 return {};
88 }
89
91 compute_normal(const Tensor<1, 3> &vector, bool normalize = false)
92 {
93 Assert(vector.norm_square() != 0.,
94 ExcMessage("The direction parameter must not be zero!"));
95 Point<3> normal;
96 if (std::abs(vector[0]) >= std::abs(vector[1]) &&
97 std::abs(vector[0]) >= std::abs(vector[2]))
98 {
99 normal[1] = -1.;
100 normal[2] = -1.;
101 normal[0] = (vector[1] + vector[2]) / vector[0];
102 }
103 else if (std::abs(vector[1]) >= std::abs(vector[0]) &&
104 std::abs(vector[1]) >= std::abs(vector[2]))
105 {
106 normal[0] = -1.;
107 normal[2] = -1.;
108 normal[1] = (vector[0] + vector[2]) / vector[1];
109 }
110 else
111 {
112 normal[0] = -1.;
113 normal[1] = -1.;
114 normal[2] = (vector[0] + vector[1]) / vector[2];
115 }
116 if (normalize)
117 normal /= normal.norm();
118 return normal;
119 }
120} // namespace internal
121
122
123
124// ============================================================
125// PolarManifold
126// ============================================================
127
128template <int dim, int spacedim>
130 : ChartManifold<dim, spacedim, spacedim>(
131 PolarManifold<dim, spacedim>::get_periodicity())
132 , center(center)
133{}
134
135
136
137template <int dim, int spacedim>
138std::unique_ptr<Manifold<dim, spacedim>>
140{
141 return std::make_unique<PolarManifold<dim, spacedim>>(*this);
142}
143
144
145
146template <int dim, int spacedim>
149{
150 Tensor<1, spacedim> periodicity;
151 // In two dimensions, theta is periodic.
152 // In three dimensions things are a little more complicated, since the only
153 // variable that is truly periodic is phi, while theta should be bounded
154 // between 0 and pi. There is currently no way to enforce this, so here we
155 // only fix periodicity for the last variable, corresponding to theta in 2d
156 // and phi in 3d.
157 periodicity[spacedim - 1] = 2 * numbers::PI;
158 return periodicity;
159}
160
161
162
163template <int dim, int spacedim>
166 const Point<spacedim> &spherical_point) const
167{
168 Assert(spherical_point[0] >= 0.0,
169 ExcMessage("Negative radius for given point."));
170 const double rho = spherical_point[0];
171 const double theta = spherical_point[1];
172
174 if (rho > 1e-10)
175 switch (spacedim)
176 {
177 case 2:
178 p[0] = rho * std::cos(theta);
179 p[1] = rho * std::sin(theta);
180 break;
181 case 3:
182 {
183 const double phi = spherical_point[2];
184 p[0] = rho * std::sin(theta) * std::cos(phi);
185 p[1] = rho * std::sin(theta) * std::sin(phi);
186 p[2] = rho * std::cos(theta);
187 break;
188 }
189 default:
190 Assert(false, ExcNotImplemented());
191 }
192 return p + center;
193}
194
195
196
197template <int dim, int spacedim>
200 const Point<spacedim> &space_point) const
201{
202 const Tensor<1, spacedim> R = space_point - center;
203 const double rho = R.norm();
204
206 p[0] = rho;
207
208 switch (spacedim)
209 {
210 case 2:
211 {
212 p[1] = std::atan2(R[1], R[0]);
213 if (p[1] < 0)
214 p[1] += 2 * numbers::PI;
215 break;
216 }
217
218 case 3:
219 {
220 const double z = R[2];
221 p[2] = std::atan2(R[1], R[0]); // phi
222 if (p[2] < 0)
223 p[2] += 2 * numbers::PI; // phi is periodic
224 p[1] = std::atan2(std::sqrt(R[0] * R[0] + R[1] * R[1]), z); // theta
225 break;
226 }
227
228 default:
229 Assert(false, ExcNotImplemented());
230 }
231 return p;
232}
233
234
235
236template <int dim, int spacedim>
239 const Point<spacedim> &spherical_point) const
240{
241 Assert(spherical_point[0] >= 0.0,
242 ExcMessage("Negative radius for given point."));
243 const double rho = spherical_point[0];
244 const double theta = spherical_point[1];
245
247 if (rho > 1e-10)
248 switch (spacedim)
249 {
250 case 2:
251 {
252 DX[0][0] = std::cos(theta);
253 DX[0][1] = -rho * std::sin(theta);
254 DX[1][0] = std::sin(theta);
255 DX[1][1] = rho * std::cos(theta);
256 break;
257 }
258
259 case 3:
260 {
261 const double phi = spherical_point[2];
262 DX[0][0] = std::sin(theta) * std::cos(phi);
263 DX[0][1] = rho * std::cos(theta) * std::cos(phi);
264 DX[0][2] = -rho * std::sin(theta) * std::sin(phi);
265
266 DX[1][0] = std::sin(theta) * std::sin(phi);
267 DX[1][1] = rho * std::cos(theta) * std::sin(phi);
268 DX[1][2] = rho * std::sin(theta) * std::cos(phi);
269
270 DX[2][0] = std::cos(theta);
271 DX[2][1] = -rho * std::sin(theta);
272 DX[2][2] = 0;
273 break;
274 }
275
276 default:
277 Assert(false, ExcNotImplemented());
278 }
279 return DX;
280}
281
282
283
284namespace
285{
286 template <int dim, int spacedim>
287 bool
288 spherical_face_is_horizontal(
290 const Point<spacedim> & manifold_center)
291 {
292 // We test whether a face is horizontal by checking that the vertices
293 // all have roughly the same distance from the center: If the
294 // maximum deviation for the distances from the vertices to the
295 // center is less than 1.e-5 of the distance between vertices (as
296 // measured by the minimum distance from any of the other vertices
297 // to the first vertex), then we call this a horizontal face.
298 constexpr unsigned int n_vertices =
300 std::array<double, n_vertices> sqr_distances_to_center;
301 std::array<double, n_vertices - 1> sqr_distances_to_first_vertex;
302 sqr_distances_to_center[0] =
303 (face->vertex(0) - manifold_center).norm_square();
304 for (unsigned int i = 1; i < n_vertices; ++i)
305 {
306 sqr_distances_to_center[i] =
307 (face->vertex(i) - manifold_center).norm_square();
308 sqr_distances_to_first_vertex[i - 1] =
309 (face->vertex(i) - face->vertex(0)).norm_square();
310 }
311 const auto minmax_sqr_distance =
312 std::minmax_element(sqr_distances_to_center.begin(),
313 sqr_distances_to_center.end());
314 const auto min_sqr_distance_to_first_vertex =
315 std::min_element(sqr_distances_to_first_vertex.begin(),
316 sqr_distances_to_first_vertex.end());
317
318 return (*minmax_sqr_distance.second - *minmax_sqr_distance.first <
319 1.e-10 * *min_sqr_distance_to_first_vertex);
320 }
321} // namespace
322
323
324
325template <int dim, int spacedim>
329 const Point<spacedim> & p) const
330{
331 // Let us first test whether we are on a "horizontal" face
332 // (tangential to the sphere). In this case, the normal vector is
333 // easy to compute since it is proportional to the vector from the
334 // center to the point 'p'.
335 if (spherical_face_is_horizontal<dim, spacedim>(face, center))
336 {
337 // So, if this is a "horizontal" face, then just compute the normal
338 // vector as the one from the center to the point 'p', adequately
339 // scaled.
340 const Tensor<1, spacedim> unnormalized_spherical_normal = p - center;
341 const Tensor<1, spacedim> normalized_spherical_normal =
342 unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
343 return normalized_spherical_normal;
344 }
345 else
346 // If it is not a horizontal face, just use the machinery of the
347 // base class.
349
350 return Tensor<1, spacedim>();
351}
352
353
354
355// ============================================================
356// SphericalManifold
357// ============================================================
358
359template <int dim, int spacedim>
362 : center(center)
363 , polar_manifold(center)
364{}
365
366
367
368template <int dim, int spacedim>
369std::unique_ptr<Manifold<dim, spacedim>>
371{
372 return std::make_unique<SphericalManifold<dim, spacedim>>(*this);
373}
374
375
376
377template <int dim, int spacedim>
380 const Point<spacedim> &p1,
381 const Point<spacedim> &p2,
382 const double w) const
383{
384 const double tol = 1e-10;
385
386 if ((p1 - p2).norm_square() < tol * tol || std::abs(w) < tol)
387 return p1;
388 else if (std::abs(w - 1.0) < tol)
389 return p2;
390
391 // If the points are one dimensional then there is no need for anything but
392 // a linear combination.
393 if (spacedim == 1)
394 return Point<spacedim>(w * p2 + (1 - w) * p1);
395
396 const Tensor<1, spacedim> v1 = p1 - center;
397 const Tensor<1, spacedim> v2 = p2 - center;
398 const double r1 = v1.norm();
399 const double r2 = v2.norm();
400
401 Assert(r1 > tol && r2 > tol,
402 ExcMessage("p1 and p2 cannot coincide with the center."));
403
404 const Tensor<1, spacedim> e1 = v1 / r1;
405 const Tensor<1, spacedim> e2 = v2 / r2;
406
407 // Find the cosine of the angle gamma described by v1 and v2.
408 const double cosgamma = e1 * e2;
409
410 // Points are collinear with the center (allow for 8*eps as a tolerance)
411 if (cosgamma < -1 + 8. * std::numeric_limits<double>::epsilon())
412 return center;
413
414 // Points are along a line, in which case e1 and e2 are essentially the same.
415 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
416 return Point<spacedim>(center + w * v2 + (1 - w) * v1);
417
418 // Find the angle sigma that corresponds to arclength equal to w. acos
419 // should never be undefined because we have ruled out the two special cases
420 // above.
421 const double sigma = w * std::acos(cosgamma);
422
423 // Normal to v1 in the plane described by v1,v2,and the origin.
424 // Since p1 and p2 do not coincide n is not zero and well defined.
425 Tensor<1, spacedim> n = v2 - (v2 * e1) * e1;
426 const double n_norm = n.norm();
427 Assert(n_norm > 0,
428 ExcInternalError("n should be different from the null vector. "
429 "Probably, this means v1==v2 or v2==0."));
430
431 n /= n_norm;
432
433 // Find the point Q along O,v1 such that
434 // P1,V,P2 has measure sigma.
435 const Tensor<1, spacedim> P = std::cos(sigma) * e1 + std::sin(sigma) * n;
436
437 // Project this point on the manifold.
438 return Point<spacedim>(center + (w * r2 + (1.0 - w) * r1) * P);
439}
440
441
442
443template <int dim, int spacedim>
446 const Point<spacedim> &p1,
447 const Point<spacedim> &p2) const
448{
449 const double tol = 1e-10;
450 (void)tol;
451
452 Assert(p1 != p2, ExcMessage("p1 and p2 should not concide."));
453
454 const Tensor<1, spacedim> v1 = p1 - center;
455 const Tensor<1, spacedim> v2 = p2 - center;
456 const double r1 = v1.norm();
457 const double r2 = v2.norm();
458
459 Assert(r1 > tol, ExcMessage("p1 cannot coincide with the center."));
460
461 Assert(r2 > tol, ExcMessage("p2 cannot coincide with the center."));
462
463 const Tensor<1, spacedim> e1 = v1 / r1;
464 const Tensor<1, spacedim> e2 = v2 / r2;
465
466 // Find the cosine of the angle gamma described by v1 and v2.
467 const double cosgamma = e1 * e2;
468
469 Assert(cosgamma > -1 + 8. * std::numeric_limits<double>::epsilon(),
470 ExcMessage("p1 and p2 cannot lie on the same diameter and be opposite "
471 "respect to the center."));
472
473 if (cosgamma > 1 - 8. * std::numeric_limits<double>::epsilon())
474 return v2 - v1;
475
476 // Normal to v1 in the plane described by v1,v2,and the origin.
477 // Since p1 and p2 do not coincide n is not zero and well defined.
478 Tensor<1, spacedim> n = v2 - (v2 * e1) * e1;
479 const double n_norm = n.norm();
480 Assert(n_norm > 0,
481 ExcInternalError("n should be different from the null vector. "
482 "Probably, this means v1==v2 or v2==0."));
483
484 n /= n_norm;
485
486 // this is the derivative of the geodesic in get_intermediate_point
487 // derived with respect to w and inserting w=0.
488 const double gamma = std::acos(cosgamma);
489 return (r2 - r1) * e1 + r1 * gamma * n;
490}
491
492
493
494template <int dim, int spacedim>
498 const Point<spacedim> & p) const
499{
500 // Let us first test whether we are on a "horizontal" face
501 // (tangential to the sphere). In this case, the normal vector is
502 // easy to compute since it is proportional to the vector from the
503 // center to the point 'p'.
504 if (spherical_face_is_horizontal<dim, spacedim>(face, center))
505 {
506 // So, if this is a "horizontal" face, then just compute the normal
507 // vector as the one from the center to the point 'p', adequately
508 // scaled.
509 const Tensor<1, spacedim> unnormalized_spherical_normal = p - center;
510 const Tensor<1, spacedim> normalized_spherical_normal =
511 unnormalized_spherical_normal / unnormalized_spherical_normal.norm();
512 return normalized_spherical_normal;
513 }
514 else
515 // If it is not a horizontal face, just use the machinery of the
516 // base class.
518
519 return Tensor<1, spacedim>();
520}
521
522
523
524template <>
525void
529{
530 Assert(false, ExcImpossibleInDim(1));
531}
532
533
534
535template <>
536void
540{
541 Assert(false, ExcImpossibleInDim(1));
542}
543
544
545
546template <int dim, int spacedim>
547void
550 typename Manifold<dim, spacedim>::FaceVertexNormals &face_vertex_normals)
551 const
552{
553 // Let us first test whether we are on a "horizontal" face
554 // (tangential to the sphere). In this case, the normal vector is
555 // easy to compute since it is proportional to the vector from the
556 // center to the point 'p'.
557 if (spherical_face_is_horizontal<dim, spacedim>(face, center))
558 {
559 // So, if this is a "horizontal" face, then just compute the normal
560 // vector as the one from the center to the point 'p', adequately
561 // scaled.
562 for (unsigned int vertex = 0;
563 vertex < GeometryInfo<spacedim>::vertices_per_face;
564 ++vertex)
565 face_vertex_normals[vertex] = face->vertex(vertex) - center;
566 }
567 else
568 Manifold<dim, spacedim>::get_normals_at_vertices(face, face_vertex_normals);
569}
570
571
572
573template <int dim, int spacedim>
574void
576 const ArrayView<const Point<spacedim>> &surrounding_points,
577 const Table<2, double> & weights,
578 ArrayView<Point<spacedim>> new_points) const
579{
580 AssertDimension(new_points.size(), weights.size(0));
581 AssertDimension(surrounding_points.size(), weights.size(1));
582
583 get_new_points(surrounding_points, make_array_view(weights), new_points);
584
585 return;
586}
587
588
589
590template <int dim, int spacedim>
593 const ArrayView<const Point<spacedim>> &vertices,
594 const ArrayView<const double> & weights) const
595{
596 // To avoid duplicating all of the logic in get_new_points, simply call it
597 // for one position.
598 Point<spacedim> new_point;
599 get_new_points(vertices,
600 weights,
601 make_array_view(&new_point, &new_point + 1));
602
603 return new_point;
604}
605
606
607
608template <int dim, int spacedim>
609void
611 const ArrayView<const Point<spacedim>> &surrounding_points,
612 const ArrayView<const double> & weights,
613 ArrayView<Point<spacedim>> new_points) const
614{
615 AssertDimension(weights.size(),
616 new_points.size() * surrounding_points.size());
617 const unsigned int weight_rows = new_points.size();
618 const unsigned int weight_columns = surrounding_points.size();
619
620 if (surrounding_points.size() == 2)
621 {
622 for (unsigned int row = 0; row < weight_rows; ++row)
623 new_points[row] =
624 get_intermediate_point(surrounding_points[0],
625 surrounding_points[1],
626 weights[row * weight_columns + 1]);
627 return;
628 }
629
630 boost::container::small_vector<std::pair<double, Tensor<1, spacedim>>, 100>
631 new_candidates(new_points.size());
632 boost::container::small_vector<Tensor<1, spacedim>, 100> directions(
633 surrounding_points.size(), Point<spacedim>());
634 boost::container::small_vector<double, 100> distances(
635 surrounding_points.size(), 0.0);
636 double max_distance = 0.;
637 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
638 {
639 directions[i] = surrounding_points[i] - center;
640 distances[i] = directions[i].norm();
641
642 if (distances[i] != 0.)
643 directions[i] /= distances[i];
644 else
645 Assert(false,
646 ExcMessage("One of the vertices coincides with the center. "
647 "This is not allowed!"));
648
649 // Check if an estimate is good enough,
650 // this is often the case for sufficiently refined meshes.
651 for (unsigned int k = 0; k < i; ++k)
652 {
653 const double squared_distance =
654 (directions[i] - directions[k]).norm_square();
655 max_distance = std::max(max_distance, squared_distance);
656 }
657 }
658
659 // Step 1: Check for some special cases, create simple linear guesses
660 // otherwise.
661 const double tolerance = 1e-10;
662 boost::container::small_vector<bool, 100> accurate_point_was_found(
663 new_points.size(), false);
664 const ArrayView<const Tensor<1, spacedim>> array_directions =
665 make_array_view(directions.begin(), directions.end());
666 const ArrayView<const double> array_distances =
667 make_array_view(distances.begin(), distances.end());
668 for (unsigned int row = 0; row < weight_rows; ++row)
669 {
670 new_candidates[row] =
671 guess_new_point(array_directions,
672 array_distances,
673 ArrayView<const double>(&weights[row * weight_columns],
674 weight_columns));
675
676 // If the candidate is the center, mark it as found to avoid entering
677 // the Newton iteration in step 2, which would crash.
678 if (new_candidates[row].first == 0.0)
679 {
680 new_points[row] = center;
681 accurate_point_was_found[row] = true;
682 continue;
683 }
684
685 // If not in 3D, just use the implementation from PolarManifold
686 // after we verified that the candidate is not the center.
687 if (spacedim < 3)
688 new_points[row] = polar_manifold.get_new_point(
689 surrounding_points,
690 ArrayView<const double>(&weights[row * weight_columns],
691 weight_columns));
692 }
693
694 // In this case, we treated the case that the candidate is the center and
695 // obtained the new locations from the PolarManifold object otherwise.
696 if (spacedim < 3)
697 return;
698
699 // If all the points are close to each other, we expect the estimate to
700 // be good enough. This tolerance was chosen such that the first iteration
701 // for a at least three time refined HyperShell mesh with radii .5 and 1.
702 // doesn't already succeed.
703 if (max_distance < 2e-2)
704 {
705 for (unsigned int row = 0; row < weight_rows; ++row)
706 new_points[row] =
707 center + new_candidates[row].first * new_candidates[row].second;
708
709 return;
710 }
711
712 // Step 2:
713 // Do more expensive Newton-style iterations to improve the estimate.
714
715 // Search for duplicate directions and merge them to minimize the cost of
716 // the get_new_point function call below.
717 boost::container::small_vector<double, 1000> merged_weights(weights.size());
718 boost::container::small_vector<Tensor<1, spacedim>, 100> merged_directions(
719 surrounding_points.size(), Point<spacedim>());
720 boost::container::small_vector<double, 100> merged_distances(
721 surrounding_points.size(), 0.0);
722
723 unsigned int n_unique_directions = 0;
724 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
725 {
726 bool found_duplicate = false;
727
728 // This inner loop is of @f$O(N^2)@f$ complexity, but
729 // surrounding_points.size() is usually at most 8 points large.
730 for (unsigned int j = 0; j < n_unique_directions; ++j)
731 {
732 const double squared_distance =
733 (directions[i] - directions[j]).norm_square();
734 if (!found_duplicate && squared_distance < 1e-28)
735 {
736 found_duplicate = true;
737 for (unsigned int row = 0; row < weight_rows; ++row)
738 merged_weights[row * weight_columns + j] +=
739 weights[row * weight_columns + i];
740 }
741 }
742
743 if (found_duplicate == false)
744 {
745 merged_directions[n_unique_directions] = directions[i];
746 merged_distances[n_unique_directions] = distances[i];
747 for (unsigned int row = 0; row < weight_rows; ++row)
748 merged_weights[row * weight_columns + n_unique_directions] =
749 weights[row * weight_columns + i];
750
751 ++n_unique_directions;
752 }
753 }
754
755 // Search for duplicate weight rows and merge them to minimize the cost of
756 // the get_new_point function call below.
757 boost::container::small_vector<unsigned int, 100> merged_weights_index(
758 new_points.size(), numbers::invalid_unsigned_int);
759 for (unsigned int row = 0; row < weight_rows; ++row)
760 {
761 for (unsigned int existing_row = 0; existing_row < row; ++existing_row)
762 {
763 bool identical_weights = true;
764
765 for (unsigned int weight_index = 0;
766 weight_index < n_unique_directions;
767 ++weight_index)
768 if (std::abs(merged_weights[row * weight_columns + weight_index] -
769 merged_weights[existing_row * weight_columns +
770 weight_index]) > tolerance)
771 {
772 identical_weights = false;
773 break;
774 }
775
776 if (identical_weights)
777 {
778 merged_weights_index[row] = existing_row;
779 break;
780 }
781 }
782 }
783
784 // Note that we only use the n_unique_directions first entries in the
785 // ArrayView
786 const ArrayView<const Tensor<1, spacedim>> array_merged_directions =
787 make_array_view(merged_directions.begin(),
788 merged_directions.begin() + n_unique_directions);
789 const ArrayView<const double> array_merged_distances =
790 make_array_view(merged_distances.begin(),
791 merged_distances.begin() + n_unique_directions);
792
793 for (unsigned int row = 0; row < weight_rows; ++row)
794 if (!accurate_point_was_found[row])
795 {
796 if (merged_weights_index[row] == numbers::invalid_unsigned_int)
797 {
798 const ArrayView<const double> array_merged_weights(
799 &merged_weights[row * weight_columns], n_unique_directions);
800 new_candidates[row].second =
801 get_new_point(array_merged_directions,
802 array_merged_distances,
803 array_merged_weights,
804 Point<spacedim>(new_candidates[row].second));
805 }
806 else
807 new_candidates[row].second =
808 new_candidates[merged_weights_index[row]].second;
809
810 new_points[row] =
811 center + new_candidates[row].first * new_candidates[row].second;
812 }
813}
814
815
816
817template <int dim, int spacedim>
818std::pair<double, Tensor<1, spacedim>>
820 const ArrayView<const Tensor<1, spacedim>> &directions,
821 const ArrayView<const double> & distances,
822 const ArrayView<const double> & weights) const
823{
824 const double tolerance = 1e-10;
825 double rho = 0.;
826 Tensor<1, spacedim> candidate;
827
828 // Perform a simple average ...
829 double total_weights = 0.;
830 for (unsigned int i = 0; i < directions.size(); ++i)
831 {
832 // if one weight is one, return its direction
833 if (std::abs(1 - weights[i]) < tolerance)
834 return std::make_pair(distances[i], directions[i]);
835
836 rho += distances[i] * weights[i];
837 candidate += directions[i] * weights[i];
838 total_weights += weights[i];
839 }
840
841 // ... and normalize if the candidate is different from the origin.
842 const double norm = candidate.norm();
843 if (norm == 0.)
844 return std::make_pair(0.0, Point<spacedim>());
845 candidate /= norm;
846 rho /= total_weights;
847
848 return std::make_pair(rho, candidate);
849}
850
851
852namespace
853{
854 template <int spacedim>
856 do_get_new_point(const ArrayView<const Tensor<1, spacedim>> & /*directions*/,
857 const ArrayView<const double> & /*distances*/,
858 const ArrayView<const double> & /*weights*/,
859 const Point<spacedim> & /*candidate_point*/)
860 {
861 Assert(false, ExcNotImplemented());
862 return {};
863 }
864
865 template <>
867 do_get_new_point(const ArrayView<const Tensor<1, 3>> &directions,
868 const ArrayView<const double> & distances,
869 const ArrayView<const double> & weights,
870 const Point<3> & candidate_point)
871 {
872 (void)distances;
873
874 AssertDimension(directions.size(), distances.size());
875 AssertDimension(directions.size(), weights.size());
876
877 Point<3> candidate = candidate_point;
878 const unsigned int n_merged_points = directions.size();
879 const double tolerance = 1e-10;
880 const int max_iterations = 10;
881
882 {
883 // If the candidate happens to coincide with a normalized
884 // direction, we return it. Otherwise, the Hessian would be singular.
885 for (unsigned int i = 0; i < n_merged_points; ++i)
886 {
887 const double squared_distance =
888 (candidate - directions[i]).norm_square();
889 if (squared_distance < tolerance * tolerance)
890 return candidate;
891 }
892
893 // check if we only have two points now, in which case we can use the
894 // get_intermediate_point function
895 if (n_merged_points == 2)
896 {
897 SphericalManifold<3, 3> unit_manifold;
898 Assert(std::abs(weights[0] + weights[1] - 1.0) < 1e-13,
899 ExcMessage("Weights do not sum up to 1"));
900 Point<3> intermediate =
901 unit_manifold.get_intermediate_point(Point<3>(directions[0]),
902 Point<3>(directions[1]),
903 weights[1]);
904 return intermediate;
905 }
906
907 Tensor<1, 3> vPerp;
908 Tensor<2, 2> Hessian;
910 Tensor<1, 2> gradlocal;
911
912 // On success we exit the loop early.
913 // Otherwise, we just take the result after max_iterations steps.
914 for (unsigned int i = 0; i < max_iterations; ++i)
915 {
916 // Step 2a: Find new descent direction
917
918 // Get local basis for the estimate candidate
919 const Tensor<1, 3> Clocalx = internal::compute_normal(candidate);
920 const Tensor<1, 3> Clocaly = cross_product_3d(candidate, Clocalx);
921
922 // For each vertices vector, compute the tangent vector from candidate
923 // towards the vertices vector -- its length is the spherical length
924 // from candidate to the vertices vector.
925 // Then compute its contribution to the Hessian.
926 gradient = 0.;
927 Hessian = 0.;
928 for (unsigned int i = 0; i < n_merged_points; ++i)
929 if (std::abs(weights[i]) > 1.e-15)
930 {
931 vPerp = internal::projected_direction(directions[i], candidate);
932 const double sinthetaSq = vPerp.norm_square();
933 const double sintheta = std::sqrt(sinthetaSq);
934 if (sintheta < tolerance)
935 {
936 Hessian[0][0] += weights[i];
937 Hessian[1][1] += weights[i];
938 }
939 else
940 {
941 const double costheta = (directions[i]) * candidate;
942 const double theta = std::atan2(sintheta, costheta);
943 const double sincthetaInv = theta / sintheta;
944
945 const double cosphi = vPerp * Clocalx;
946 const double sinphi = vPerp * Clocaly;
947
948 gradlocal[0] = cosphi;
949 gradlocal[1] = sinphi;
950 gradient += (weights[i] * sincthetaInv) * gradlocal;
951
952 const double wt = weights[i] / sinthetaSq;
953 const double sinphiSq = sinphi * sinphi;
954 const double cosphiSq = cosphi * cosphi;
955 const double tt = sincthetaInv * costheta;
956 const double offdiag = cosphi * sinphi * wt * (1.0 - tt);
957 Hessian[0][0] += wt * (cosphiSq + tt * sinphiSq);
958 Hessian[0][1] += offdiag;
959 Hessian[1][0] += offdiag;
960 Hessian[1][1] += wt * (sinphiSq + tt * cosphiSq);
961 }
962 }
963
964 Assert(determinant(Hessian) > tolerance, ExcInternalError());
965
966 const Tensor<2, 2> inverse_Hessian = invert(Hessian);
967
968 const Tensor<1, 2> xDisplocal = inverse_Hessian * gradient;
969 const Tensor<1, 3> xDisp =
970 xDisplocal[0] * Clocalx + xDisplocal[1] * Clocaly;
971
972 // Step 2b: rotate candidate in direction xDisp for a new candidate.
973 const Point<3> candidateOld = candidate;
974 candidate =
975 Point<3>(internal::apply_exponential_map(candidate, xDisp));
976
977 // Step 2c: return the new candidate if we didn't move
978 if ((candidate - candidateOld).norm_square() < tolerance * tolerance)
979 break;
980 }
981 }
982 return candidate;
983 }
984} // namespace
985
986
987
988template <int dim, int spacedim>
991 const ArrayView<const Tensor<1, spacedim>> &,
994 const Point<spacedim> &) const
995{
996 Assert(false, ExcNotImplemented());
997 return {};
998}
999
1000
1001
1002template <>
1005 const ArrayView<const Tensor<1, 3>> &directions,
1006 const ArrayView<const double> & distances,
1007 const ArrayView<const double> & weights,
1008 const Point<3> & candidate_point) const
1009{
1010 return do_get_new_point(directions, distances, weights, candidate_point);
1011}
1012
1013
1014
1015template <>
1018 const ArrayView<const Tensor<1, 3>> &directions,
1019 const ArrayView<const double> & distances,
1020 const ArrayView<const double> & weights,
1021 const Point<3> & candidate_point) const
1022{
1023 return do_get_new_point(directions, distances, weights, candidate_point);
1024}
1025
1026
1027
1028template <>
1031 const ArrayView<const Tensor<1, 3>> &directions,
1032 const ArrayView<const double> & distances,
1033 const ArrayView<const double> & weights,
1034 const Point<3> & candidate_point) const
1035{
1036 return do_get_new_point(directions, distances, weights, candidate_point);
1037}
1038
1039
1040
1041// ============================================================
1042// CylindricalManifold
1043// ============================================================
1044template <int dim, int spacedim>
1046 const double tolerance)
1047 : CylindricalManifold<dim, spacedim>(Point<spacedim>::unit_vector(axis),
1048 Point<spacedim>(),
1049 tolerance)
1050{
1051 // do not use static_assert to make dimension-independent programming
1052 // easier.
1053 Assert(spacedim == 3,
1054 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1055}
1056
1057
1058
1059template <int dim, int spacedim>
1061 const Tensor<1, spacedim> &direction,
1062 const Point<spacedim> & point_on_axis,
1063 const double tolerance)
1064 : ChartManifold<dim, spacedim, 3>(Tensor<1, 3>({0, 2. * numbers::PI, 0}))
1065 , normal_direction(internal::compute_normal(direction, true))
1066 , direction(direction / direction.norm())
1067 , point_on_axis(point_on_axis)
1068 , tolerance(tolerance)
1069 , dxn(cross_product_3d(this->direction, normal_direction))
1070{
1071 // do not use static_assert to make dimension-independent programming
1072 // easier.
1073 Assert(spacedim == 3,
1074 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1075}
1076
1077
1078
1079template <int dim, int spacedim>
1080std::unique_ptr<Manifold<dim, spacedim>>
1082{
1083 return std::make_unique<CylindricalManifold<dim, spacedim>>(*this);
1084}
1085
1086
1087
1088template <int dim, int spacedim>
1091 const ArrayView<const Point<spacedim>> &surrounding_points,
1092 const ArrayView<const double> & weights) const
1093{
1094 Assert(spacedim == 3,
1095 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1096
1097 // First check if the average in space lies on the axis.
1098 Point<spacedim> middle;
1099 double average_length = 0.;
1100 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
1101 {
1102 middle += surrounding_points[i] * weights[i];
1103 average_length += surrounding_points[i].square() * weights[i];
1104 }
1105 middle -= point_on_axis;
1106 const double lambda = middle * direction;
1107
1108 if ((middle - direction * lambda).square() < tolerance * average_length)
1109 return point_on_axis + direction * lambda;
1110 else // If not, using the ChartManifold should yield valid results.
1111 return ChartManifold<dim, spacedim, 3>::get_new_point(surrounding_points,
1112 weights);
1113}
1114
1115
1116
1117template <int dim, int spacedim>
1120 const Point<spacedim> &space_point) const
1121{
1122 Assert(spacedim == 3,
1123 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1124
1125 // First find the projection of the given point to the axis.
1126 const Tensor<1, spacedim> normalized_point = space_point - point_on_axis;
1127 const double lambda = normalized_point * direction;
1128 const Point<spacedim> projection = point_on_axis + direction * lambda;
1129 const Tensor<1, spacedim> p_diff = space_point - projection;
1130
1131 // Then compute the angle between the projection direction and
1132 // another vector orthogonal to the direction vector.
1133 const double phi = Physics::VectorRelations::signed_angle(normal_direction,
1134 p_diff,
1135 /*axis=*/direction);
1136
1137 // Return distance from the axis, angle and signed distance on the axis.
1138 return Point<3>(p_diff.norm(), phi, lambda);
1139}
1140
1141
1142
1143template <int dim, int spacedim>
1146 const Point<3> &chart_point) const
1147{
1148 Assert(spacedim == 3,
1149 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1150
1151 // Rotate the orthogonal direction by the given angle
1152 const double sine_r = std::sin(chart_point(1)) * chart_point(0);
1153 const double cosine_r = std::cos(chart_point(1)) * chart_point(0);
1154 const Tensor<1, spacedim> intermediate =
1155 normal_direction * cosine_r + dxn * sine_r;
1156
1157 // Finally, put everything together.
1158 return point_on_axis + direction * chart_point(2) + intermediate;
1159}
1160
1161
1162
1163template <int dim, int spacedim>
1166 const Point<3> &chart_point) const
1167{
1168 Assert(spacedim == 3,
1169 ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
1170
1171 Tensor<2, 3> derivatives;
1172
1173 // Rotate the orthogonal direction by the given angle
1174 const double sine = std::sin(chart_point(1));
1175 const double cosine = std::cos(chart_point(1));
1176 const Tensor<1, spacedim> intermediate =
1177 normal_direction * cosine + dxn * sine;
1178
1179 // avoid compiler warnings
1180 constexpr int s0 = 0 % spacedim;
1181 constexpr int s1 = 1 % spacedim;
1182 constexpr int s2 = 2 % spacedim;
1183
1184 // derivative w.r.t the radius
1185 derivatives[s0][s0] = intermediate[s0];
1186 derivatives[s1][s0] = intermediate[s1];
1187 derivatives[s2][s0] = intermediate[s2];
1188
1189 // derivatives w.r.t the angle
1190 derivatives[s0][s1] = -normal_direction[s0] * sine + dxn[s0] * cosine;
1191 derivatives[s1][s1] = -normal_direction[s1] * sine + dxn[s1] * cosine;
1192 derivatives[s2][s1] = -normal_direction[s2] * sine + dxn[s2] * cosine;
1193
1194 // derivatives w.r.t the direction of the axis
1195 derivatives[s0][s2] = direction[s0];
1196 derivatives[s1][s2] = direction[s1];
1197 derivatives[s2][s2] = direction[s2];
1198
1199 return derivatives;
1200}
1201
1202
1203
1204// ============================================================
1205// EllipticalManifold
1206// ============================================================
1207template <int dim, int spacedim>
1209 const Point<spacedim> & center,
1210 const Tensor<1, spacedim> &major_axis_direction,
1211 const double eccentricity)
1212 : ChartManifold<dim, spacedim, spacedim>(
1213 EllipticalManifold<dim, spacedim>::get_periodicity())
1214 , direction(major_axis_direction)
1215 , center(center)
1216 , cosh_u(1.0 / eccentricity)
1217 , sinh_u(std::sqrt(cosh_u * cosh_u - 1.0))
1218{
1219 // Throws an exception if dim!=2 || spacedim!=2.
1220 Assert(dim == 2 && spacedim == 2, ExcNotImplemented());
1221 // Throws an exception if eccentricity is not in range.
1222 Assert(std::signbit(cosh_u * cosh_u - 1.0) == false,
1223 ExcMessage(
1224 "Invalid eccentricity: It must satisfy 0 < eccentricity < 1."));
1225 const double direction_norm = direction.norm();
1226 Assert(direction_norm != 0.0,
1227 ExcMessage(
1228 "Invalid major axis direction vector: Null vector not allowed."));
1229 direction /= direction_norm;
1230}
1231
1232
1233
1234template <int dim, int spacedim>
1235std::unique_ptr<Manifold<dim, spacedim>>
1237{
1238 return std::make_unique<EllipticalManifold<dim, spacedim>>(*this);
1239}
1240
1241
1242
1243template <int dim, int spacedim>
1246{
1247 Tensor<1, spacedim> periodicity;
1248 // The second elliptical coordinate is periodic, while the first is not.
1249 // Enforce periodicity on the last variable.
1250 periodicity[spacedim - 1] = 2.0 * numbers::PI;
1251 return periodicity;
1252}
1253
1254
1255
1256template <int dim, int spacedim>
1259{
1260 Assert(false, ExcNotImplemented());
1261 return {};
1262}
1263
1264
1265
1266template <>
1269{
1270 const double cs = std::cos(chart_point[1]);
1271 const double sn = std::sin(chart_point[1]);
1272 // Coordinates in the reference frame (i.e. major axis direction is
1273 // x-axis)
1274 const double x = chart_point[0] * cosh_u * cs;
1275 const double y = chart_point[0] * sinh_u * sn;
1276 // Rotate them according to the major axis direction
1277 const Point<2> p(direction[0] * x - direction[1] * y,
1278 direction[1] * x + direction[0] * y);
1279 return p + center;
1280}
1281
1282
1283
1284template <int dim, int spacedim>
1287{
1288 Assert(false, ExcNotImplemented());
1289 return {};
1290}
1291
1292
1293
1294template <>
1297{
1298 // Moving space_point in the reference coordinate system.
1299 const double x0 = space_point[0] - center[0];
1300 const double y0 = space_point[1] - center[1];
1301 const double x = direction[0] * x0 + direction[1] * y0;
1302 const double y = -direction[1] * x0 + direction[0] * y0;
1303 const double pt0 =
1304 std::sqrt((x * x) / (cosh_u * cosh_u) + (y * y) / (sinh_u * sinh_u));
1305 // If the radius is exactly zero, the point coincides with the origin.
1306 if (pt0 == 0.0)
1307 {
1308 return center;
1309 }
1310 double cos_eta = x / (pt0 * cosh_u);
1311 if (cos_eta < -1.0)
1312 {
1313 cos_eta = -1.0;
1314 }
1315 if (cos_eta > 1.0)
1316 {
1317 cos_eta = 1.0;
1318 }
1319 const double eta = std::acos(cos_eta);
1320 const double pt1 = (std::signbit(y) ? 2.0 * numbers::PI - eta : eta);
1321 return {pt0, pt1};
1322}
1323
1324
1325
1326template <int dim, int spacedim>
1329 const Point<spacedim> &) const
1330{
1331 Assert(false, ExcNotImplemented());
1332 return {};
1333}
1334
1335
1336
1337template <>
1340 const Point<2> &chart_point) const
1341{
1342 const double cs = std::cos(chart_point[1]);
1343 const double sn = std::sin(chart_point[1]);
1344 Tensor<2, 2> dX;
1345 dX[0][0] = cosh_u * cs;
1346 dX[0][1] = -chart_point[0] * cosh_u * sn;
1347 dX[1][0] = sinh_u * sn;
1348 dX[1][1] = chart_point[0] * sinh_u * cs;
1349
1350 // rotate according to the major axis direction
1352 {{+direction[0], -direction[1]}, {direction[1], direction[0]}}};
1353
1354 return rot * dX;
1355}
1356
1357
1358
1359// ============================================================
1360// FunctionManifold
1361// ============================================================
1362template <int dim, int spacedim, int chartdim>
1364 const Function<chartdim> & push_forward_function,
1365 const Function<spacedim> & pull_back_function,
1366 const Tensor<1, chartdim> &periodicity,
1367 const double tolerance)
1368 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1369 , const_map()
1370 , push_forward_function(&push_forward_function)
1371 , pull_back_function(&pull_back_function)
1372 , tolerance(tolerance)
1373 , owns_pointers(false)
1374 , finite_difference_step(0)
1375{
1376 AssertDimension(push_forward_function.n_components, spacedim);
1377 AssertDimension(pull_back_function.n_components, chartdim);
1378}
1379
1380
1381
1382template <int dim, int spacedim, int chartdim>
1384 std::unique_ptr<Function<chartdim>> push_forward,
1385 std::unique_ptr<Function<spacedim>> pull_back,
1386 const Tensor<1, chartdim> & periodicity,
1387 const double tolerance)
1388 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1389 , const_map()
1390 , push_forward_function(push_forward.release())
1391 , pull_back_function(pull_back.release())
1392 , tolerance(tolerance)
1393 , owns_pointers(true)
1394 , finite_difference_step(0)
1395{
1396 AssertDimension(push_forward_function->n_components, spacedim);
1397 AssertDimension(pull_back_function->n_components, chartdim);
1398}
1399
1400
1401
1402template <int dim, int spacedim, int chartdim>
1404 const std::string push_forward_expression,
1405 const std::string pull_back_expression,
1406 const Tensor<1, chartdim> & periodicity,
1407 const typename FunctionParser<spacedim>::ConstMap const_map,
1408 const std::string chart_vars,
1409 const std::string space_vars,
1410 const double tolerance,
1411 const double h)
1412 : ChartManifold<dim, spacedim, chartdim>(periodicity)
1413 , const_map(const_map)
1414 , tolerance(tolerance)
1415 , owns_pointers(true)
1416 , push_forward_expression(push_forward_expression)
1417 , pull_back_expression(pull_back_expression)
1418 , chart_vars(chart_vars)
1419 , space_vars(space_vars)
1420 , finite_difference_step(h)
1421{
1422 FunctionParser<chartdim> *pf = new FunctionParser<chartdim>(spacedim, 0.0, h);
1423 FunctionParser<spacedim> *pb = new FunctionParser<spacedim>(chartdim, 0.0, h);
1427 pull_back_function = pb;
1428}
1429
1430
1431
1432template <int dim, int spacedim, int chartdim>
1434{
1435 if (owns_pointers == true)
1436 {
1437 const Function<chartdim> *pf = push_forward_function;
1438 push_forward_function = nullptr;
1439 delete pf;
1440
1441 const Function<spacedim> *pb = pull_back_function;
1442 pull_back_function = nullptr;
1443 delete pb;
1444 }
1445}
1446
1447
1448
1449template <int dim, int spacedim, int chartdim>
1450std::unique_ptr<Manifold<dim, spacedim>>
1452{
1453 // This manifold can be constructed either by providing an expression for the
1454 // push forward and the pull back charts, or by providing two Function
1455 // objects. In the first case, the push_forward and pull_back functions are
1456 // created internally in FunctionManifold, and destroyed when this object is
1457 // deleted. In the second case, the function objects are destroyed if they
1458 // are passed as pointers upon construction.
1459 // We need to make sure that our cloned object is constructed in the
1460 // same way this class was constructed, and that its internal Function
1461 // pointers point either to the same Function objects used to construct this
1462 // function or that the newly generated manifold creates internally the
1463 // push_forward and pull_back functions using the same expressions that were
1464 // used to construct this class.
1465 if (!(push_forward_expression.empty() && pull_back_expression.empty()))
1466 {
1467 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1468 push_forward_expression,
1469 pull_back_expression,
1470 this->get_periodicity(),
1471 const_map,
1472 chart_vars,
1473 space_vars,
1474 tolerance,
1475 finite_difference_step);
1476 }
1477 else
1478 {
1479 return std::make_unique<FunctionManifold<dim, spacedim, chartdim>>(
1480 *push_forward_function,
1481 *pull_back_function,
1482 this->get_periodicity(),
1483 tolerance);
1484 }
1485}
1486
1487
1488
1489template <int dim, int spacedim, int chartdim>
1492 const Point<chartdim> &chart_point) const
1493{
1494 Vector<double> pf(spacedim);
1495 Point<spacedim> result;
1496 push_forward_function->vector_value(chart_point, pf);
1497 for (unsigned int i = 0; i < spacedim; ++i)
1498 result[i] = pf[i];
1499
1500#ifdef DEBUG
1501 Vector<double> pb(chartdim);
1502 pull_back_function->vector_value(result, pb);
1503 for (unsigned int i = 0; i < chartdim; ++i)
1504 Assert(
1505 (chart_point.norm() > tolerance &&
1506 (std::abs(pb[i] - chart_point[i]) < tolerance * chart_point.norm())) ||
1507 (std::abs(pb[i] - chart_point[i]) < tolerance),
1508 ExcMessage(
1509 "The push forward is not the inverse of the pull back! Bailing out."));
1510#endif
1511
1512 return result;
1513}
1514
1515
1516
1517template <int dim, int spacedim, int chartdim>
1520 const Point<chartdim> &chart_point) const
1521{
1523 for (unsigned int i = 0; i < spacedim; ++i)
1524 {
1525 const auto gradient = push_forward_function->gradient(chart_point, i);
1526 for (unsigned int j = 0; j < chartdim; ++j)
1527 DF[i][j] = gradient[j];
1528 }
1529 return DF;
1530}
1531
1532
1533
1534template <int dim, int spacedim, int chartdim>
1537 const Point<spacedim> &space_point) const
1538{
1539 Vector<double> pb(chartdim);
1540 Point<chartdim> result;
1541 pull_back_function->vector_value(space_point, pb);
1542 for (unsigned int i = 0; i < chartdim; ++i)
1543 result[i] = pb[i];
1544 return result;
1545}
1546
1547
1548
1549// ============================================================
1550// TorusManifold
1551// ============================================================
1552template <int dim>
1555{
1556 double x = p(0);
1557 double z = p(1);
1558 double y = p(2);
1559 double phi = std::atan2(y, x);
1560 double theta = std::atan2(z, std::sqrt(x * x + y * y) - R);
1561 double w = std::sqrt(std::pow(y - std::sin(phi) * R, 2.0) +
1562 std::pow(x - std::cos(phi) * R, 2.0) + 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 QIteraded
1693 // quadrature with no subdivisions.
1694 std::vector<Point<dim>> unit_points =
1695 QIterated<dim>(QTrapez<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:699
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:999
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
Definition: manifold_lib.h:557
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
Definition: manifold_lib.h:565
SmartPointer< const Function< spacedim >, FunctionManifold< dim, spacedim, chartdim > > pull_back_function
Definition: manifold_lib.h:728
const FunctionParser< spacedim >::ConstMap const_map
Definition: manifold_lib.h:714
const std::string chart_vars
Definition: manifold_lib.h:759
const std::string pull_back_expression
Definition: manifold_lib.h:754
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
Definition: manifold_lib.h:721
const std::string space_vars
Definition: manifold_lib.h:764
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
Definition: manifold_lib.h:749
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:311
Definition: point.h:111
PolarManifold(const Point< spacedim > center=Point< spacedim >())
static Tensor< 1, spacedim > get_periodicity()
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
virtual DerivativeForm< 1, spacedim, spacedim > push_forward_gradient(const Point< spacedim > &chart_point) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
virtual Point< spacedim > push_forward(const Point< spacedim > &chart_point) const override
virtual Point< spacedim > pull_back(const Point< spacedim > &space_point) const override
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
Definition: tensor.h:503
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()
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
Point< 3 > center
Point< 3 > vertices[4]
Point< 2 > second
Definition: grid_out.cc:4604
Point< 2 > first
Definition: grid_out.cc:4603
const unsigned int v1
Definition: grid_tools.cc:1000
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition: tria.h:1355
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)
Definition: manifold_lib.cc:73
Tensor< 1, 3 > apply_exponential_map(const Tensor< 1, 3 > &u, const Tensor< 1, 3 > &dir)
Definition: manifold_lib.cc:53
static constexpr double invalid_pull_back_coordinate
Definition: manifold_lib.cc:47
Point< spacedim > compute_normal(const Tensor< 1, spacedim > &, bool=false)
Definition: manifold_lib.cc:84
static constexpr double PI
Definition: numbers.h:233
static const unsigned int invalid_unsigned_int
Definition: types.h:201
const types::manifold_id flat_manifold_id
Definition: types.h:269
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 > pow(const ::VectorizedArray< Number, width > &, const Number p)
::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)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2660
const ::Triangulation< dim, spacedim > & tria