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