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