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