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