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