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