Reference documentation for deal.II version Git f9e0250abc 2021-03-01 16:46:02 +0100
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
manifold.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/table.h>
17 #include <deal.II/base/tensor.h>
18 
19 #include <deal.II/fe/fe_q.h>
20 
21 #include <deal.II/grid/manifold.h>
22 #include <deal.II/grid/tria.h>
25 
27 #include <boost/container/small_vector.hpp>
29 
30 #include <cmath>
31 #include <memory>
32 
34 
35 using namespace Manifolds;
36 
37 /* -------------------------- Manifold --------------------- */
38 template <int dim, int spacedim>
41  const ArrayView<const Point<spacedim>> &,
42  const Point<spacedim> &) const
43 {
44  Assert(false, ExcPureFunctionCalled());
45  return {};
46 }
47 
48 
49 
50 template <int dim, int spacedim>
53  const Point<spacedim> &p2,
54  const double w) const
55 {
56  const std::array<Point<spacedim>, 2> vertices{{p1, p2}};
57  return project_to_manifold(make_array_view(vertices.begin(), vertices.end()),
58  w * p2 + (1 - w) * p1);
59 }
60 
61 
62 
63 template <int dim, int spacedim>
66  const ArrayView<const Point<spacedim>> &surrounding_points,
67  const ArrayView<const double> & weights) const
68 {
69  const double tol = 1e-10;
70  const unsigned int n_points = surrounding_points.size();
71 
72  Assert(n_points > 0, ExcMessage("There should be at least one point."));
73 
74  Assert(n_points == weights.size(),
75  ExcMessage(
76  "There should be as many surrounding points as weights given."));
77 
78  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
79  tol,
80  ExcMessage("The weights for the individual points should sum to 1!"));
81 
82  // First sort points in the order of their weights. This is done to
83  // produce unique points even if get_intermediate_points is not
84  // associative (as for the SphericalManifold).
85  boost::container::small_vector<unsigned int, 100> permutation(n_points);
86  std::iota(permutation.begin(), permutation.end(), 0u);
87  std::sort(permutation.begin(),
88  permutation.end(),
89  [&weights](const std::size_t a, const std::size_t b) {
90  return weights[a] < weights[b];
91  });
92 
93  // Now loop over points in the order of their associated weight
94  Point<spacedim> p = surrounding_points[permutation[0]];
95  double w = weights[permutation[0]];
96 
97  for (unsigned int i = 1; i < n_points; ++i)
98  {
99  double weight = 0.0;
100  if (std::abs(weights[permutation[i]] + w) < tol)
101  weight = 0.0;
102  else
103  weight = w / (weights[permutation[i]] + w);
104 
105  if (std::abs(weight) > 1e-14)
106  {
107  p = get_intermediate_point(p,
108  surrounding_points[permutation[i]],
109  1.0 - weight);
110  }
111  else
112  {
113  p = surrounding_points[permutation[i]];
114  }
115  w += weights[permutation[i]];
116  }
117 
118  return p;
119 }
120 
121 
122 
123 template <int dim, int spacedim>
124 void
126  const ArrayView<const Point<spacedim>> &surrounding_points,
127  const Table<2, double> & weights,
128  ArrayView<Point<spacedim>> new_points) const
129 {
130  AssertDimension(surrounding_points.size(), weights.size(1));
131 
132  for (unsigned int row = 0; row < weights.size(0); ++row)
133  {
134  new_points[row] =
135  get_new_point(make_array_view(surrounding_points.begin(),
136  surrounding_points.end()),
137  make_array_view(weights, row));
138  }
139 }
140 
141 
142 
143 template <>
146  const Point<2> & p) const
147 {
148  const int spacedim = 2;
149 
150  // get the tangent vector from the point 'p' in the direction of the further
151  // one of the two vertices that make up the face of this 2d cell
152  const Tensor<1, spacedim> tangent =
153  ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
154  -get_tangent_vector(p, face->vertex(0)) :
155  get_tangent_vector(p, face->vertex(1)));
156 
157  // then rotate it by 90 degrees
158  const Tensor<1, spacedim> normal = cross_product_2d(tangent);
159  return normal / normal.norm();
160 }
161 
162 
163 
164 template <>
167  const Point<3> & p) const
168 {
169  const int spacedim = 3;
170 
171  const std::array<Point<spacedim>, 4> vertices{
172  {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
173  const std::array<double, 4> distances{{vertices[0].distance(p),
174  vertices[1].distance(p),
175  vertices[2].distance(p),
176  vertices[3].distance(p)}};
177  const double max_distance = std::max(std::max(distances[0], distances[1]),
178  std::max(distances[2], distances[3]));
179 
180  // We need to find two tangential vectors to the given point p, but we do
181  // not know how the point is oriented against the face. We guess the two
182  // directions by assuming a flat topology and take the two directions that
183  // indicate the angle closest to a perpendicular one (i.e., cos(theta) close
184  // to zero). We start with an invalid value but the loops should always find
185  // a value.
186  double abs_cos_angle = std::numeric_limits<double>::max();
187  unsigned int first_index = numbers::invalid_unsigned_int,
188  second_index = numbers::invalid_unsigned_int;
189  for (unsigned int i = 0; i < 3; ++i)
190  if (distances[i] > 1e-8 * max_distance)
191  for (unsigned int j = i + 1; j < 4; ++j)
192  if (distances[j] > 1e-8 * max_distance)
193  {
194  const double new_angle = (p - vertices[i]) * (p - vertices[j]) /
195  (distances[i] * distances[j]);
196  // multiply by factor 0.999 to bias the search in a way that
197  // avoids trouble with roundoff
198  if (std::abs(new_angle) < 0.999 * abs_cos_angle)
199  {
200  abs_cos_angle = std::abs(new_angle);
201  first_index = i;
202  second_index = j;
203  }
204  }
205  Assert(first_index != numbers::invalid_unsigned_int,
206  ExcMessage("The search for possible directions did not succeed."));
207 
208  // Compute tangents and normal for selected vertices
209  Tensor<1, spacedim> t1 = get_tangent_vector(p, vertices[first_index]);
210  Tensor<1, spacedim> t2 = get_tangent_vector(p, vertices[second_index]);
211  Tensor<1, spacedim> normal = cross_product_3d(t1, t2);
212 
213  Assert(
215  t1.norm_square() * t2.norm_square(),
216  ExcMessage(
217  "Manifold::normal_vector was unable to find a suitable combination "
218  "of vertices to compute a normal on this face. We chose the secants "
219  "that are as orthogonal as possible, but tangents appear to be "
220  "linearly dependent. Check for distorted faces in your triangulation."));
221 
222  // Now figure out if we need to flip the direction, we do this by comparing
223  // to a reference normal that would be the correct result if all vertices
224  // would lie in a plane
225  const Tensor<1, spacedim> rt1 = vertices[3] - vertices[0];
226  const Tensor<1, spacedim> rt2 = vertices[2] - vertices[1];
227  const Tensor<1, spacedim> reference_normal = cross_product_3d(rt1, rt2);
228 
229  if (reference_normal * normal < 0.0)
230  normal *= -1.0;
231 
232  return normal / normal.norm();
233 }
234 
235 
236 
237 template <int dim, int spacedim>
240  const typename Triangulation<dim, spacedim>::face_iterator & /*face*/,
241  const Point<spacedim> & /*p*/) const
242 {
243  Assert(false, ExcPureFunctionCalled());
244  return Tensor<1, spacedim>();
245 }
246 
247 
248 
249 template <>
250 void
253  FaceVertexNormals & n) const
254 {
255  n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
256  n[1] =
257  -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
258 
259  for (unsigned int i = 0; i < 2; ++i)
260  {
261  Assert(n[i].norm() != 0,
262  ExcInternalError("Something went wrong. The "
263  "computed normals have "
264  "zero length."));
265  n[i] /= n[i].norm();
266  }
267 }
268 
269 
270 
271 template <>
272 void
275  FaceVertexNormals & n) const
276 {
277  n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
278  get_tangent_vector(face->vertex(0), face->vertex(2)));
279 
280  n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
281  get_tangent_vector(face->vertex(1), face->vertex(0)));
282 
283  n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
284  get_tangent_vector(face->vertex(2), face->vertex(3)));
285 
286  n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
287  get_tangent_vector(face->vertex(3), face->vertex(1)));
288 
289  for (unsigned int i = 0; i < 4; ++i)
290  {
291  Assert(n[i].norm() != 0,
292  ExcInternalError("Something went wrong. The "
293  "computed normals have "
294  "zero length."));
295  n[i] /= n[i].norm();
296  }
297 }
298 
299 
300 
301 template <int dim, int spacedim>
302 void
304  const typename Triangulation<dim, spacedim>::face_iterator &face,
305  FaceVertexNormals & n) const
306 {
307  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
308  {
309  n[v] = normal_vector(face, face->vertex(v));
310  n[v] /= n[v].norm();
311  }
312 }
313 
314 
315 
316 template <int dim, int spacedim>
319  const typename Triangulation<dim, spacedim>::line_iterator &line) const
320 {
321  const auto points_weights = get_default_points_and_weights(line);
322  return get_new_point(make_array_view(points_weights.first.begin(),
323  points_weights.first.end()),
324  make_array_view(points_weights.second.begin(),
325  points_weights.second.end()));
326 }
327 
328 
329 
330 template <int dim, int spacedim>
333  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
334 {
335  const auto points_weights = get_default_points_and_weights(quad);
336  return get_new_point(make_array_view(points_weights.first.begin(),
337  points_weights.first.end()),
338  make_array_view(points_weights.second.begin(),
339  points_weights.second.end()));
340 }
341 
342 
343 
344 template <int dim, int spacedim>
347  const typename Triangulation<dim, spacedim>::face_iterator &face) const
348 {
349  Assert(dim > 1, ExcImpossibleInDim(dim));
350 
351  switch (dim)
352  {
353  case 2:
354  return get_new_point_on_line(face);
355  case 3:
356  return get_new_point_on_quad(face);
357  }
358 
359  return {};
360 }
361 
362 
363 
364 template <int dim, int spacedim>
367  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
368 {
369  switch (dim)
370  {
371  case 1:
372  return get_new_point_on_line(cell);
373  case 2:
374  return get_new_point_on_quad(cell);
375  case 3:
376  return get_new_point_on_hex(cell);
377  }
378 
379  return {};
380 }
381 
382 
383 
384 template <>
385 Point<1>
388 {
389  Assert(false, ExcImpossibleInDim(1));
390  return {};
391 }
392 
393 
394 
395 template <>
396 Point<2>
399 {
400  Assert(false, ExcImpossibleInDim(1));
401  return {};
402 }
403 
404 
405 
406 template <>
407 Point<3>
410 {
411  Assert(false, ExcImpossibleInDim(1));
412  return {};
413 }
414 
415 
416 
417 template <>
418 Point<1>
421 {
422  Assert(false, ExcImpossibleInDim(1));
423  return {};
424 }
425 
426 
427 
428 template <>
429 Point<2>
432 {
433  Assert(false, ExcImpossibleInDim(1));
434  return {};
435 }
436 
437 
438 
439 template <>
440 Point<3>
443 {
444  Assert(false, ExcImpossibleInDim(1));
445  return {};
446 }
447 
448 
449 
450 template <int dim, int spacedim>
453  const typename Triangulation<dim, spacedim>::hex_iterator & /*hex*/) const
454 {
455  Assert(false, ExcImpossibleInDim(dim));
456  return {};
457 }
458 
459 
460 
461 template <>
462 Point<3>
464  const Triangulation<3, 3>::hex_iterator &hex) const
465 {
466  const auto points_weights = get_default_points_and_weights(hex, true);
467  return get_new_point(make_array_view(points_weights.first.begin(),
468  points_weights.first.end()),
469  make_array_view(points_weights.second.begin(),
470  points_weights.second.end()));
471 }
472 
473 
474 
475 template <int dim, int spacedim>
478  const Point<spacedim> &x2) const
479 {
480  const double epsilon = 1e-8;
481 
482  const std::array<Point<spacedim>, 2> points{{x1, x2}};
483  const std::array<double, 2> weights{{epsilon, 1.0 - epsilon}};
484  const Point<spacedim> neighbor_point =
485  get_new_point(make_array_view(points.begin(), points.end()),
486  make_array_view(weights.begin(), weights.end()));
487  return (neighbor_point - x1) / epsilon;
488 }
489 
490 /* -------------------------- FlatManifold --------------------- */
491 
492 namespace internal
493 {
494  namespace
495  {
497  normalized_alternating_product(const Tensor<1, 3> (&)[1])
498  {
499  // we get here from FlatManifold<2,3>::normal_vector, but
500  // the implementation below is bogus for this case anyway
501  // (see the assert at the beginning of that function).
502  Assert(false, ExcNotImplemented());
503  return {};
504  }
505 
506 
507 
509  normalized_alternating_product(const Tensor<1, 3> (&basis_vectors)[2])
510  {
511  Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
512  return tmp / tmp.norm();
513  }
514 
515  } // namespace
516 } // namespace internal
517 
518 template <int dim, int spacedim>
520  const Tensor<1, spacedim> &periodicity,
521  const double tolerance)
522  : periodicity(periodicity)
523  , tolerance(tolerance)
524 {}
525 
526 
527 
528 template <int dim, int spacedim>
529 std::unique_ptr<Manifold<dim, spacedim>>
531 {
532  return std::make_unique<FlatManifold<dim, spacedim>>(periodicity, tolerance);
533 }
534 
535 
536 
537 template <int dim, int spacedim>
540  const ArrayView<const Point<spacedim>> &surrounding_points,
541  const ArrayView<const double> & weights) const
542 {
543  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
544  1e-10,
545  ExcMessage("The weights for the individual points should sum to 1!"));
546 
547  Point<spacedim> p;
548 
549  // if there is no periodicity, use a shortcut
551  {
552  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
553  p += surrounding_points[i] * weights[i];
554  }
555  else
556  {
558 
559  for (unsigned int d = 0; d < spacedim; ++d)
560  if (periodicity[d] > 0)
561  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
562  {
563  minP[d] = std::min(minP[d], surrounding_points[i][d]);
564  Assert((surrounding_points[i][d] <
565  periodicity[d] + tolerance * periodicity[d]) ||
566  (surrounding_points[i][d] >=
567  -tolerance * periodicity[d]),
568  ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
569  }
570 
571  // compute the weighted average point, possibly taking into account
572  // periodicity
573  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
574  {
575  Point<spacedim> dp;
576  for (unsigned int d = 0; d < spacedim; ++d)
577  if (periodicity[d] > 0)
578  dp[d] =
579  ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0 ?
580  -periodicity[d] :
581  0.0);
582 
583  p += (surrounding_points[i] + dp) * weights[i];
584  }
585 
586  // if necessary, also adjust the weighted point by the periodicity
587  for (unsigned int d = 0; d < spacedim; ++d)
588  if (periodicity[d] > 0)
589  if (p[d] < 0)
590  p[d] += periodicity[d];
591  }
592 
593  return project_to_manifold(surrounding_points, p);
594 }
595 
596 
597 
598 template <int dim, int spacedim>
599 void
601  const ArrayView<const Point<spacedim>> &surrounding_points,
602  const Table<2, double> & weights,
603  ArrayView<Point<spacedim>> new_points) const
604 {
605  AssertDimension(surrounding_points.size(), weights.size(1));
606  if (weights.size(0) == 0)
607  return;
608 
609  const std::size_t n_points = surrounding_points.size();
610 
612  for (unsigned int d = 0; d < spacedim; ++d)
613  if (periodicity[d] > 0)
614  for (unsigned int i = 0; i < n_points; ++i)
615  {
616  minP[d] = std::min(minP[d], surrounding_points[i][d]);
617  Assert((surrounding_points[i][d] <
618  periodicity[d] + tolerance * periodicity[d]) ||
619  (surrounding_points[i][d] >= -tolerance * periodicity[d]),
620  ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
621  }
622 
623  // check whether periodicity shifts some of the points. Only do this if
624  // necessary to avoid memory allocation
625  const Point<spacedim> *surrounding_points_start = surrounding_points.data();
626 
627  boost::container::small_vector<Point<spacedim>, 200> modified_points;
628  bool adjust_periodicity = false;
629  for (unsigned int d = 0; d < spacedim; ++d)
630  if (periodicity[d] > 0)
631  for (unsigned int i = 0; i < n_points; ++i)
632  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
633  {
634  adjust_periodicity = true;
635  break;
636  }
637  if (adjust_periodicity == true)
638  {
639  modified_points.resize(surrounding_points.size());
640  std::copy(surrounding_points.begin(),
641  surrounding_points.end(),
642  modified_points.begin());
643  for (unsigned int d = 0; d < spacedim; ++d)
644  if (periodicity[d] > 0)
645  for (unsigned int i = 0; i < n_points; ++i)
646  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
647  modified_points[i][d] -= periodicity[d];
648  surrounding_points_start = modified_points.data();
649  }
650 
651  // Now perform the interpolation
652  for (unsigned int row = 0; row < weights.size(0); ++row)
653  {
654  Assert(
655  std::abs(
656  std::accumulate(&weights(row, 0), &weights(row, 0) + n_points, 0.0) -
657  1.0) < 1e-10,
658  ExcMessage("The weights for the individual points should sum to 1!"));
659  Point<spacedim> new_point;
660  for (unsigned int p = 0; p < n_points; ++p)
661  new_point += surrounding_points_start[p] * weights(row, p);
662 
663  // if necessary, also adjust the weighted point by the periodicity
664  for (unsigned int d = 0; d < spacedim; ++d)
665  if (periodicity[d] > 0)
666  if (new_point[d] < 0)
667  new_point[d] += periodicity[d];
668 
669  // TODO should this use surrounding_points_start or surrounding_points?
670  // The older version used surrounding_points
671  new_points[row] =
672  project_to_manifold(make_array_view(surrounding_points.begin(),
673  surrounding_points.end()),
674  new_point);
675  }
676 }
677 
678 
679 
680 template <int dim, int spacedim>
683  const ArrayView<const Point<spacedim>> & /*vertices*/,
684  const Point<spacedim> &candidate) const
685 {
686  return candidate;
687 }
688 
689 
690 
691 template <int dim, int spacedim>
692 const Tensor<1, spacedim> &
694 {
695  return periodicity;
696 }
697 
698 
699 
700 template <int dim, int spacedim>
703  const Point<spacedim> &x2) const
704 {
705  Tensor<1, spacedim> direction = x2 - x1;
706 
707  // see if we have to take into account periodicity. if so, we need
708  // to make sure that if a distance in one coordinate direction
709  // is larger than half of the box length, then go the other way
710  // around (i.e., via the periodic box)
711  for (unsigned int d = 0; d < spacedim; ++d)
712  if (periodicity[d] > tolerance)
713  {
714  if (direction[d] < -periodicity[d] / 2)
715  direction[d] += periodicity[d];
716  else if (direction[d] > periodicity[d] / 2)
717  direction[d] -= periodicity[d];
718  }
719 
720  return direction;
721 }
722 
723 
724 
725 template <>
726 void
730 {
731  Assert(false, ExcImpossibleInDim(1));
732 }
733 
734 
735 
736 template <>
737 void
741 {
742  Assert(false, ExcNotImplemented());
743 }
744 
745 
746 
747 template <>
748 void
752 {
753  Assert(false, ExcNotImplemented());
754 }
755 
756 
757 
758 template <>
759 void
762  Manifold<2, 2>::FaceVertexNormals & face_vertex_normals) const
763 {
764  const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
765  for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_face;
766  ++vertex)
767  // compute normals from tangent
768  face_vertex_normals[vertex] = Point<2>(tangent[1], -tangent[0]);
769 }
770 
771 
772 
773 template <>
774 void
776  const Triangulation<2, 3>::face_iterator & /*face*/,
777  Manifold<2, 3>::FaceVertexNormals & /*face_vertex_normals*/) const
778 {
779  Assert(false, ExcNotImplemented());
780 }
781 
782 
783 
784 template <>
785 void
788  Manifold<3, 3>::FaceVertexNormals & face_vertex_normals) const
789 {
790  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
791 
792  static const unsigned int neighboring_vertices[4][2] = {{1, 2},
793  {3, 0},
794  {0, 3},
795  {2, 1}};
796  for (unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
797  {
798  // first define the two tangent vectors at the vertex by using the
799  // two lines radiating away from this vertex
800  const Tensor<1, 3> tangents[2] = {
801  face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
802  face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
803 
804  // then compute the normal by taking the cross product. since the
805  // normal is not required to be normalized, no problem here
806  face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
807  }
808 }
809 
810 
811 
812 template <>
815  const Point<1> &) const
816 {
817  Assert(false, ExcNotImplemented());
818  return {};
819 }
820 
821 
822 
823 template <>
826  const Point<2> &) const
827 {
828  Assert(false, ExcNotImplemented());
829  return {};
830 }
831 
832 
833 
834 template <>
837  const Point<3> &) const
838 {
839  Assert(false, ExcNotImplemented());
840  return {};
841 }
842 
843 
844 
845 template <>
849  const Point<2> & p) const
850 {
851  // In 2d, a face is just a straight line and
852  // we can use the 'standard' implementation.
853  return Manifold<2, 2>::normal_vector(face, p);
854 }
855 
856 
857 
858 template <int dim, int spacedim>
861  const typename Triangulation<dim, spacedim>::face_iterator &face,
862  const Point<spacedim> & p) const
863 {
864  // I don't think the implementation below will work when dim!=spacedim;
865  // in fact, I believe that we don't even have enough information here,
866  // because we would need to know not only about the tangent vectors
867  // of the face, but also of the cell, to compute the normal vector.
868  // Someone will have to think about this some more.
869  Assert(dim == spacedim, ExcNotImplemented());
870 
871  // in order to find out what the normal vector is, we first need to
872  // find the reference coordinates of the point p on the given face,
873  // or at least the reference coordinates of the closest point on the
874  // face
875  //
876  // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
877  // where F(xi) is the mapping. this algorithm is implemented in
878  // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
879  // while we need it for faces here. it's also implemented in somewhat
880  // more generality there using the machinery of the MappingQ1 class
881  // while we really only need it for a specific case here
882  //
883  // in any case, the iteration we use here is a Gauss-Newton's iteration with
884  // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
885  // where
886  // J(xi) = (grad F(xi))^T (F(xi)-p)
887  // and
888  // H(xi) = [grad F(xi)]^T [grad F(xi)]
889  // In all this,
890  // F(xi) = sum_v vertex[v] phi_v(xi)
891  // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
892 
893  // we start with the point xi=1/2, xi=(1/2,1/2), ...
894  const unsigned int facedim = dim - 1;
895 
896  Point<facedim> xi;
897 
898  const auto face_reference_cell = face->reference_cell();
899 
900  if (face_reference_cell == ReferenceCells::get_hypercube<facedim>())
901  {
902  for (unsigned int i = 0; i < facedim; ++i)
903  xi[i] = 1. / 2;
904  }
905  else
906  {
907  for (unsigned int i = 0; i < facedim; ++i)
908  xi[i] = 1. / 3;
909  }
910 
911  const double eps = 1e-12;
912  Tensor<1, spacedim> grad_F[facedim];
913  unsigned int iteration = 0;
914  while (true)
915  {
917  for (const unsigned int v : face->vertex_indices())
918  F +=
919  face->vertex(v) * face_reference_cell.d_linear_shape_function(xi, v);
920 
921  for (unsigned int i = 0; i < facedim; ++i)
922  {
923  grad_F[i] = 0;
924  for (const unsigned int v : face->vertex_indices())
925  grad_F[i] +=
926  face->vertex(v) *
927  face_reference_cell.d_linear_shape_function_gradient(xi, v)[i];
928  }
929 
931  for (unsigned int i = 0; i < facedim; ++i)
932  for (unsigned int j = 0; j < spacedim; ++j)
933  J[i] += grad_F[i][j] * (F - p)[j];
934 
936  for (unsigned int i = 0; i < facedim; ++i)
937  for (unsigned int j = 0; j < facedim; ++j)
938  for (unsigned int k = 0; k < spacedim; ++k)
939  H[i][j] += grad_F[i][k] * grad_F[j][k];
940 
941  const Tensor<1, facedim> delta_xi = -invert(H) * J;
942  xi += delta_xi;
943  ++iteration;
944 
945  Assert(iteration < 10,
946  ExcMessage("The Newton iteration to find the reference point "
947  "did not converge in 10 iterations. Do you have a "
948  "deformed cell? (See the glossary for a definition "
949  "of what a deformed cell is. You may want to output "
950  "the vertices of your cell."));
951 
952  // It turns out that the check in reference coordinates with an absolute
953  // tolerance can cause a convergence failure of the Newton method as
954  // seen in tests/manifold/flat_manifold_09.cc. To work around this, also
955  // use a convergence check in world coordinates. This check has to be
956  // relative to the size of the face of course. Here we decided to use
957  // diameter because it works for non-planar faces and is cheap to
958  // compute:
959  const double normalized_delta_world = (F - p).norm() / face->diameter();
960 
961  if (delta_xi.norm() < eps || normalized_delta_world < eps)
962  break;
963  }
964 
965  // so now we have the reference coordinates xi of the point p.
966  // we then have to compute the normal vector, which we can do
967  // by taking the (normalize) alternating product of all the tangent
968  // vectors given by grad_F
969  return internal::normalized_alternating_product(grad_F);
970 }
971 
972 
973 /* -------------------------- ChartManifold --------------------- */
974 template <int dim, int spacedim, int chartdim>
977  : sub_manifold(periodicity)
978 {}
979 
980 
981 
982 template <int dim, int spacedim, int chartdim>
985  const Point<spacedim> &p1,
986  const Point<spacedim> &p2,
987  const double w) const
988 {
989  const std::array<Point<spacedim>, 2> points{{p1, p2}};
990  const std::array<double, 2> weights{{1. - w, w}};
991  return get_new_point(make_array_view(points.begin(), points.end()),
992  make_array_view(weights.begin(), weights.end()));
993 }
994 
995 
996 
997 template <int dim, int spacedim, int chartdim>
1000  const ArrayView<const Point<spacedim>> &surrounding_points,
1001  const ArrayView<const double> & weights) const
1002 {
1003  const std::size_t n_points = surrounding_points.size();
1004 
1005  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1006 
1007  for (unsigned int i = 0; i < n_points; ++i)
1008  chart_points[i] = pull_back(surrounding_points[i]);
1009 
1010  const Point<chartdim> p_chart = sub_manifold.get_new_point(
1011  make_array_view(chart_points.begin(), chart_points.end()), weights);
1012 
1013  return push_forward(p_chart);
1014 }
1015 
1016 
1017 
1018 template <int dim, int spacedim, int chartdim>
1019 void
1021  const ArrayView<const Point<spacedim>> &surrounding_points,
1022  const Table<2, double> & weights,
1023  ArrayView<Point<spacedim>> new_points) const
1024 {
1025  Assert(weights.size(0) > 0, ExcEmptyObject());
1026  AssertDimension(surrounding_points.size(), weights.size(1));
1027 
1028  const std::size_t n_points = surrounding_points.size();
1029 
1030  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1031  for (std::size_t i = 0; i < n_points; ++i)
1032  chart_points[i] = pull_back(surrounding_points[i]);
1033 
1034  boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1035  weights.size(0));
1037  make_array_view(chart_points.begin(), chart_points.end()),
1038  weights,
1039  make_array_view(new_points_on_chart.begin(), new_points_on_chart.end()));
1040 
1041  for (std::size_t row = 0; row < weights.size(0); ++row)
1042  new_points[row] = push_forward(new_points_on_chart[row]);
1043 }
1044 
1045 
1046 
1047 template <int dim, int spacedim, int chartdim>
1050  const Point<chartdim> &) const
1051 {
1052  // function must be implemented in a derived class to be usable,
1053  // as discussed in this function's documentation
1054  Assert(false, ExcPureFunctionCalled());
1056 }
1057 
1058 
1059 
1060 template <int dim, int spacedim, int chartdim>
1063  const Point<spacedim> &x1,
1064  const Point<spacedim> &x2) const
1065 {
1068 
1069  // ensure that the chart is not singular by asserting that its
1070  // derivative has a positive determinant. we need to make this
1071  // comparison relative to the size of the derivative. since the
1072  // determinant is the product of chartdim factors, take the
1073  // chartdim-th root of it in comparing against the size of the
1074  // derivative
1075  Assert(std::pow(std::abs(F_prime.determinant()), 1. / chartdim) >=
1076  1e-12 * F_prime.norm(),
1077  ExcMessage(
1078  "The derivative of a chart function must not be singular."));
1079 
1080  const Tensor<1, chartdim> delta =
1082 
1083  Tensor<1, spacedim> result;
1084  for (unsigned int i = 0; i < spacedim; ++i)
1085  result[i] += F_prime[i] * delta;
1086 
1087  return result;
1088 }
1089 
1090 
1091 
1092 template <int dim, int spacedim, int chartdim>
1093 const Tensor<1, chartdim> &
1095 {
1096  return sub_manifold.get_periodicity();
1097 }
1098 
1099 // explicit instantiations
1100 #include "manifold.inst"
1101 
static ::ExceptionBase & ExcPureFunctionCalled()
static const unsigned int invalid_unsigned_int
Definition: types.h:196
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
Definition: manifold.cc:519
std::pair< std::array< Point< MeshIteratorType::AccessorType::space_dimension >, n_default_points_per_cell< MeshIteratorType >)>, std::array< double, n_default_points_per_cell< MeshIteratorType >)> > get_default_points_and_weights(const MeshIteratorType &iterator, const bool with_interpolation=false)
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
Definition: manifold.cc:860
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1448
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const Tensor< 1, spacedim > & get_periodicity() const
Definition: manifold.cc:693
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:452
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1496
iterator end() const
Definition: array_view.h:588
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
Definition: manifold.cc:65
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Number determinant() const
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2443
std::size_t size() const
Definition: array_view.h:570
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:407
static ::ExceptionBase & ExcMessage(std::string arg1)
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:318
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:477
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1049
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
Definition: manifold.cc:40
#define Assert(cond, exc)
Definition: exceptions.h:1466
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
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:693
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1094
Point< 3 > vertices[4]
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
Definition: manifold.cc:682
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold.cc:530
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1472
const double tolerance
Definition: manifold.h:808
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:52
const FlatManifold< chartdim, chartdim > sub_manifold
Definition: manifold.h:1085
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
Definition: manifold.cc:1020
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1062
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const
Definition: manifold.cc:125
size_type size(const unsigned int i) const
Definition: tensor.h:449
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:702
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:999
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:393
T min(const T &t, const MPI_Comm &mpi_communicator)
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:332
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:539
static ::ExceptionBase & ExcEmptyObject()
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, typename Manifold< dim, spacedim >::FaceVertexNormals &face_vertex_normals) const override
virtual Point< spacedim > push_forward(const Point< chartdim > &chart_point) const =0
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:307
static ::ExceptionBase & ExcNotImplemented()
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:239
iterator begin() const
Definition: array_view.h:579
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
numbers::NumberTraits< Number >::real_type norm() const
virtual void get_new_points(const ArrayView< const Point< spacedim >> &surrounding_points, const Table< 2, double > &weights, ArrayView< Point< spacedim >> new_points) const override
Definition: manifold.cc:600
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2418
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:975
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:346
void copy(const T *begin, const T *end, U *dest)
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:303
virtual Point< chartdim > pull_back(const Point< spacedim > &space_point) const =0
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:984
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:794
static ::ExceptionBase & ExcInternalError()
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:366