Reference documentation for deal.II version GIT 196b8d1939 2022-07-03 13:35:01+00:00
\(\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 - 2022 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>
23 #include <deal.II/grid/tria.h>
26 
28 #include <boost/container/small_vector.hpp>
30 
31 #include <cmath>
32 #include <memory>
33 
35 
36 /* -------------------------- Manifold --------------------- */
37 template <int dim, int spacedim>
40  const ArrayView<const Point<spacedim>> &,
41  const Point<spacedim> &) const
42 {
43  Assert(false, ExcPureFunctionCalled());
44  return {};
45 }
46 
47 
48 
49 template <int dim, int spacedim>
52  const Point<spacedim> &p2,
53  const double w) const
54 {
55  const std::array<Point<spacedim>, 2> vertices{{p1, p2}};
56  return project_to_manifold(make_array_view(vertices.begin(), vertices.end()),
57  w * p2 + (1 - w) * p1);
58 }
59 
60 
61 
62 template <int dim, int spacedim>
65  const ArrayView<const Point<spacedim>> &surrounding_points,
66  const ArrayView<const double> & weights) const
67 {
68  const double tol = 1e-10;
69  const unsigned int n_points = surrounding_points.size();
70 
71  Assert(n_points > 0, ExcMessage("There should be at least one point."));
72 
73  Assert(n_points == weights.size(),
74  ExcMessage(
75  "There should be as many surrounding points as weights given."));
76 
77  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
78  tol,
79  ExcMessage("The weights for the individual points should sum to 1!"));
80 
81  // First sort points in the order of their weights. This is done to
82  // produce unique points even if get_intermediate_points is not
83  // associative (as for the SphericalManifold).
84  boost::container::small_vector<unsigned int, 100> permutation(n_points);
85  std::iota(permutation.begin(), permutation.end(), 0u);
86  std::sort(permutation.begin(),
87  permutation.end(),
88  [&weights](const std::size_t a, const std::size_t b) {
89  return weights[a] < weights[b];
90  });
91 
92  // Now loop over points in the order of their associated weight
93  Point<spacedim> p = surrounding_points[permutation[0]];
94  double w = weights[permutation[0]];
95 
96  for (unsigned int i = 1; i < n_points; ++i)
97  {
98  double weight = 0.0;
99  if (std::abs(weights[permutation[i]] + w) < tol)
100  weight = 0.0;
101  else
102  weight = w / (weights[permutation[i]] + w);
103 
104  if (std::abs(weight) > 1e-14)
105  {
106  p = get_intermediate_point(p,
107  surrounding_points[permutation[i]],
108  1.0 - weight);
109  }
110  else
111  {
112  p = surrounding_points[permutation[i]];
113  }
114  w += weights[permutation[i]];
115  }
116 
117  return p;
118 }
119 
120 
121 
122 template <int dim, int spacedim>
123 void
125  const ArrayView<const Point<spacedim>> &surrounding_points,
126  const Table<2, double> & weights,
127  ArrayView<Point<spacedim>> new_points) const
128 {
129  AssertDimension(surrounding_points.size(), weights.size(1));
130 
131  for (unsigned int row = 0; row < weights.size(0); ++row)
132  {
133  new_points[row] =
134  get_new_point(make_array_view(surrounding_points.begin(),
135  surrounding_points.end()),
136  make_array_view(weights, row));
137  }
138 }
139 
140 
141 
142 template <>
145  const Point<2> & p) const
146 {
147  const int spacedim = 2;
148 
149  // get the tangent vector from the point 'p' in the direction of the further
150  // one of the two vertices that make up the face of this 2d cell
151  const Tensor<1, spacedim> tangent =
152  ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
153  -get_tangent_vector(p, face->vertex(0)) :
154  get_tangent_vector(p, face->vertex(1)));
155 
156  // then rotate it by 90 degrees
157  const Tensor<1, spacedim> normal = cross_product_2d(tangent);
158  return normal / normal.norm();
159 }
160 
161 
162 
163 template <>
166  const Point<3> & p) const
167 {
168  const int spacedim = 3;
169 
170  const std::array<Point<spacedim>, 4> vertices{
171  {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
172  const std::array<double, 4> distances{{vertices[0].distance(p),
173  vertices[1].distance(p),
174  vertices[2].distance(p),
175  vertices[3].distance(p)}};
176  const double max_distance = std::max(std::max(distances[0], distances[1]),
177  std::max(distances[2], distances[3]));
178 
179  // We need to find two tangential vectors to the given point p, but we do
180  // not know how the point is oriented against the face. We guess the two
181  // directions by assuming a flat topology and take the two directions that
182  // indicate the angle closest to a perpendicular one (i.e., cos(theta) close
183  // to zero). We start with an invalid value but the loops should always find
184  // a value.
185  double abs_cos_angle = std::numeric_limits<double>::max();
186  unsigned int first_index = numbers::invalid_unsigned_int,
187  second_index = numbers::invalid_unsigned_int;
188  for (unsigned int i = 0; i < 3; ++i)
189  if (distances[i] > 1e-8 * max_distance)
190  for (unsigned int j = i + 1; j < 4; ++j)
191  if (distances[j] > 1e-8 * max_distance)
192  {
193  const double new_angle = (p - vertices[i]) * (p - vertices[j]) /
194  (distances[i] * distances[j]);
195  // multiply by factor 0.999 to bias the search in a way that
196  // avoids trouble with roundoff
197  if (std::abs(new_angle) < 0.999 * abs_cos_angle)
198  {
199  abs_cos_angle = std::abs(new_angle);
200  first_index = i;
201  second_index = j;
202  }
203  }
204  Assert(first_index != numbers::invalid_unsigned_int,
205  ExcMessage("The search for possible directions did not succeed."));
206 
207  // Compute tangents and normal for selected vertices
208  Tensor<1, spacedim> t1 = get_tangent_vector(p, vertices[first_index]);
209  Tensor<1, spacedim> t2 = get_tangent_vector(p, vertices[second_index]);
210  Tensor<1, spacedim> normal = cross_product_3d(t1, t2);
211 
212  Assert(
214  t1.norm_square() * t2.norm_square(),
215  ExcMessage(
216  "Manifold::normal_vector was unable to find a suitable combination "
217  "of vertices to compute a normal on this face. We chose the secants "
218  "that are as orthogonal as possible, but tangents appear to be "
219  "linearly dependent. Check for distorted faces in your triangulation."));
220 
221  // Now figure out if we need to flip the direction, we do this by comparing
222  // to a reference normal that would be the correct result if all vertices
223  // would lie in a plane
224  const Tensor<1, spacedim> rt1 = vertices[3] - vertices[0];
225  const Tensor<1, spacedim> rt2 = vertices[2] - vertices[1];
226  const Tensor<1, spacedim> reference_normal = cross_product_3d(rt1, rt2);
227 
228  if (reference_normal * normal < 0.0)
229  normal *= -1.0;
230 
231  return normal / normal.norm();
232 }
233 
234 
235 
236 template <int dim, int spacedim>
239  const typename Triangulation<dim, spacedim>::face_iterator & /*face*/,
240  const Point<spacedim> & /*p*/) const
241 {
242  Assert(false, ExcPureFunctionCalled());
243  return Tensor<1, spacedim>();
244 }
245 
246 
247 
248 template <>
249 void
252  FaceVertexNormals & n) const
253 {
254  n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
255  n[1] =
256  -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
257 
258  for (unsigned int i = 0; i < 2; ++i)
259  {
260  Assert(n[i].norm() != 0,
261  ExcInternalError("Something went wrong. The "
262  "computed normals have "
263  "zero length."));
264  n[i] /= n[i].norm();
265  }
266 }
267 
268 
269 
270 template <>
271 void
274  FaceVertexNormals & n) const
275 {
276  n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
277  get_tangent_vector(face->vertex(0), face->vertex(2)));
278 
279  n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
280  get_tangent_vector(face->vertex(1), face->vertex(0)));
281 
282  n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
283  get_tangent_vector(face->vertex(2), face->vertex(3)));
284 
285  n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
286  get_tangent_vector(face->vertex(3), face->vertex(1)));
287 
288  for (unsigned int i = 0; i < 4; ++i)
289  {
290  Assert(n[i].norm() != 0,
291  ExcInternalError("Something went wrong. The "
292  "computed normals have "
293  "zero length."));
294  n[i] /= n[i].norm();
295  }
296 }
297 
298 
299 
300 template <int dim, int spacedim>
301 void
303  const typename Triangulation<dim, spacedim>::face_iterator &face,
304  FaceVertexNormals & n) const
305 {
306  for (unsigned int v = 0; v < face->reference_cell().n_vertices(); ++v)
307  {
308  n[v] = normal_vector(face, face->vertex(v));
309  n[v] /= n[v].norm();
310  }
311 }
312 
313 
314 
315 template <int dim, int spacedim>
318  const typename Triangulation<dim, spacedim>::line_iterator &line) const
319 {
320  const auto points_weights = Manifolds::get_default_points_and_weights(line);
321  return get_new_point(make_array_view(points_weights.first.begin(),
322  points_weights.first.end()),
323  make_array_view(points_weights.second.begin(),
324  points_weights.second.end()));
325 }
326 
327 
328 
329 template <int dim, int spacedim>
332  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
333 {
334  const auto points_weights = Manifolds::get_default_points_and_weights(quad);
335  return get_new_point(make_array_view(points_weights.first.begin(),
336  points_weights.first.end()),
337  make_array_view(points_weights.second.begin(),
338  points_weights.second.end()));
339 }
340 
341 
342 
343 template <int dim, int spacedim>
346  const typename Triangulation<dim, spacedim>::face_iterator &face) const
347 {
348  Assert(dim > 1, ExcImpossibleInDim(dim));
349 
350  switch (dim)
351  {
352  case 2:
353  return get_new_point_on_line(face);
354  case 3:
355  return get_new_point_on_quad(face);
356  }
357 
358  return {};
359 }
360 
361 
362 
363 template <int dim, int spacedim>
366  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
367 {
368  switch (dim)
369  {
370  case 1:
371  return get_new_point_on_line(cell);
372  case 2:
373  return get_new_point_on_quad(cell);
374  case 3:
375  return get_new_point_on_hex(cell);
376  }
377 
378  return {};
379 }
380 
381 
382 
383 template <>
384 Point<1>
387 {
388  Assert(false, ExcImpossibleInDim(1));
389  return {};
390 }
391 
392 
393 
394 template <>
395 Point<2>
398 {
399  Assert(false, ExcImpossibleInDim(1));
400  return {};
401 }
402 
403 
404 
405 template <>
406 Point<3>
409 {
410  Assert(false, ExcImpossibleInDim(1));
411  return {};
412 }
413 
414 
415 
416 template <>
417 Point<1>
420 {
421  Assert(false, ExcImpossibleInDim(1));
422  return {};
423 }
424 
425 
426 
427 template <>
428 Point<2>
431 {
432  Assert(false, ExcImpossibleInDim(1));
433  return {};
434 }
435 
436 
437 
438 template <>
439 Point<3>
442 {
443  Assert(false, ExcImpossibleInDim(1));
444  return {};
445 }
446 
447 
448 
449 template <int dim, int spacedim>
452  const typename Triangulation<dim, spacedim>::hex_iterator & /*hex*/) const
453 {
454  Assert(false, ExcImpossibleInDim(dim));
455  return {};
456 }
457 
458 
459 
460 template <>
461 Point<3>
463  const Triangulation<3, 3>::hex_iterator &hex) const
464 {
465  const auto points_weights =
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] <
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] <
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  // We're in 2d. Faces are edges:
766  for (const unsigned int vertex : ReferenceCells::Line.vertex_indices())
767  // compute normals from tangent
768  face_vertex_normals[vertex] = Tensor<1, 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>
976  const Tensor<1, chartdim> &periodicity)
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);
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));
1036  sub_manifold.get_new_points(
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 {
1067  push_forward_gradient(pull_back(x1));
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 =
1081  sub_manifold.get_tangent_vector(pull_back(x1), pull_back(x2));
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 
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:699
iterator begin() const
Definition: array_view.h:585
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:699
iterator end() const
Definition: array_view.h:594
std::size_t size() const
Definition: array_view.h:576
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:999
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:975
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1062
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1049
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1094
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 Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:984
Number determinant() const
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
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 > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
const Tensor< 1, spacedim > & get_periodicity() const
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const override
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:794
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
const double tolerance
Definition: manifold.h:808
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:124
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
Definition: manifold.cc:64
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
Definition: manifold.cc:39
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:302
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:331
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:307
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:451
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:238
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:51
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:345
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:477
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:365
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:317
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Definition: tensor.h:503
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2635
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:2660
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1293
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcPureFunctionCalled()
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1490
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1466
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
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)
Point< dim > push_forward(const TopoDS_Shape &in_shape, const double u, const double v)
Definition: utilities.cc:833
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 1, dim, Number > pull_back(const Tensor< 1, dim, Number > &v, const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Line
void copy(const T *begin, const T *end, U *dest)
static const unsigned int invalid_unsigned_int
Definition: types.h:201
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
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:2660