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