Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
\(\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 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/table.h>
17 #include <deal.II/base/tensor.h>
18 
19 #include <deal.II/fe/fe_q.h>
20 
21 #include <deal.II/grid/manifold.h>
22 #include <deal.II/grid/tria.h>
25 
26 #include <boost/container/small_vector.hpp>
27 
28 #include <cmath>
29 #include <memory>
30 
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(
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::make_unique<FlatManifold<dim, spacedim>>(periodicity, tolerance);
531 }
532 
533 
534 
535 template <int dim, int spacedim>
538  const ArrayView<const Point<spacedim>> &surrounding_points,
539  const ArrayView<const double> & weights) const
540 {
541  Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
542  1e-10,
543  ExcMessage("The weights for the individual points should sum to 1!"));
544 
545  Point<spacedim> p;
546 
547  // if there is no periodicity, use a shortcut
549  {
550  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
551  p += surrounding_points[i] * weights[i];
552  }
553  else
554  {
556 
557  for (unsigned int d = 0; d < spacedim; ++d)
558  if (periodicity[d] > 0)
559  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
560  {
561  minP[d] = std::min(minP[d], surrounding_points[i][d]);
562  Assert((surrounding_points[i][d] <
563  periodicity[d] + tolerance * periodicity[d]) ||
564  (surrounding_points[i][d] >=
565  -tolerance * periodicity[d]),
566  ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
567  }
568 
569  // compute the weighted average point, possibly taking into account
570  // periodicity
571  for (unsigned int i = 0; i < surrounding_points.size(); ++i)
572  {
573  Point<spacedim> dp;
574  for (unsigned int d = 0; d < spacedim; ++d)
575  if (periodicity[d] > 0)
576  dp[d] =
577  ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0 ?
578  -periodicity[d] :
579  0.0);
580 
581  p += (surrounding_points[i] + dp) * weights[i];
582  }
583 
584  // if necessary, also adjust the weighted point by the periodicity
585  for (unsigned int d = 0; d < spacedim; ++d)
586  if (periodicity[d] > 0)
587  if (p[d] < 0)
588  p[d] += periodicity[d];
589  }
590 
591  return project_to_manifold(surrounding_points, p);
592 }
593 
594 
595 
596 template <int dim, int spacedim>
597 void
599  const ArrayView<const Point<spacedim>> &surrounding_points,
600  const Table<2, double> & weights,
601  ArrayView<Point<spacedim>> new_points) const
602 {
603  AssertDimension(surrounding_points.size(), weights.size(1));
604  if (weights.size(0) == 0)
605  return;
606 
607  const std::size_t n_points = surrounding_points.size();
608 
610  for (unsigned int d = 0; d < spacedim; ++d)
611  if (periodicity[d] > 0)
612  for (unsigned int i = 0; i < n_points; ++i)
613  {
614  minP[d] = std::min(minP[d], surrounding_points[i][d]);
615  Assert((surrounding_points[i][d] <
616  periodicity[d] + tolerance * periodicity[d]) ||
617  (surrounding_points[i][d] >= -tolerance * periodicity[d]),
618  ExcPeriodicBox(d, surrounding_points[i], periodicity[i]));
619  }
620 
621  // check whether periodicity shifts some of the points. Only do this if
622  // necessary to avoid memory allocation
623  const Point<spacedim> *surrounding_points_start = surrounding_points.data();
624 
625  boost::container::small_vector<Point<spacedim>, 200> modified_points;
626  bool adjust_periodicity = false;
627  for (unsigned int d = 0; d < spacedim; ++d)
628  if (periodicity[d] > 0)
629  for (unsigned int i = 0; i < n_points; ++i)
630  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
631  {
632  adjust_periodicity = true;
633  break;
634  }
635  if (adjust_periodicity == true)
636  {
637  modified_points.resize(surrounding_points.size());
638  std::copy(surrounding_points.begin(),
639  surrounding_points.end(),
640  modified_points.begin());
641  for (unsigned int d = 0; d < spacedim; ++d)
642  if (periodicity[d] > 0)
643  for (unsigned int i = 0; i < n_points; ++i)
644  if ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0)
645  modified_points[i][d] -= periodicity[d];
646  surrounding_points_start = modified_points.data();
647  }
648 
649  // Now perform the interpolation
650  for (unsigned int row = 0; row < weights.size(0); ++row)
651  {
652  Assert(
653  std::abs(
654  std::accumulate(&weights(row, 0), &weights(row, 0) + n_points, 0.0) -
655  1.0) < 1e-10,
656  ExcMessage("The weights for the individual points should sum to 1!"));
657  Point<spacedim> new_point;
658  for (unsigned int p = 0; p < n_points; ++p)
659  new_point += surrounding_points_start[p] * weights(row, p);
660 
661  // if necessary, also adjust the weighted point by the periodicity
662  for (unsigned int d = 0; d < spacedim; ++d)
663  if (periodicity[d] > 0)
664  if (new_point[d] < 0)
665  new_point[d] += periodicity[d];
666 
667  // TODO should this use surrounding_points_start or surrounding_points?
668  // The older version used surrounding_points
669  new_points[row] =
670  project_to_manifold(make_array_view(surrounding_points.begin(),
671  surrounding_points.end()),
672  new_point);
673  }
674 }
675 
676 
677 
678 template <int dim, int spacedim>
681  const ArrayView<const Point<spacedim>> & /*vertices*/,
682  const Point<spacedim> &candidate) const
683 {
684  return candidate;
685 }
686 
687 
688 
689 template <int dim, int spacedim>
690 const Tensor<1, spacedim> &
692 {
693  return periodicity;
694 }
695 
696 
697 
698 template <int dim, int spacedim>
701  const Point<spacedim> &x2) const
702 {
703  Tensor<1, spacedim> direction = x2 - x1;
704 
705  // see if we have to take into account periodicity. if so, we need
706  // to make sure that if a distance in one coordinate direction
707  // is larger than half of the box length, then go the other way
708  // around (i.e., via the periodic box)
709  for (unsigned int d = 0; d < spacedim; ++d)
710  if (periodicity[d] > tolerance)
711  {
712  if (direction[d] < -periodicity[d] / 2)
713  direction[d] += periodicity[d];
714  else if (direction[d] > periodicity[d] / 2)
715  direction[d] -= periodicity[d];
716  }
717 
718  return direction;
719 }
720 
721 
722 
723 template <>
724 void
728 {
729  Assert(false, ExcImpossibleInDim(1));
730 }
731 
732 
733 
734 template <>
735 void
739 {
740  Assert(false, ExcNotImplemented());
741 }
742 
743 
744 
745 template <>
746 void
750 {
751  Assert(false, ExcNotImplemented());
752 }
753 
754 
755 
756 template <>
757 void
760  Manifold<2, 2>::FaceVertexNormals & face_vertex_normals) const
761 {
762  const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
763  for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_face;
764  ++vertex)
765  // compute normals from tangent
766  face_vertex_normals[vertex] = Point<2>(tangent[1], -tangent[0]);
767 }
768 
769 
770 
771 template <>
772 void
774  const Triangulation<2, 3>::face_iterator & /*face*/,
775  Manifold<2, 3>::FaceVertexNormals & /*face_vertex_normals*/) const
776 {
777  Assert(false, ExcNotImplemented());
778 }
779 
780 
781 
782 template <>
783 void
786  Manifold<3, 3>::FaceVertexNormals & face_vertex_normals) const
787 {
788  const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
789 
790  static const unsigned int neighboring_vertices[4][2] = {{1, 2},
791  {3, 0},
792  {0, 3},
793  {2, 1}};
794  for (unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
795  {
796  // first define the two tangent vectors at the vertex by using the
797  // two lines radiating away from this vertex
798  const Tensor<1, 3> tangents[2] = {
799  face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
800  face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
801 
802  // then compute the normal by taking the cross product. since the
803  // normal is not required to be normalized, no problem here
804  face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
805  }
806 }
807 
808 
809 
810 template <>
813  const Point<1> &) const
814 {
815  Assert(false, ExcNotImplemented());
816  return {};
817 }
818 
819 
820 
821 template <>
824  const Point<2> &) const
825 {
826  Assert(false, ExcNotImplemented());
827  return {};
828 }
829 
830 
831 
832 template <>
835  const Point<3> &) const
836 {
837  Assert(false, ExcNotImplemented());
838  return {};
839 }
840 
841 
842 
843 template <>
847  const Point<2> & p) const
848 {
849  // In 2d, a face is just a straight line and
850  // we can use the 'standard' implementation.
851  return Manifold<2, 2>::normal_vector(face, p);
852 }
853 
854 
855 
856 template <int dim, int spacedim>
859  const typename Triangulation<dim, spacedim>::face_iterator &face,
860  const Point<spacedim> & p) const
861 {
862  // I don't think the implementation below will work when dim!=spacedim;
863  // in fact, I believe that we don't even have enough information here,
864  // because we would need to know not only about the tangent vectors
865  // of the face, but also of the cell, to compute the normal vector.
866  // Someone will have to think about this some more.
867  Assert(dim == spacedim, ExcNotImplemented());
868 
869  // in order to find out what the normal vector is, we first need to
870  // find the reference coordinates of the point p on the given face,
871  // or at least the reference coordinates of the closest point on the
872  // face
873  //
874  // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
875  // where F(xi) is the mapping. this algorithm is implemented in
876  // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
877  // while we need it for faces here. it's also implemented in somewhat
878  // more generality there using the machinery of the MappingQ1 class
879  // while we really only need it for a specific case here
880  //
881  // in any case, the iteration we use here is a Gauss-Newton's iteration with
882  // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
883  // where
884  // J(xi) = (grad F(xi))^T (F(xi)-p)
885  // and
886  // H(xi) = [grad F(xi)]^T [grad F(xi)]
887  // In all this,
888  // F(xi) = sum_v vertex[v] phi_v(xi)
889  // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
890 
891  // we start with the point xi=1/2, xi=(1/2,1/2), ...
892  const unsigned int facedim = dim - 1;
893 
894  Point<facedim> xi;
895  for (unsigned int i = 0; i < facedim; ++i)
896  xi[i] = 1. / 2;
897 
898  const double eps = 1e-12;
899  Tensor<1, spacedim> grad_F[facedim];
900  unsigned int iteration = 0;
901  while (true)
902  {
904  for (const unsigned int v : GeometryInfo<facedim>::vertex_indices())
905  F += face->vertex(v) *
907 
908  for (unsigned int i = 0; i < facedim; ++i)
909  {
910  grad_F[i] = 0;
911  for (const unsigned int v : GeometryInfo<facedim>::vertex_indices())
912  grad_F[i] +=
913  face->vertex(v) *
915  }
916 
918  for (unsigned int i = 0; i < facedim; ++i)
919  for (unsigned int j = 0; j < spacedim; ++j)
920  J[i] += grad_F[i][j] * (F - p)[j];
921 
923  for (unsigned int i = 0; i < facedim; ++i)
924  for (unsigned int j = 0; j < facedim; ++j)
925  for (unsigned int k = 0; k < spacedim; ++k)
926  H[i][j] += grad_F[i][k] * grad_F[j][k];
927 
928  const Tensor<1, facedim> delta_xi = -invert(H) * J;
929  xi += delta_xi;
930  ++iteration;
931 
932  Assert(iteration < 10,
933  ExcMessage("The Newton iteration to find the reference point "
934  "did not converge in 10 iterations. Do you have a "
935  "deformed cell? (See the glossary for a definition "
936  "of what a deformed cell is. You may want to output "
937  "the vertices of your cell."));
938 
939  // It turns out that the check in reference coordinates with an absolute
940  // tolerance can cause a convergence failure of the Newton method as
941  // seen in tests/manifold/flat_manifold_09.cc. To work around this, also
942  // use a convergence check in world coordinates. This check has to be
943  // relative to the size of the face of course. Here we decided to use
944  // diameter because it works for non-planar faces and is cheap to
945  // compute:
946  const double normalized_delta_world = (F - p).norm() / face->diameter();
947 
948  if (delta_xi.norm() < eps || normalized_delta_world < eps)
949  break;
950  }
951 
952  // so now we have the reference coordinates xi of the point p.
953  // we then have to compute the normal vector, which we can do
954  // by taking the (normalize) alternating product of all the tangent
955  // vectors given by grad_F
956  return internal::normalized_alternating_product(grad_F);
957 }
958 
959 
960 /* -------------------------- ChartManifold --------------------- */
961 template <int dim, int spacedim, int chartdim>
964  : sub_manifold(periodicity)
965 {}
966 
967 
968 
969 template <int dim, int spacedim, int chartdim>
972  const Point<spacedim> &p1,
973  const Point<spacedim> &p2,
974  const double w) const
975 {
976  const std::array<Point<spacedim>, 2> points{{p1, p2}};
977  const std::array<double, 2> weights{{1. - w, w}};
978  return get_new_point(make_array_view(points.begin(), points.end()),
979  make_array_view(weights.begin(), weights.end()));
980 }
981 
982 
983 
984 template <int dim, int spacedim, int chartdim>
987  const ArrayView<const Point<spacedim>> &surrounding_points,
988  const ArrayView<const double> & weights) const
989 {
990  const std::size_t n_points = surrounding_points.size();
991 
992  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
993 
994  for (unsigned int i = 0; i < n_points; ++i)
995  chart_points[i] = pull_back(surrounding_points[i]);
996 
998  make_array_view(chart_points.begin(), chart_points.end()), weights);
999 
1000  return push_forward(p_chart);
1001 }
1002 
1003 
1004 
1005 template <int dim, int spacedim, int chartdim>
1006 void
1008  const ArrayView<const Point<spacedim>> &surrounding_points,
1009  const Table<2, double> & weights,
1010  ArrayView<Point<spacedim>> new_points) const
1011 {
1012  Assert(weights.size(0) > 0, ExcEmptyObject());
1013  AssertDimension(surrounding_points.size(), weights.size(1));
1014 
1015  const std::size_t n_points = surrounding_points.size();
1016 
1017  boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1018  for (std::size_t i = 0; i < n_points; ++i)
1019  chart_points[i] = pull_back(surrounding_points[i]);
1020 
1021  boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1022  weights.size(0));
1024  make_array_view(chart_points.begin(), chart_points.end()),
1025  weights,
1026  make_array_view(new_points_on_chart.begin(), new_points_on_chart.end()));
1027 
1028  for (std::size_t row = 0; row < weights.size(0); ++row)
1029  new_points[row] = push_forward(new_points_on_chart[row]);
1030 }
1031 
1032 
1033 
1034 template <int dim, int spacedim, int chartdim>
1037  const Point<chartdim> &) const
1038 {
1039  // function must be implemented in a derived class to be usable,
1040  // as discussed in this function's documentation
1041  Assert(false, ExcPureFunctionCalled());
1043 }
1044 
1045 
1046 
1047 template <int dim, int spacedim, int chartdim>
1050  const Point<spacedim> &x1,
1051  const Point<spacedim> &x2) const
1052 {
1055 
1056  // ensure that the chart is not singular by asserting that its
1057  // derivative has a positive determinant. we need to make this
1058  // comparison relative to the size of the derivative. since the
1059  // determinant is the product of chartdim factors, take the
1060  // chartdim-th root of it in comparing against the size of the
1061  // derivative
1062  Assert(std::pow(std::abs(F_prime.determinant()), 1. / chartdim) >=
1063  1e-12 * F_prime.norm(),
1064  ExcMessage(
1065  "The derivative of a chart function must not be singular."));
1066 
1067  const Tensor<1, chartdim> delta =
1069 
1070  Tensor<1, spacedim> result;
1071  for (unsigned int i = 0; i < spacedim; ++i)
1072  result[i] += F_prime[i] * delta;
1073 
1074  return result;
1075 }
1076 
1077 
1078 
1079 template <int dim, int spacedim, int chartdim>
1080 const Tensor<1, chartdim> &
1082 {
1083  return sub_manifold.get_periodicity();
1084 }
1085 
1086 // explicit instantiations
1087 #include "manifold.inst"
1088 
static ::ExceptionBase & ExcPureFunctionCalled()
static const unsigned int invalid_unsigned_int
Definition: types.h:196
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:858
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1623
typename IteratorSelector::line_iterator line_iterator
Definition: tria.h:1444
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const Tensor< 1, spacedim > & get_periodicity() const
Definition: manifold.cc:691
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:1492
iterator end() const
Definition: array_view.h:560
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:472
static Tensor< 1, dim > d_linear_shape_function_gradient(const Point< dim > &xi, const unsigned int i)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Number determinant() const
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
Definition: tensor.h:2441
std::size_t size() const
Definition: array_view.h:542
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)
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:1036
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:1466
static ::ExceptionBase & ExcPeriodicBox(int arg1, Point< spacedim > arg2, double arg3)
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:665
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
const Tensor< 1, chartdim > & get_periodicity() const
Definition: manifold.cc:1081
Point< 3 > vertices[4]
virtual Point< spacedim > project_to_manifold(const ArrayView< const Point< spacedim >> &points, const Point< spacedim > &candidate) const override
Definition: manifold.cc:680
numbers::NumberTraits< Number >::real_type norm() const
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
virtual std::unique_ptr< Manifold< dim, spacedim > > clone() const override
Definition: manifold.cc:528
typename IteratorSelector::quad_iterator quad_iterator
Definition: tria.h:1468
const double tolerance
Definition: manifold.h:808
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:1085
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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:1007
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:1049
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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:448
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition: manifold.cc:700
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim >> &surrounding_points, const ArrayView< const double > &weights) const override
Definition: manifold.cc:986
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
T min(const T &t, const MPI_Comm &mpi_communicator)
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:537
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:307
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:551
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
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
Definition: manifold.cc:598
constexpr Tensor< 1, dim, Number > cross_product_2d(const Tensor< 1, dim, Number > &src)
Definition: tensor.h:2416
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition: manifold.cc:962
Point< spacedim > get_new_point_on_face(const typename Triangulation< dim, spacedim >::face_iterator &face) const
Definition: manifold.cc:344
void copy(const T *begin, const T *end, U *dest)
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
T max(const T &t, const MPI_Comm &mpi_communicator)
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition: manifold.cc:971
const Tensor< 1, spacedim > periodicity
Definition: manifold.h:794
static ::ExceptionBase & ExcInternalError()
Point< spacedim > get_new_point_on_cell(const typename Triangulation< dim, spacedim >::cell_iterator &cell) const
Definition: manifold.cc:364