Reference documentation for deal.II version GIT 4e92747666 2022-12-02 08:40:02+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 <limits>
33 #include <memory>
34 
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
211  Tensor<1, spacedim> normal;
212 
213  bool done = false;
214  std::vector<bool> tested_vertices(vertices.size(), false);
215  tested_vertices[first_index] = true;
216  tested_vertices[second_index] = true;
217 
218  do
219  {
220  // Compute tangents and normal for selected vertices
221  t1 = get_tangent_vector(p, vertices[first_index]);
222  t2 = get_tangent_vector(p, vertices[second_index]);
223  normal = cross_product_3d(t1, t2);
224 
225  // if normal is zero, try some other combination of vertices
226  if (normal.norm_square() < 1e4 * std::numeric_limits<double>::epsilon() *
227  t1.norm_square() * t2.norm_square())
228  {
229  // See if we have tested all vertices already
230  auto first_false =
231  std::find(tested_vertices.begin(), tested_vertices.end(), false);
232  if (first_false == tested_vertices.end())
233  {
234  done = true;
235  }
236  else
237  {
238  *first_false = true;
239  second_index = first_false - tested_vertices.begin();
240  }
241  }
242  else
243  {
244  done = true;
245  }
246  }
247  while (!done);
248 
249  Assert(
251  t1.norm_square() * t2.norm_square(),
252  ExcMessage(
253  "Manifold::normal_vector was unable to find a suitable combination "
254  "of vertices to compute a normal on this face. We chose the secants "
255  "that are as orthogonal as possible, but tangents appear to be "
256  "linearly dependent. Check for distorted faces in your triangulation."));
257 
258  // Now figure out if we need to flip the direction, we do this by comparing
259  // to a reference normal that would be the correct result if all vertices
260  // would lie in a plane
261  const Tensor<1, spacedim> rt1 = vertices[3] - vertices[0];
262  const Tensor<1, spacedim> rt2 = vertices[2] - vertices[1];
263  const Tensor<1, spacedim> reference_normal = cross_product_3d(rt1, rt2);
264 
265  if (reference_normal * normal < 0.0)
266  normal *= -1.0;
267 
268  return normal / normal.norm();
269 }
270 
271 
272 
273 template <int dim, int spacedim>
276  const typename Triangulation<dim, spacedim>::face_iterator & /*face*/,
277  const Point<spacedim> & /*p*/) const
278 {
279  Assert(false, ExcPureFunctionCalled());
280  return Tensor<1, spacedim>();
281 }
282 
283 
284 
285 template <>
286 void
289  FaceVertexNormals & n) const
290 {
291  n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
292  n[1] =
293  -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
294 
295  for (unsigned int i = 0; i < 2; ++i)
296  {
297  Assert(n[i].norm() != 0,
298  ExcInternalError("Something went wrong. The "
299  "computed normals have "
300  "zero length."));
301  n[i] /= n[i].norm();
302  }
303 }
304 
305 
306 
307 template <>
308 void
311  FaceVertexNormals & n) const
312 {
313  n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
314  get_tangent_vector(face->vertex(0), face->vertex(2)));
315 
316  n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
317  get_tangent_vector(face->vertex(1), face->vertex(0)));
318 
319  n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
320  get_tangent_vector(face->vertex(2), face->vertex(3)));
321 
322  n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
323  get_tangent_vector(face->vertex(3), face->vertex(1)));
324 
325  for (unsigned int i = 0; i < 4; ++i)
326  {
327  Assert(n[i].norm() != 0,
328  ExcInternalError("Something went wrong. The "
329  "computed normals have "
330  "zero length."));
331  n[i] /= n[i].norm();
332  }
333 }
334 
335 
336 
337 template <int dim, int spacedim>
338 void
340  const typename Triangulation<dim, spacedim>::face_iterator &face,
341  FaceVertexNormals & n) const
342 {
343  for (unsigned int v = 0; v < face->reference_cell().n_vertices(); ++v)
344  {
345  n[v] = normal_vector(face, face->vertex(v));
346  n[v] /= n[v].norm();
347  }
348 }
349 
350 
351 
352 template <int dim, int spacedim>
355  const typename Triangulation<dim, spacedim>::line_iterator &line) const
356 {
357  const auto points_weights = Manifolds::get_default_points_and_weights(line);
358  return get_new_point(make_array_view(points_weights.first.begin(),
359  points_weights.first.end()),
360  make_array_view(points_weights.second.begin(),
361  points_weights.second.end()));
362 }
363 
364 
365 
366 template <int dim, int spacedim>
369  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
370 {
371  const auto points_weights = Manifolds::get_default_points_and_weights(quad);
372  return get_new_point(make_array_view(points_weights.first.begin(),
373  points_weights.first.end()),
374  make_array_view(points_weights.second.begin(),
375  points_weights.second.end()));
376 }
377 
378 
379 
380 template <int dim, int spacedim>
383  const typename Triangulation<dim, spacedim>::face_iterator &face) const
384 {
385  Assert(dim > 1, ExcImpossibleInDim(dim));
386 
387  switch (dim)
388  {
389  case 2:
390  return get_new_point_on_line(face);
391  case 3:
392  return get_new_point_on_quad(face);
393  }
394 
395  return {};
396 }
397 
398 
399 
400 template <int dim, int spacedim>
403  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
404 {
405  switch (dim)
406  {
407  case 1:
408  return get_new_point_on_line(cell);
409  case 2:
410  return get_new_point_on_quad(cell);
411  case 3:
412  return get_new_point_on_hex(cell);
413  }
414 
415  return {};
416 }
417 
418 
419 
420 template <>
421 Point<1>
424 {
425  Assert(false, ExcImpossibleInDim(1));
426  return {};
427 }
428 
429 
430 
431 template <>
432 Point<2>
435 {
436  Assert(false, ExcImpossibleInDim(1));
437  return {};
438 }
439 
440 
441 
442 template <>
443 Point<3>
446 {
447  Assert(false, ExcImpossibleInDim(1));
448  return {};
449 }
450 
451 
452 
453 template <>
454 Point<1>
457 {
458  Assert(false, ExcImpossibleInDim(1));
459  return {};
460 }
461 
462 
463 
464 template <>
465 Point<2>
468 {
469  Assert(false, ExcImpossibleInDim(1));
470  return {};
471 }
472 
473 
474 
475 template <>
476 Point<3>
479 {
480  Assert(false, ExcImpossibleInDim(1));
481  return {};
482 }
483 
484 
485 
486 template <int dim, int spacedim>
489  const typename Triangulation<dim, spacedim>::hex_iterator & /*hex*/) const
490 {
491  Assert(false, ExcImpossibleInDim(dim));
492  return {};
493 }
494 
495 
496 
497 template <>
498 Point<3>
500  const Triangulation<3, 3>::hex_iterator &hex) const
501 {
502  const auto points_weights =
504  return get_new_point(make_array_view(points_weights.first.begin(),
505  points_weights.first.end()),
506  make_array_view(points_weights.second.begin(),
507  points_weights.second.end()));
508 }
509 
510 
511 
512 template <int dim, int spacedim>
515  const Point<spacedim> &x2) const
516 {
517  const double epsilon = 1e-8;
518 
519  const std::array<Point<spacedim>, 2> points{{x1, x2}};
520  const std::array<double, 2> weights{{epsilon, 1.0 - epsilon}};
521  const Point<spacedim> neighbor_point =
522  get_new_point(make_array_view(points.begin(), points.end()),
523  make_array_view(weights.begin(), weights.end()));
524  return (neighbor_point - x1) / epsilon;
525 }
526 
527 /* -------------------------- FlatManifold --------------------- */
528 
529 namespace internal
530 {
531  namespace
532  {
534  normalized_alternating_product(const Tensor<1, 3> (&)[1])
535  {
536  // we get here from FlatManifold<2,3>::normal_vector, but
537  // the implementation below is bogus for this case anyway
538  // (see the assert at the beginning of that function).
540  return {};
541  }
542 
543 
544 
546  normalized_alternating_product(const Tensor<1, 3> (&basis_vectors)[2])
547  {
548  Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
549  return tmp / tmp.norm();
550  }
551 
552  } // namespace
553 } // namespace internal
554 
555 template <int dim, int spacedim>
557  const Tensor<1, spacedim> &periodicity,
558  const double tolerance)
559  : periodicity(periodicity)
560  , tolerance(tolerance)
561 {}
562 
563 
564 
565 template <int dim, int spacedim>
566 std::unique_ptr<Manifold<dim, spacedim>>
568 {
569  return std::make_unique<FlatManifold<dim, spacedim>>(periodicity, tolerance);
570 }
571 
572 
573 
574 template <int dim, int spacedim>
577  const ArrayView<const Point<spacedim>> &surrounding_points,
578  const ArrayView<const double> & weights) const
579 {
580  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
581  1e-10,
582  ExcMessage("The weights for the individual points should sum to 1!"));
583 
584  Point<spacedim> p;
585 
586  // if there is no periodicity, use a shortcut
588  {
589  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
590  p += surrounding_points[i] * weights[i];
591  }
592  else
593  {
595 
596  for (unsigned int d = 0; d < spacedim; ++d)
597  if (periodicity[d] > 0)
598  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
599  {
600  minP[d] = std::min(minP[d], surrounding_points[i][d]);
601  Assert((surrounding_points[i][d] <
603  (surrounding_points[i][d] >=
604  -tolerance * periodicity[d]),
605  ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
606  }
607 
608  // compute the weighted average point, possibly taking into account
609  // periodicity
610  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
611  {
612  Point<spacedim> dp;
613  for (unsigned int d = 0; d < spacedim; ++d)
614  if (periodicity[d] > 0)
615  dp[d] =
616  ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0 ?
618  0.0);
619 
620  p += (surrounding_points[i] + dp) * weights[i];
621  }
622 
623  // if necessary, also adjust the weighted point by the periodicity
624  for (unsigned int d = 0; d < spacedim; ++d)
625  if (periodicity[d] > 0)
626  if (p[d] < 0)
627  p[d] += periodicity[d];
628  }
629 
630  return project_to_manifold(surrounding_points, p);
631 }
632 
633 
634 
635 template <int dim, int spacedim>
636 void
638  const ArrayView<const Point<spacedim>> &surrounding_points,
639  const Table<2, double> & weights,
640  ArrayView<Point<spacedim>> new_points) const
641 {
642  AssertDimension(surrounding_points.size(), weights.size(1));
643  if (weights.size(0) == 0)
644  return;
645 
646  const std::size_t n_points = surrounding_points.size();
647 
649  for (unsigned int d = 0; d < spacedim; ++d)
650  if (periodicity[d] > 0)
651  for (unsigned int i = 0; i < n_points; ++i)
652  {
653  minP[d] = std::min(minP[d], surrounding_points[i][d]);
654  Assert((surrounding_points[i][d] <
656  (surrounding_points[i][d] >= -tolerance * periodicity[d]),
657  ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
658  }
659 
660  // check whether periodicity shifts some of the points. Only do this if
661  // necessary to avoid memory allocation
662  const Point<spacedim> *surrounding_points_start = surrounding_points.data();
663 
664  boost::container::small_vector<Point<spacedim>, 200> modified_points;
665  bool adjust_periodicity = false;
666  for (unsigned int d = 0; d < spacedim; ++d)
667  if (periodicity[d] > 0)
668  for (unsigned int i = 0; i < n_points; ++i)
669  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
670  {
671  adjust_periodicity = true;
672  break;
673  }
674  if (adjust_periodicity == true)
675  {
676  modified_points.resize(surrounding_points.size());
677  std::copy(surrounding_points.begin(),
678  surrounding_points.end(),
679  modified_points.begin());
680  for (unsigned int d = 0; d < spacedim; ++d)
681  if (periodicity[d] > 0)
682  for (unsigned int i = 0; i < n_points; ++i)
683  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
684  modified_points[i][d] -= periodicity[d];
685  surrounding_points_start = modified_points.data();
686  }
687 
688  // Now perform the interpolation
689  for (unsigned int row = 0; row < weights.size(0); ++row)
690  {
691  Assert(
692  std::abs(
693  std::accumulate(&weights(row, 0), &weights(row, 0) + n_points, 0.0) -
694  1.0) < 1e-10,
695  ExcMessage("The weights for the individual points should sum to 1!"));
696  Point<spacedim> new_point;
697  for (unsigned int p = 0; p < n_points; ++p)
698  new_point += surrounding_points_start[p] * weights(row, p);
699 
700  // if necessary, also adjust the weighted point by the periodicity
701  for (unsigned int d = 0; d < spacedim; ++d)
702  if (periodicity[d] > 0)
703  if (new_point[d] < 0)
704  new_point[d] += periodicity[d];
705 
706  // TODO should this use surrounding_points_start or surrounding_points?
707  // The older version used surrounding_points
708  new_points[row] =
709  project_to_manifold(make_array_view(surrounding_points.begin(),
710  surrounding_points.end()),
711  new_point);
712  }
713 }
714 
715 
716 
717 template <int dim, int spacedim>
720  const ArrayView<const Point<spacedim>> & /*vertices*/,
721  const Point<spacedim> &candidate) const
722 {
723  return candidate;
724 }
725 
726 
727 
728 template <int dim, int spacedim>
729 const Tensor<1, spacedim> &
731 {
732  return periodicity;
733 }
734 
735 
736 
737 template <int dim, int spacedim>
740  const Point<spacedim> &x2) const
741 {
742  Tensor<1, spacedim> direction = x2 - x1;
743 
744  // see if we have to take into account periodicity. if so, we need
745  // to make sure that if a distance in one coordinate direction
746  // is larger than half of the box length, then go the other way
747  // around (i.e., via the periodic box)
748  for (unsigned int d = 0; d < spacedim; ++d)
749  if (periodicity[d] > tolerance)
750  {
751  if (direction[d] < -periodicity[d] / 2)
752  direction[d] += periodicity[d];
753  else if (direction[d] > periodicity[d] / 2)
754  direction[d] -= periodicity[d];
755  }
756 
757  return direction;
758 }
759 
760 
761 
762 template <>
763 void
767 {
768  Assert(false, ExcImpossibleInDim(1));
769 }
770 
771 
772 
773 template <>
774 void
778 {
779  Assert(false, ExcNotImplemented());
780 }
781 
782 
783 
784 template <>
785 void
789 {
790  Assert(false, ExcNotImplemented());
791 }
792 
793 
794 
795 template <>
796 void
799  Manifold<2, 2>::FaceVertexNormals & face_vertex_normals) const
800 {
801  const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
802  // We're in 2d. Faces are edges:
803  for (const unsigned int vertex : ReferenceCells::Line.vertex_indices())
804  // compute normals from tangent
805  face_vertex_normals[vertex] = Tensor<1, 2>({tangent[1], -tangent[0]});
806 }
807 
808 
809 
810 template <>
811 void
813  const Triangulation<2, 3>::face_iterator & /*face*/,
814  Manifold<2, 3>::FaceVertexNormals & /*face_vertex_normals*/) const
815 {
816  Assert(false, ExcNotImplemented());
817 }
818 
819 
820 
821 template <>
822 void
825  Manifold<3, 3>::FaceVertexNormals & face_vertex_normals) const
826 {
827  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
828 
829  static const unsigned int neighboring_vertices[4][2] = {{1, 2},
830  {3, 0},
831  {0, 3},
832  {2, 1}};
833  for (unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
834  {
835  // first define the two tangent vectors at the vertex by using the
836  // two lines radiating away from this vertex
837  const Tensor<1, 3> tangents[2] = {
838  face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
839  face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
840 
841  // then compute the normal by taking the cross product. since the
842  // normal is not required to be normalized, no problem here
843  face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
844  }
845 }
846 
847 
848 
849 template <>
852  const Point<1> &) const
853 {
854  Assert(false, ExcNotImplemented());
855  return {};
856 }
857 
858 
859 
860 template <>
863  const Point<2> &) const
864 {
865  Assert(false, ExcNotImplemented());
866  return {};
867 }
868 
869 
870 
871 template <>
874  const Point<3> &) const
875 {
876  Assert(false, ExcNotImplemented());
877  return {};
878 }
879 
880 
881 
882 template <>
886  const Point<2> & p) const
887 {
888  // In 2d, a face is just a straight line and
889  // we can use the 'standard' implementation.
890  return Manifold<2, 2>::normal_vector(face, p);
891 }
892 
893 
894 
895 template <int dim, int spacedim>
898  const typename Triangulation<dim, spacedim>::face_iterator &face,
899  const Point<spacedim> & p) const
900 {
901  // I don't think the implementation below will work when dim!=spacedim;
902  // in fact, I believe that we don't even have enough information here,
903  // because we would need to know not only about the tangent vectors
904  // of the face, but also of the cell, to compute the normal vector.
905  // Someone will have to think about this some more.
906  Assert(dim == spacedim, ExcNotImplemented());
907 
908  // in order to find out what the normal vector is, we first need to
909  // find the reference coordinates of the point p on the given face,
910  // or at least the reference coordinates of the closest point on the
911  // face
912  //
913  // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
914  // where F(xi) is the mapping. this algorithm is implemented in
915  // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
916  // while we need it for faces here. it's also implemented in somewhat
917  // more generality there using the machinery of the MappingQ1 class
918  // while we really only need it for a specific case here
919  //
920  // in any case, the iteration we use here is a Gauss-Newton's iteration with
921  // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
922  // where
923  // J(xi) = (grad F(xi))^T (F(xi)-p)
924  // and
925  // H(xi) = [grad F(xi)]^T [grad F(xi)]
926  // In all this,
927  // F(xi) = sum_v vertex[v] phi_v(xi)
928  // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
929 
930  // we start with the point xi=1/2, xi=(1/2,1/2), ...
931  const unsigned int facedim = dim - 1;
932 
933  Point<facedim> xi;
934 
935  const auto face_reference_cell = face->reference_cell();
936 
937  if (face_reference_cell == ReferenceCells::get_hypercube<facedim>())
938  {
939  for (unsigned int i = 0; i < facedim; ++i)
940  xi[i] = 1. / 2;
941  }
942  else
943  {
944  for (unsigned int i = 0; i < facedim; ++i)
945  xi[i] = 1. / 3;
946  }
947 
948  const double eps = 1e-12;
949  Tensor<1, spacedim> grad_F[facedim];
950  unsigned int iteration = 0;
951  while (true)
952  {
954  for (const unsigned int v : face->vertex_indices())
955  F +=
956  face->vertex(v) * face_reference_cell.d_linear_shape_function(xi, v);
957 
958  for (unsigned int i = 0; i < facedim; ++i)
959  {
960  grad_F[i] = 0;
961  for (const unsigned int v : face->vertex_indices())
962  grad_F[i] +=
963  face->vertex(v) *
964  face_reference_cell.d_linear_shape_function_gradient(xi, v)[i];
965  }
966 
968  for (unsigned int i = 0; i < facedim; ++i)
969  for (unsigned int j = 0; j < spacedim; ++j)
970  J[i] += grad_F[i][j] * (F - p)[j];
971 
973  for (unsigned int i = 0; i < facedim; ++i)
974  for (unsigned int j = 0; j < facedim; ++j)
975  for (unsigned int k = 0; k < spacedim; ++k)
976  H[i][j] += grad_F[i][k] * grad_F[j][k];
977 
978  const Tensor<1, facedim> delta_xi = -invert(H) * J;
979  xi += delta_xi;
980  ++iteration;
981 
982  Assert(iteration < 10,
983  ExcMessage("The Newton iteration to find the reference point "
984  "did not converge in 10 iterations. Do you have a "
985  "deformed cell? (See the glossary for a definition "
986  "of what a deformed cell is. You may want to output "
987  "the vertices of your cell."));
988 
989  // It turns out that the check in reference coordinates with an absolute
990  // tolerance can cause a convergence failure of the Newton method as
991  // seen in tests/manifold/flat_manifold_09.cc. To work around this, also
992  // use a convergence check in world coordinates. This check has to be
993  // relative to the size of the face of course. Here we decided to use
994  // diameter because it works for non-planar faces and is cheap to
995  // compute:
996  const double normalized_delta_world = (F - p).norm() / face->diameter();
997 
998  if (delta_xi.norm() < eps || normalized_delta_world < eps)
999  break;
1000  }
1001 
1002  // so now we have the reference coordinates xi of the point p.
1003  // we then have to compute the normal vector, which we can do
1004  // by taking the (normalize) alternating product of all the tangent
1005  // vectors given by grad_F
1006  return internal::normalized_alternating_product(grad_F);
1007 }
1008 
1010 /* -------------------------- ChartManifold --------------------- */
1011 template <int dim, int spacedim, int chartdim>
1013  const Tensor<1, chartdim> &periodicity)
1014  : sub_manifold(periodicity)
1015 {}
1016 
1017 
1018 
1019 template <int dim, int spacedim, int chartdim>
1022  const Point<spacedim> &p1,
1023  const Point<spacedim> &p2,
1024  const double w) const
1025 {
1026  const std::array<Point<spacedim>, 2> points{{p1, p2}};
1027  const std::array<double, 2> weights{{1. - w, w}};
1028  return get_new_point(make_array_view(points.begin(), points.end()),
1029  make_array_view(weights.begin(), weights.end()));
1030 }
1031 
1032 
1033 
1034 template <int dim, int spacedim, int chartdim>
1037  const ArrayView<const Point<spacedim>> &surrounding_points,
1038  const ArrayView<const double> & weights) const
1039 {
1040  const std::size_t n_points = surrounding_points.size();
1041 
1042  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1043 
1044  for (unsigned int i = 0; i < n_points; ++i)
1045  chart_points[i] = pull_back(surrounding_points[i]);
1046 
1047  const Point<chartdim> p_chart = sub_manifold.get_new_point(
1048  make_array_view(chart_points.begin(), chart_points.end()), weights);
1049 
1050  return push_forward(p_chart);
1051 }
1052 
1053 
1054 
1055 template <int dim, int spacedim, int chartdim>
1056 void
1058  const ArrayView<const Point<spacedim>> &surrounding_points,
1059  const Table<2, double> & weights,
1060  ArrayView<Point<spacedim>> new_points) const
1061 {
1062  Assert(weights.size(0) > 0, ExcEmptyObject());
1063  AssertDimension(surrounding_points.size(), weights.size(1));
1064 
1065  const std::size_t n_points = surrounding_points.size();
1066 
1067  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1068  for (std::size_t i = 0; i < n_points; ++i)
1069  chart_points[i] = pull_back(surrounding_points[i]);
1070 
1071  boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1072  weights.size(0));
1073  sub_manifold.get_new_points(
1074  make_array_view(chart_points.begin(), chart_points.end()),
1075  weights,
1076  make_array_view(new_points_on_chart.begin(), new_points_on_chart.end()));
1077 
1078  for (std::size_t row = 0; row < weights.size(0); ++row)
1079  new_points[row] = push_forward(new_points_on_chart[row]);
1080 }
1081 
1082 
1083 
1084 template <int dim, int spacedim, int chartdim>
1087  const Point<chartdim> &) const
1088 {
1089  // function must be implemented in a derived class to be usable,
1090  // as discussed in this function's documentation
1091  Assert(false, ExcPureFunctionCalled());
1093 }
1094 
1095 
1096 
1097 template <int dim, int spacedim, int chartdim>
1100  const Point<spacedim> &x1,
1101  const Point<spacedim> &x2) const
1102 {
1104  push_forward_gradient(pull_back(x1));
1105 
1106  // ensure that the chart is not singular by asserting that its
1107  // derivative has a positive determinant. we need to make this
1108  // comparison relative to the size of the derivative. since the
1109  // determinant is the product of chartdim factors, take the
1110  // chartdim-th root of it in comparing against the size of the
1111  // derivative
1112  Assert(std::pow(std::abs(F_prime.determinant()), 1. / chartdim) >=
1113  1e-12 * F_prime.norm(),
1114  ExcMessage(
1115  "The derivative of a chart function must not be singular."));
1116 
1117  const Tensor<1, chartdim> delta =
1118  sub_manifold.get_tangent_vector(pull_back(x1), pull_back(x2));
1119 
1120  Tensor<1, spacedim> result;
1121  for (unsigned int i = 0; i < spacedim; ++i)
1122  result[i] += F_prime[i] * delta;
1123 
1124  return result;
1125 }
1126 
1127 
1128 
1129 template <int dim, int spacedim, int chartdim>
1130 const Tensor<1, chartdim> &
1132 {
1133  return sub_manifold.get_periodicity();
1134 }
1135 
1136 // explicit instantiations
1137 #include "manifold.inst"
1138 
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:1036
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:1012
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1099
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1086
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1131
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:1057
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:1021
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:797
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:811
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
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
Definition: manifold.cc:65
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &surrounding_points, const Point< spacedim > &candidate) const
Definition: manifold.cc:40
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:339
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:368
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:488
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:275
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:52
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:382
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:514
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:402
virtual Point< spacedim > get_new_point_on_line(const typename Triangulation< dim, spacedim >::line_iterator &line) const
Definition: manifold.cc:354
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
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
Point< 3 > vertices[4]
unsigned int vertex_indices[2]
Definition: grid_tools.cc:1341
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcPureFunctionCalled()
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcMessage(std::string arg1)
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1503
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1479
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:832
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:206
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)