Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 ((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  p = get_intermediate_point(p,
105  surrounding_points[permutation[i]],
106  1.0 - weight);
107  w += weights[permutation[i]];
108  }
109 
110  return p;
111 }
112 
113 
114 
115 template <int dim, int spacedim>
116 void
118  const ArrayView<const Point<spacedim>> &surrounding_points,
119  const Table<2, double> & weights,
120  ArrayView<Point<spacedim>> new_points) const
121 {
122  AssertDimension(surrounding_points.size(), weights.size(1));
123 
124  for (unsigned int row = 0; row < weights.size(0); ++row)
125  {
126  new_points[row] =
127  get_new_point(make_array_view(surrounding_points.begin(),
128  surrounding_points.end()),
129  make_array_view(weights, row));
130  }
131 }
132 
133 
134 
135 template <>
138  const Point<2> & p) const
139 {
140  const int spacedim = 2;
141 
142  // get the tangent vector from the point 'p' in the direction of the further
143  // one of the two vertices that make up the face of this 2d cell
144  const Tensor<1, spacedim> tangent =
145  ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
146  -get_tangent_vector(p, face->vertex(0)) :
147  get_tangent_vector(p, face->vertex(1)));
148 
149  // then rotate it by 90 degrees
150  const Tensor<1, spacedim> normal = cross_product_2d(tangent);
151  return normal / normal.norm();
152 }
153 
154 
155 
156 template <>
159  const Point<3> & p) const
160 {
161  const int spacedim = 3;
162 
163  const std::array<Point<spacedim>, 4> vertices{
164  {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
165  const std::array<double, 4> distances{{vertices[0].distance(p),
166  vertices[1].distance(p),
167  vertices[2].distance(p),
168  vertices[3].distance(p)}};
169  const double max_distance = std::max(std::max(distances[0], distances[1]),
170  std::max(distances[2], distances[3]));
171 
172  // We need to find two tangential vectors to the given point p, but we do
173  // not know how the point is oriented against the face. We guess the two
174  // directions by assuming a flat topology and take the two directions that
175  // indicate the angle closest to a perpendicular one (i.e., cos(theta) close
176  // to zero). We start with an invalid value but the loops should always find
177  // a value.
178  double abs_cos_angle = std::numeric_limits<double>::max();
179  unsigned int first_index = numbers::invalid_unsigned_int,
180  second_index = numbers::invalid_unsigned_int;
181  for (unsigned int i = 0; i < 3; ++i)
182  if (distances[i] > 1e-8 * max_distance)
183  for (unsigned int j = i + 1; j < 4; ++j)
184  if (distances[j] > 1e-8 * max_distance)
185  {
186  const double new_angle = (p - vertices[i]) * (p - vertices[j]) /
187  (distances[i] * distances[j]);
188  // multiply by factor 0.999 to bias the search in a way that
189  // avoids trouble with roundoff
190  if (std::abs(new_angle) < 0.999 * abs_cos_angle)
191  {
192  abs_cos_angle = std::abs(new_angle);
193  first_index = i;
194  second_index = j;
195  }
196  }
197  Assert(first_index != numbers::invalid_unsigned_int,
198  ExcMessage("The search for possible directions did not succeed."));
199 
200  // Compute tangents and normal for selected vertices
201  Tensor<1, spacedim> t1 = get_tangent_vector(p, vertices[first_index]);
202  Tensor<1, spacedim> t2 = get_tangent_vector(p, vertices[second_index]);
203  Tensor<1, spacedim> normal = cross_product_3d(t1, t2);
204 
205  Assert(
206  normal.norm_square() > 1e4 * std::numeric_limits<double>::epsilon() *
207  t1.norm_square() * t2.norm_square(),
208  ExcMessage(
209  "Manifold::normal_vector was unable to find a suitable combination "
210  "of vertices to compute a normal on this face. We chose the secants "
211  "that are as orthogonal as possible, but tangents appear to be "
212  "linearly dependent. Check for distorted faces in your triangulation."));
213 
214  // Now figure out if we need to flip the direction, we do this by comparing
215  // to a reference normal that would be the correct result if all vertices
216  // would lie in a plane
217  const Tensor<1, spacedim> rt1 = vertices[3] - vertices[0];
218  const Tensor<1, spacedim> rt2 = vertices[2] - vertices[1];
219  const Tensor<1, spacedim> reference_normal = cross_product_3d(rt1, rt2);
220 
221  if (reference_normal * normal < 0.0)
222  normal *= -1.0;
223 
224  return normal / normal.norm();
225 }
226 
227 
228 
229 template <int dim, int spacedim>
232  const typename Triangulation<dim, spacedim>::face_iterator & /*face*/,
233  const Point<spacedim> & /*p*/) const
234 {
235  Assert(false, ExcPureFunctionCalled());
236  return Tensor<1, spacedim>();
237 }
238 
239 
240 
241 template <>
242 void
245  FaceVertexNormals & n) const
246 {
247  n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
248  n[1] =
249  -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
250 
251  for (unsigned int i = 0; i < 2; ++i)
252  {
253  Assert(n[i].norm() != 0,
254  ExcInternalError("Something went wrong. The "
255  "computed normals have "
256  "zero length."));
257  n[i] /= n[i].norm();
258  }
259 }
260 
261 
262 
263 template <>
264 void
267  FaceVertexNormals & n) const
268 {
269  n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
270  get_tangent_vector(face->vertex(0), face->vertex(2)));
271 
272  n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
273  get_tangent_vector(face->vertex(1), face->vertex(0)));
274 
275  n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
276  get_tangent_vector(face->vertex(2), face->vertex(3)));
277 
278  n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
279  get_tangent_vector(face->vertex(3), face->vertex(1)));
280 
281  for (unsigned int i = 0; i < 4; ++i)
282  {
283  Assert(n[i].norm() != 0,
284  ExcInternalError("Something went wrong. The "
285  "computed normals have "
286  "zero length."));
287  n[i] /= n[i].norm();
288  }
289 }
290 
291 
292 
293 template <int dim, int spacedim>
294 void
296  const typename Triangulation<dim, spacedim>::face_iterator &face,
297  FaceVertexNormals & n) const
298 {
299  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
300  {
301  n[v] = normal_vector(face, face->vertex(v));
302  n[v] /= n[v].norm();
303  }
304 }
305 
306 
307 
308 template <int dim, int spacedim>
311  const typename Triangulation<dim, spacedim>::line_iterator &line) const
312 {
313  const auto points_weights = get_default_points_and_weights(line);
314  return get_new_point(make_array_view(points_weights.first.begin(),
315  points_weights.first.end()),
316  make_array_view(points_weights.second.begin(),
317  points_weights.second.end()));
318 }
319 
320 
321 
322 template <int dim, int spacedim>
325  const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
326 {
327  const auto points_weights = get_default_points_and_weights(quad);
328  return get_new_point(make_array_view(points_weights.first.begin(),
329  points_weights.first.end()),
330  make_array_view(points_weights.second.begin(),
331  points_weights.second.end()));
332 }
333 
334 
335 
336 template <int dim, int spacedim>
339  const typename Triangulation<dim, spacedim>::face_iterator &face) const
340 {
341  Assert(dim > 1, ExcImpossibleInDim(dim));
342 
343  switch (dim)
344  {
345  case 2:
346  return get_new_point_on_line(face);
347  case 3:
348  return get_new_point_on_quad(face);
349  }
350 
351  return Point<spacedim>();
352 }
353 
354 
355 
356 template <int dim, int spacedim>
359  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
360 {
361  switch (dim)
362  {
363  case 1:
364  return get_new_point_on_line(cell);
365  case 2:
366  return get_new_point_on_quad(cell);
367  case 3:
368  return get_new_point_on_hex(cell);
369  }
370 
371  return Point<spacedim>();
372 }
373 
374 
375 
376 template <>
377 Point<1>
380 {
381  Assert(false, ExcImpossibleInDim(1));
382  return {};
383 }
384 
385 
386 
387 template <>
388 Point<2>
391 {
392  Assert(false, ExcImpossibleInDim(1));
393  return {};
394 }
395 
396 
397 
398 template <>
399 Point<3>
402 {
403  Assert(false, ExcImpossibleInDim(1));
404  return {};
405 }
406 
407 
408 
409 template <>
410 Point<1>
413 {
414  Assert(false, ExcImpossibleInDim(1));
415  return {};
416 }
417 
418 
419 
420 template <>
421 Point<2>
424 {
425  Assert(false, ExcImpossibleInDim(1));
426  return {};
427 }
428 
429 
430 
431 template <>
432 Point<3>
435 {
436  Assert(false, ExcImpossibleInDim(1));
437  return {};
438 }
439 
440 
441 
442 template <int dim, int spacedim>
445  const typename Triangulation<dim, spacedim>::hex_iterator & /*hex*/) const
446 {
447  Assert(false, ExcImpossibleInDim(dim));
448  return Point<spacedim>();
449 }
450 
451 
452 
453 template <>
454 Point<3>
456  const Triangulation<3, 3>::hex_iterator &hex) const
457 {
458  const auto points_weights = get_default_points_and_weights(hex, true);
459  return get_new_point(make_array_view(points_weights.first.begin(),
460  points_weights.first.end()),
461  make_array_view(points_weights.second.begin(),
462  points_weights.second.end()));
463 }
464 
465 
466 
467 template <int dim, int spacedim>
470  const Point<spacedim> &x2) const
471 {
472  const double epsilon = 1e-8;
473 
474  const std::array<Point<spacedim>, 2> points{{x1, x2}};
475  const std::array<double, 2> weights{{epsilon, 1.0 - epsilon}};
476  const Point<spacedim> neighbor_point =
477  get_new_point(make_array_view(points.begin(), points.end()),
478  make_array_view(weights.begin(), weights.end()));
479  return (neighbor_point - x1) / epsilon;
480 }
481 
482 /* -------------------------- FlatManifold --------------------- */
483 
484 namespace internal
485 {
486  namespace
487  {
489  normalized_alternating_product(const Tensor<1, 3> (&)[1])
490  {
491  // we get here from FlatManifold<2,3>::normal_vector, but
492  // the implementation below is bogus for this case anyway
493  // (see the assert at the beginning of that function).
494  Assert(false, ExcNotImplemented());
495  return {};
496  }
497 
498 
499 
501  normalized_alternating_product(const Tensor<1, 3> (&basis_vectors)[2])
502  {
503  Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
504  return tmp / tmp.norm();
505  }
506 
507  } // namespace
508 } // namespace internal
509 
510 template <int dim, int spacedim>
512  const Tensor<1, spacedim> &periodicity,
513  const double tolerance)
514  : periodicity(periodicity)
515  , tolerance(tolerance)
516 {}
517 
518 
519 
520 template <int dim, int spacedim>
521 std::unique_ptr<Manifold<dim, spacedim>>
523 {
524  return std_cxx14::make_unique<FlatManifold<dim, spacedim>>(periodicity,
525  tolerance);
526 }
527 
528 
529 
530 template <int dim, int spacedim>
533  const ArrayView<const Point<spacedim>> &surrounding_points,
534  const ArrayView<const double> & weights) const
535 {
536  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
537  1e-10,
538  ExcMessage("The weights for the individual points should sum to 1!"));
539 
540  Point<spacedim> p;
541 
542  // if there is no periodicity, use a shortcut
543  if (periodicity == Tensor<1, spacedim>())
544  {
545  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
546  p += surrounding_points[i] * weights[i];
547  }
548  else
549  {
550  Tensor<1, spacedim> minP = periodicity;
551 
552  for (unsigned int d = 0; d < spacedim; ++d)
553  if (periodicity[d] > 0)
554  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
555  {
556  minP[d] = std::min(minP[d], surrounding_points[i][d]);
557  Assert((surrounding_points[i][d] <
558  periodicity[d] + tolerance * periodicity[d]) ||
559  (surrounding_points[i][d] >=
560  -tolerance * periodicity[d]),
561  ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
562  }
563 
564  // compute the weighted average point, possibly taking into account
565  // periodicity
566  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
567  {
568  Point<spacedim> dp;
569  for (unsigned int d = 0; d < spacedim; ++d)
570  if (periodicity[d] > 0)
571  dp[d] =
572  ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0 ?
573  -periodicity[d] :
574  0.0);
575 
576  p += (surrounding_points[i] + dp) * weights[i];
577  }
578 
579  // if necessary, also adjust the weighted point by the periodicity
580  for (unsigned int d = 0; d < spacedim; ++d)
581  if (periodicity[d] > 0)
582  if (p[d] < 0)
583  p[d] += periodicity[d];
584  }
585 
586  return project_to_manifold(surrounding_points, p);
587 }
588 
589 
590 
591 template <int dim, int spacedim>
592 void
594  const ArrayView<const Point<spacedim>> &surrounding_points,
595  const Table<2, double> & weights,
596  ArrayView<Point<spacedim>> new_points) const
597 {
598  AssertDimension(surrounding_points.size(), weights.size(1));
599  if (weights.size(0) == 0)
600  return;
601 
602  const std::size_t n_points = surrounding_points.size();
603 
604  Tensor<1, spacedim> minP = periodicity;
605  for (unsigned int d = 0; d < spacedim; ++d)
606  if (periodicity[d] > 0)
607  for (unsigned int i = 0; i < n_points; ++i)
608  {
609  minP[d] = std::min(minP[d], surrounding_points[i][d]);
610  Assert((surrounding_points[i][d] <
611  periodicity[d] + tolerance * periodicity[d]) ||
612  (surrounding_points[i][d] >= -tolerance * periodicity[d]),
613  ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
614  }
615 
616  // check whether periodicity shifts some of the points. Only do this if
617  // necessary to avoid memory allocation
618  const Point<spacedim> *surrounding_points_start = surrounding_points.data();
619 
620  boost::container::small_vector<Point<spacedim>, 200> modified_points;
621  bool adjust_periodicity = false;
622  for (unsigned int d = 0; d < spacedim; ++d)
623  if (periodicity[d] > 0)
624  for (unsigned int i = 0; i < n_points; ++i)
625  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
626  {
627  adjust_periodicity = true;
628  break;
629  }
630  if (adjust_periodicity == true)
631  {
632  modified_points.resize(surrounding_points.size());
633  std::copy(surrounding_points.begin(),
634  surrounding_points.end(),
635  modified_points.begin());
636  for (unsigned int d = 0; d < spacedim; ++d)
637  if (periodicity[d] > 0)
638  for (unsigned int i = 0; i < n_points; ++i)
639  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
640  modified_points[i][d] -= periodicity[d];
641  surrounding_points_start = modified_points.data();
642  }
643 
644  // Now perform the interpolation
645  for (unsigned int row = 0; row < weights.size(0); ++row)
646  {
647  Assert(
648  std::abs(
649  std::accumulate(&weights(row, 0), &weights(row, 0) + n_points, 0.0) -
650  1.0) < 1e-10,
651  ExcMessage("The weights for the individual points should sum to 1!"));
652  Point<spacedim> new_point;
653  for (unsigned int p = 0; p < n_points; ++p)
654  new_point += surrounding_points_start[p] * weights(row, p);
655 
656  // if necessary, also adjust the weighted point by the periodicity
657  for (unsigned int d = 0; d < spacedim; ++d)
658  if (periodicity[d] > 0)
659  if (new_point[d] < 0)
660  new_point[d] += periodicity[d];
661 
662  // TODO should this use surrounding_points_start or surrounding_points?
663  // The older version used surrounding_points
664  new_points[row] =
665  project_to_manifold(make_array_view(surrounding_points.begin(),
666  surrounding_points.end()),
667  new_point);
668  }
669 }
670 
671 
672 
673 template <int dim, int spacedim>
676  const ArrayView<const Point<spacedim>> & /*vertices*/,
677  const Point<spacedim> &candidate) const
678 {
679  return candidate;
680 }
681 
682 
683 
684 template <int dim, int spacedim>
685 const Tensor<1, spacedim> &
687 {
688  return periodicity;
689 }
690 
691 
692 
693 template <int dim, int spacedim>
696  const Point<spacedim> &x2) const
697 {
698  Tensor<1, spacedim> direction = x2 - x1;
699 
700  // see if we have to take into account periodicity. if so, we need
701  // to make sure that if a distance in one coordinate direction
702  // is larger than half of the box length, then go the other way
703  // around (i.e., via the periodic box)
704  for (unsigned int d = 0; d < spacedim; ++d)
705  if (periodicity[d] > tolerance)
706  {
707  if (direction[d] < -periodicity[d] / 2)
708  direction[d] += periodicity[d];
709  else if (direction[d] > periodicity[d] / 2)
710  direction[d] -= periodicity[d];
711  }
712 
713  return direction;
714 }
715 
716 
717 
718 template <>
719 void
723 {
724  Assert(false, ExcImpossibleInDim(1));
725 }
726 
727 
728 
729 template <>
730 void
734 {
735  Assert(false, ExcNotImplemented());
736 }
737 
738 
739 
740 template <>
741 void
745 {
746  Assert(false, ExcNotImplemented());
747 }
748 
749 
750 
751 template <>
752 void
755  Manifold<2, 2>::FaceVertexNormals & face_vertex_normals) const
756 {
757  const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
758  for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_face;
759  ++vertex)
760  // compute normals from tangent
761  face_vertex_normals[vertex] = Point<2>(tangent[1], -tangent[0]);
762 }
763 
764 
765 
766 template <>
767 void
769  const Triangulation<2, 3>::face_iterator & /*face*/,
770  Manifold<2, 3>::FaceVertexNormals & /*face_vertex_normals*/) const
771 {
772  Assert(false, ExcNotImplemented());
773 }
774 
775 
776 
777 template <>
778 void
781  Manifold<3, 3>::FaceVertexNormals & face_vertex_normals) const
782 {
783  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
784 
785  static const unsigned int neighboring_vertices[4][2] = {{1, 2},
786  {3, 0},
787  {0, 3},
788  {2, 1}};
789  for (unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
790  {
791  // first define the two tangent vectors at the vertex by using the
792  // two lines radiating away from this vertex
793  const Tensor<1, 3> tangents[2] = {
794  face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
795  face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
796 
797  // then compute the normal by taking the cross product. since the
798  // normal is not required to be normalized, no problem here
799  face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
800  }
801 }
802 
803 
804 
805 template <>
808  const Point<1> &) const
809 {
810  Assert(false, ExcNotImplemented());
811  return {};
812 }
813 
814 
815 
816 template <>
819  const Point<2> &) const
820 {
821  Assert(false, ExcNotImplemented());
822  return {};
823 }
824 
825 
826 
827 template <>
830  const Point<3> &) const
831 {
832  Assert(false, ExcNotImplemented());
833  return {};
834 }
835 
836 
837 
838 template <>
842  const Point<2> & p) const
843 {
844  // In 2d, a face is just a straight line and
845  // we can use the 'standard' implementation.
846  return Manifold<2, 2>::normal_vector(face, p);
847 }
848 
849 
850 
851 template <int dim, int spacedim>
854  const typename Triangulation<dim, spacedim>::face_iterator &face,
855  const Point<spacedim> & p) const
856 {
857  // I don't think the implementation below will work when dim!=spacedim;
858  // in fact, I believe that we don't even have enough information here,
859  // because we would need to know not only about the tangent vectors
860  // of the face, but also of the cell, to compute the normal vector.
861  // Someone will have to think about this some more.
862  Assert(dim == spacedim, ExcNotImplemented());
863 
864  // in order to find out what the normal vector is, we first need to
865  // find the reference coordinates of the point p on the given face,
866  // or at least the reference coordinates of the closest point on the
867  // face
868  //
869  // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
870  // where F(xi) is the mapping. this algorithm is implemented in
871  // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
872  // while we need it for faces here. it's also implemented in somewhat
873  // more generality there using the machinery of the MappingQ1 class
874  // while we really only need it for a specific case here
875  //
876  // in any case, the iteration we use here is a Gauss-Newton's iteration with
877  // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
878  // where
879  // J(xi) = (grad F(xi))^T (F(xi)-p)
880  // and
881  // H(xi) = [grad F(xi)]^T [grad F(xi)]
882  // In all this,
883  // F(xi) = sum_v vertex[v] phi_v(xi)
884  // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
885 
886  // we start with the point xi=1/2, xi=(1/2,1/2), ...
887  const unsigned int facedim = dim - 1;
888 
889  Point<facedim> xi;
890  for (unsigned int i = 0; i < facedim; ++i)
891  xi[i] = 1. / 2;
892 
893  const double eps = 1e-12;
894  Tensor<1, spacedim> grad_F[facedim];
895  unsigned int iteration = 0;
896  while (true)
897  {
898  Point<spacedim> F;
899  for (unsigned int v = 0; v < GeometryInfo<facedim>::vertices_per_cell;
900  ++v)
901  F += face->vertex(v) *
903 
904  for (unsigned int i = 0; i < facedim; ++i)
905  {
906  grad_F[i] = 0;
907  for (unsigned int v = 0; v < GeometryInfo<facedim>::vertices_per_cell;
908  ++v)
909  grad_F[i] +=
910  face->vertex(v) *
912  }
913 
915  for (unsigned int i = 0; i < facedim; ++i)
916  for (unsigned int j = 0; j < spacedim; ++j)
917  J[i] += grad_F[i][j] * (F - p)[j];
918 
920  for (unsigned int i = 0; i < facedim; ++i)
921  for (unsigned int j = 0; j < facedim; ++j)
922  for (unsigned int k = 0; k < spacedim; ++k)
923  H[i][j] += grad_F[i][k] * grad_F[j][k];
924 
925  const Tensor<1, facedim> delta_xi = -invert(H) * J;
926  xi += delta_xi;
927  ++iteration;
928 
929  Assert(iteration < 10,
930  ExcMessage("The Newton iteration to find the reference point "
931  "did not converge in 10 iterations. Do you have a "
932  "deformed cell? (See the glossary for a definition "
933  "of what a deformed cell is. You may want to output "
934  "the vertices of your cell."));
935 
936  // It turns out that the check in reference coordinates with an absolute
937  // tolerance can cause a convergence failure of the Newton method as
938  // seen in tests/manifold/flat_manifold_09.cc. To work around this, also
939  // use a convergence check in world coordinates. This check has to be
940  // relative to the size of the face of course. Here we decided to use
941  // diameter because it works for non-planar faces and is cheap to
942  // compute:
943  const double normalized_delta_world = (F - p).norm() / face->diameter();
944 
945  if (delta_xi.norm() < eps || normalized_delta_world < eps)
946  break;
947  }
948 
949  // so now we have the reference coordinates xi of the point p.
950  // we then have to compute the normal vector, which we can do
951  // by taking the (normalize) alternating product of all the tangent
952  // vectors given by grad_F
953  return internal::normalized_alternating_product(grad_F);
954 }
955 
956 
957 /* -------------------------- ChartManifold --------------------- */
958 template <int dim, int spacedim, int chartdim>
960  const Tensor<1, chartdim> &periodicity)
961  : sub_manifold(periodicity)
962 {}
963 
964 
965 
966 template <int dim, int spacedim, int chartdim>
969  const Point<spacedim> &p1,
970  const Point<spacedim> &p2,
971  const double w) const
972 {
973  const std::array<Point<spacedim>, 2> points{{p1, p2}};
974  const std::array<double, 2> weights{{1. - w, w}};
975  return get_new_point(make_array_view(points.begin(), points.end()),
976  make_array_view(weights.begin(), weights.end()));
977 }
978 
979 
980 
981 template <int dim, int spacedim, int chartdim>
984  const ArrayView<const Point<spacedim>> &surrounding_points,
985  const ArrayView<const double> & weights) const
986 {
987  const std::size_t n_points = surrounding_points.size();
988 
989  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
990 
991  for (unsigned int i = 0; i < n_points; ++i)
992  chart_points[i] = pull_back(surrounding_points[i]);
993 
994  const Point<chartdim> p_chart = sub_manifold.get_new_point(
995  make_array_view(chart_points.begin(), chart_points.end()), weights);
996 
997  return push_forward(p_chart);
998 }
999 
1000 
1001 
1002 template <int dim, int spacedim, int chartdim>
1003 void
1005  const ArrayView<const Point<spacedim>> &surrounding_points,
1006  const Table<2, double> & weights,
1007  ArrayView<Point<spacedim>> new_points) const
1008 {
1009  Assert(weights.size(0) > 0, ExcEmptyObject());
1010  AssertDimension(surrounding_points.size(), weights.size(1));
1011 
1012  const std::size_t n_points = surrounding_points.size();
1013 
1014  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1015  for (std::size_t i = 0; i < n_points; ++i)
1016  chart_points[i] = pull_back(surrounding_points[i]);
1017 
1018  boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1019  weights.size(0));
1020  sub_manifold.get_new_points(
1021  make_array_view(chart_points.begin(), chart_points.end()),
1022  weights,
1023  make_array_view(new_points_on_chart.begin(), new_points_on_chart.end()));
1024 
1025  for (std::size_t row = 0; row < weights.size(0); ++row)
1026  new_points[row] = push_forward(new_points_on_chart[row]);
1027 }
1028 
1029 
1030 
1031 template <int dim, int spacedim, int chartdim>
1034  const Point<chartdim> &) const
1035 {
1036  // function must be implemented in a derived class to be usable,
1037  // as discussed in this function's documentation
1038  Assert(false, ExcPureFunctionCalled());
1040 }
1041 
1042 
1043 
1044 template <int dim, int spacedim, int chartdim>
1047  const Point<spacedim> &x1,
1048  const Point<spacedim> &x2) const
1049 {
1051  push_forward_gradient(pull_back(x1));
1052 
1053  // ensure that the chart is not singular by asserting that its
1054  // derivative has a positive determinant. we need to make this
1055  // comparison relative to the size of the derivative. since the
1056  // determinant is the product of chartdim factors, take the
1057  // chartdim-th root of it in comparing against the size of the
1058  // derivative
1059  Assert(std::pow(std::abs(F_prime.determinant()), 1. / chartdim) >=
1060  1e-12 * F_prime.norm(),
1061  ExcMessage(
1062  "The derivative of a chart function must not be singular."));
1063 
1064  const Tensor<1, chartdim> delta =
1065  sub_manifold.get_tangent_vector(pull_back(x1), pull_back(x2));
1066 
1067  Tensor<1, spacedim> result;
1068  for (unsigned int i = 0; i < spacedim; ++i)
1069  result[i] += F_prime[i] * delta;
1070 
1071  return result;
1072 }
1073 
1074 
1075 
1076 template <int dim, int spacedim, int chartdim>
1077 const Tensor<1, chartdim> &
1079 {
1080  return sub_manifold.get_periodicity();
1081 }
1082 
1083 // explicit instantiations
1084 #include "manifold.inst"
1085 
1086 DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcPureFunctionCalled()
static const unsigned int invalid_unsigned_int
Definition: types.h:173
FlatManifold(const Tensor< 1, spacedim > &periodicity=Tensor< 1, spacedim >(), const double tolerance=1e-10)
Definition: manifold.cc:511
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:853
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1592
SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
const Tensor< 1, spacedim > & get_periodicity() const
Definition: manifold.cc:686
virtual Point< spacedim > get_new_point_on_hex(const typename Triangulation< dim, spacedim >::hex_iterator &hex) const
Definition: manifold.cc:444
typename IteratorSelector::hex_iterator hex_iterator
Definition: tria.h:1640
iterator end() const
Definition: array_view.h:489
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const
Definition: manifold.cc:63
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
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:1318
Number determinant() const
std::size_t size() const
Definition: array_view.h:471
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:310
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const
Definition: manifold.cc:469
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition: manifold.cc:1033
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:1407
numbers::NumberTraits< Number >::real_type norm_square() const
Definition: tensor.h:1326
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1078
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
Definition: manifold.cc:675
numbers::NumberTraits< Number >::real_type norm() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold.cc:522
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1616
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const
Definition: manifold.cc:50
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
Tensor< 1, dim, Number > cross_product_3d(const Tensor< 1, dim, Number > &src1, const Tensor< 1, dim, Number > &src2)
Definition: tensor.h:2100
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:1004
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1046
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:117
size_type size(const unsigned int i) const
static double d_linear_shape_function(const Point< dim > &xi, const unsigned int i)
Definition: mpi.h:90
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:695
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:983
virtual Point< spacedim > get_new_point_on_quad(const typename Triangulation< dim, spacedim >::quad_iterator &quad) const
Definition: manifold.cc:324
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:532
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
std::array< Tensor< 1, spacedim >, GeometryInfo< dim >::vertices_per_face > FaceVertexNormals
Definition: manifold.h:350
static ::ExceptionBase & ExcNotImplemented()
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:594
virtual Tensor< 1, spacedim > normal_vector(const typename Triangulation< dim, spacedim >::face_iterator &face, const Point< spacedim > &p) const
Definition: manifold.cc:231
iterator begin() const
Definition: array_view.h:480
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:593
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:959
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:338
virtual void get_normals_at_vertices(const typename Triangulation< dim, spacedim >::face_iterator &face, FaceVertexNormals &face_vertex_normals) const
Definition: manifold.cc:295
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:968
static ::ExceptionBase & ExcInternalError()
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:358