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