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