deal.II version GIT relicensing-1936-ge55124254b 2024-10-04 21:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
manifold.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2014 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/table.h>
16#include <deal.II/base/tensor.h>
18
19#include <deal.II/fe/fe_q.h>
20
23#include <deal.II/grid/tria.h>
26
27#include <boost/container/small_vector.hpp>
28
29#include <cmath>
30#include <limits>
31#include <memory>
32
34
35/* -------------------------- Manifold --------------------- */
36#ifndef DOXYGEN
37
38template <int dim, int spacedim>
41 const ArrayView<const Point<spacedim>> &,
42 const Point<spacedim> &) const
43{
45 return {};
46}
47
48
49
50template <int dim, int spacedim>
53 const Point<spacedim> &p2,
54 const double w) const
55{
56 const std::array<Point<spacedim>, 2> vertices{{p1, p2}};
57 return project_to_manifold(make_array_view(vertices), 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(surrounding_points, make_array_view(weights, row));
135 }
136}
137
138
139
140template <>
143 const Point<2> &p) const
144{
145 const int spacedim = 2;
146
147 // get the tangent vector from the point 'p' in the direction of the further
148 // one of the two vertices that make up the face of this 2d cell
149 const Tensor<1, spacedim> tangent =
150 ((p - face->vertex(0)).norm_square() > (p - face->vertex(1)).norm_square() ?
151 -get_tangent_vector(p, face->vertex(0)) :
152 get_tangent_vector(p, face->vertex(1)));
153
154 // then rotate it by 90 degrees
155 const Tensor<1, spacedim> normal = cross_product_2d(tangent);
156 return normal / normal.norm();
157}
158
159
160
161template <>
164 const Point<3> &p) const
165{
166 const int spacedim = 3;
167
168 const std::array<Point<spacedim>, 4> vertices{
169 {face->vertex(0), face->vertex(1), face->vertex(2), face->vertex(3)}};
170 const std::array<double, 4> distances{{vertices[0].distance(p),
171 vertices[1].distance(p),
172 vertices[2].distance(p),
173 vertices[3].distance(p)}};
174 const double max_distance = std::max(std::max(distances[0], distances[1]),
175 std::max(distances[2], distances[3]));
176
177 // We need to find two tangential vectors to the given point p, but we do
178 // not know how the point is oriented against the face. We guess the two
179 // directions by assuming a flat topology and take the two directions that
180 // indicate the angle closest to a perpendicular one (i.e., cos(theta) close
181 // to zero). We start with an invalid value but the loops should always find
182 // a value.
183 double abs_cos_angle = std::numeric_limits<double>::max();
184 unsigned int first_index = numbers::invalid_unsigned_int,
185 second_index = numbers::invalid_unsigned_int;
186 for (unsigned int i = 0; i < 3; ++i)
187 if (distances[i] > 1e-8 * max_distance)
188 for (unsigned int j = i + 1; j < 4; ++j)
189 if (distances[j] > 1e-8 * max_distance)
190 {
191 const double new_angle = (p - vertices[i]) * (p - vertices[j]) /
192 (distances[i] * distances[j]);
193 // multiply by factor 0.999 to bias the search in a way that
194 // avoids trouble with roundoff
195 if (std::abs(new_angle) < 0.999 * abs_cos_angle)
196 {
197 abs_cos_angle = std::abs(new_angle);
198 first_index = i;
199 second_index = j;
200 }
201 }
203 ExcMessage("The search for possible directions did not succeed."));
204
205 // Compute tangents and normal for selected vertices
208 Tensor<1, spacedim> normal;
209
210 bool done = false;
211 std::vector<bool> tested_vertices(vertices.size(), false);
212 tested_vertices[first_index] = true;
213 tested_vertices[second_index] = true;
214
215 do
216 {
217 // Compute tangents and normal for selected vertices
218 t1 = get_tangent_vector(p, vertices[first_index]);
219 t2 = get_tangent_vector(p, vertices[second_index]);
220 normal = cross_product_3d(t1, t2);
221
222 // if normal is zero, try some other combination of vertices
223 if (normal.norm_square() < 1e4 * std::numeric_limits<double>::epsilon() *
224 t1.norm_square() * t2.norm_square())
225 {
226 // See if we have tested all vertices already
227 auto first_false =
228 std::find(tested_vertices.begin(), tested_vertices.end(), false);
229 if (first_false == tested_vertices.end())
230 {
231 done = true;
232 }
233 else
234 {
235 *first_false = true;
236 second_index = first_false - tested_vertices.begin();
237 }
238 }
239 else
240 {
241 done = true;
242 }
243 }
244 while (!done);
245
246 Assert(
247 normal.norm_square() > 1e4 * std::numeric_limits<double>::epsilon() *
248 t1.norm_square() * t2.norm_square(),
250 "Manifold::normal_vector was unable to find a suitable combination "
251 "of vertices to compute a normal on this face. We chose the secants "
252 "that are as orthogonal as possible, but tangents appear to be "
253 "linearly dependent. Check for distorted faces in your triangulation."));
254
255 // Now figure out if we need to flip the direction, we do this by comparing
256 // to a reference normal that would be the correct result if all vertices
257 // would lie in a plane
258 const Tensor<1, spacedim> rt1 = vertices[3] - vertices[0];
259 const Tensor<1, spacedim> rt2 = vertices[2] - vertices[1];
260 const Tensor<1, spacedim> reference_normal = cross_product_3d(rt1, rt2);
261
262 if (reference_normal * normal < 0.0)
263 normal *= -1.0;
264
265 return normal / normal.norm();
266}
267
268
269
270template <int dim, int spacedim>
273 const typename Triangulation<dim, spacedim>::face_iterator & /*face*/,
274 const Point<spacedim> & /*p*/) const
275{
277 return Tensor<1, spacedim>();
278}
279
280
281
282template <>
283void
286 FaceVertexNormals &n) const
287{
288 n[0] = cross_product_2d(get_tangent_vector(face->vertex(0), face->vertex(1)));
289 n[1] =
290 -cross_product_2d(get_tangent_vector(face->vertex(1), face->vertex(0)));
291
292 for (unsigned int i = 0; i < 2; ++i)
293 {
294 Assert(n[i].norm() != 0,
295 ExcInternalError("Something went wrong. The "
296 "computed normals have "
297 "zero length."));
298 n[i] /= n[i].norm();
299 }
300}
301
302
303
304template <>
305void
308 FaceVertexNormals &n) const
309{
310 n[0] = cross_product_3d(get_tangent_vector(face->vertex(0), face->vertex(1)),
311 get_tangent_vector(face->vertex(0), face->vertex(2)));
312
313 n[1] = cross_product_3d(get_tangent_vector(face->vertex(1), face->vertex(3)),
314 get_tangent_vector(face->vertex(1), face->vertex(0)));
315
316 n[2] = cross_product_3d(get_tangent_vector(face->vertex(2), face->vertex(0)),
317 get_tangent_vector(face->vertex(2), face->vertex(3)));
318
319 n[3] = cross_product_3d(get_tangent_vector(face->vertex(3), face->vertex(2)),
320 get_tangent_vector(face->vertex(3), face->vertex(1)));
321
322 for (unsigned int i = 0; i < 4; ++i)
323 {
324 Assert(n[i].norm() != 0,
325 ExcInternalError("Something went wrong. The "
326 "computed normals have "
327 "zero length."));
328 n[i] /= n[i].norm();
329 }
330}
331
332
333
334template <int dim, int spacedim>
335void
338 FaceVertexNormals &n) const
339{
340 for (unsigned int v = 0; v < face->reference_cell().n_vertices(); ++v)
341 {
342 n[v] = normal_vector(face, face->vertex(v));
343 n[v] /= n[v].norm();
344 }
345}
346
347
348
349template <int dim, int spacedim>
352 const typename Triangulation<dim, spacedim>::line_iterator &line) const
353{
354 const auto points_weights = Manifolds::get_default_points_and_weights(line);
355 return get_new_point(make_array_view(points_weights.first),
356 make_array_view(points_weights.second));
357}
358
359
360
361template <int dim, int spacedim>
364 const typename Triangulation<dim, spacedim>::quad_iterator &quad) const
365{
366 const auto points_weights = Manifolds::get_default_points_and_weights(quad);
367 return get_new_point(make_array_view(points_weights.first),
368 make_array_view(points_weights.second));
369}
370
371
372
373template <int dim, int spacedim>
376 const typename Triangulation<dim, spacedim>::face_iterator &face) const
377{
378 Assert(dim > 1, ExcImpossibleInDim(dim));
379
380 switch (dim)
381 {
382 case 2:
383 return get_new_point_on_line(face);
384 case 3:
385 return get_new_point_on_quad(face);
386 }
387
388 return {};
389}
390
391
392
393template <int dim, int spacedim>
396 const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
397{
398 switch (dim)
399 {
400 case 1:
401 return get_new_point_on_line(cell);
402 case 2:
403 return get_new_point_on_quad(cell);
404 case 3:
405 return get_new_point_on_hex(cell);
406 }
407
408 return {};
409}
410
411
412
413template <>
417{
418 Assert(false, ExcImpossibleInDim(1));
419 return {};
420}
421
422
423
424template <>
428{
429 Assert(false, ExcImpossibleInDim(1));
430 return {};
431}
432
433
434
435template <>
439{
440 Assert(false, ExcImpossibleInDim(1));
441 return {};
442}
443
444
445
446template <>
450{
451 Assert(false, ExcImpossibleInDim(1));
452 return {};
453}
454
455
456
457template <>
461{
462 Assert(false, ExcImpossibleInDim(1));
463 return {};
464}
465
466
467
468template <>
472{
473 Assert(false, ExcImpossibleInDim(1));
474 return {};
475}
476
477
478
479template <int dim, int spacedim>
482 const typename Triangulation<dim, spacedim>::hex_iterator & /*hex*/) const
483{
484 Assert(false, ExcImpossibleInDim(dim));
485 return {};
486}
487
488
489
490template <>
493 const Triangulation<3, 3>::hex_iterator &hex) const
494{
495 const auto points_weights =
497 return get_new_point(make_array_view(points_weights.first),
498 make_array_view(points_weights.second));
499}
500
501
502
503template <int dim, int spacedim>
506 const Point<spacedim> &x2) const
507{
508 const double epsilon = 1e-8;
509
510 const std::array<Point<spacedim>, 2> points{{x1, x2}};
511 const std::array<double, 2> weights{{epsilon, 1.0 - epsilon}};
512 const Point<spacedim> neighbor_point =
513 get_new_point(make_array_view(points), make_array_view(weights));
514 return (neighbor_point - x1) / epsilon;
515}
516#endif
517/* -------------------------- FlatManifold --------------------- */
518
519namespace internal
520{
521 namespace
522 {
524 normalized_alternating_product(const Tensor<1, 3> (&)[1])
525 {
526 // we get here from FlatManifold<2,3>::normal_vector, but
527 // the implementation below is bogus for this case anyway
528 // (see the assert at the beginning of that function).
530 return {};
531 }
532
533
534
536 normalized_alternating_product(const Tensor<1, 3> (&basis_vectors)[2])
537 {
538 Tensor<1, 3> tmp = cross_product_3d(basis_vectors[0], basis_vectors[1]);
539 return tmp / tmp.norm();
540 }
541
542 } // namespace
543} // namespace internal
544
545#ifndef DOXYGEN
546
547template <int dim, int spacedim>
549 const Tensor<1, spacedim> &periodicity,
550 const double tolerance)
551 : periodicity(periodicity)
552 , tolerance(tolerance)
553{}
554
555
556
557template <int dim, int spacedim>
558std::unique_ptr<Manifold<dim, spacedim>>
560{
561 return std::make_unique<FlatManifold<dim, spacedim>>(periodicity, tolerance);
562}
563
564
565
566template <int dim, int spacedim>
569 const ArrayView<const Point<spacedim>> &surrounding_points,
570 const ArrayView<const double> &weights) const
571{
572 Assert(std::abs(std::accumulate(weights.begin(), weights.end(), 0.0) - 1.0) <
573 1e-10,
574 ExcMessage("The weights for the new point should sum to 1!"));
575
577
578 // if there is no periodicity, use a shortcut
579 if (periodicity == Tensor<1, spacedim>())
580 {
581 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
582 p += surrounding_points[i] * weights[i];
583 }
584 else
585 {
586 Tensor<1, spacedim> minP = periodicity;
587
588 for (unsigned int d = 0; d < spacedim; ++d)
589 if (periodicity[d] > 0)
590 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
591 {
592 minP[d] = std::min(minP[d], surrounding_points[i][d]);
593 Assert((surrounding_points[i][d] <
594 periodicity[d] + tolerance * periodicity[d]) ||
595 (surrounding_points[i][d] >=
596 -tolerance * periodicity[d]),
597 ExcPeriodicBox(d, surrounding_points[i], periodicity[d]));
598 }
599
600 // compute the weighted average point, possibly taking into account
601 // periodicity
602 for (unsigned int i = 0; i < surrounding_points.size(); ++i)
603 {
605 for (unsigned int d = 0; d < spacedim; ++d)
606 if (periodicity[d] > 0)
607 dp[d] =
608 ((surrounding_points[i][d] - minP[d]) > periodicity[d] / 2.0 ?
609 -periodicity[d] :
610 0.0);
611
612 p += (surrounding_points[i] + dp) * weights[i];
613 }
614
615 // if necessary, also adjust the weighted point by the periodicity
616 for (unsigned int d = 0; d < spacedim; ++d)
617 if (periodicity[d] > 0)
618 if (p[d] < 0)
619 p[d] += periodicity[d];
620 }
621
622 return project_to_manifold(surrounding_points, p);
623}
624
625
626
627template <int dim, int spacedim>
628void
630 const ArrayView<const Point<spacedim>> &surrounding_points,
631 const Table<2, double> &weights,
632 ArrayView<Point<spacedim>> new_points) const
633{
634 AssertDimension(surrounding_points.size(), weights.size(1));
635 if (weights.size(0) == 0)
636 return;
637 AssertDimension(new_points.size(), weights.size(0));
638
639 const std::size_t n_points = surrounding_points.size();
640
641 // if there is no periodicity, use an optimized implementation with
642 // VectorizedArray, otherwise go to the get_new_point function for adjusting
643 // the domain
644 if (periodicity == Tensor<1, spacedim>())
645 {
646 for (unsigned int row = 0; row < weights.size(0); ++row)
647 Assert(std::abs(std::accumulate(&weights(row, 0),
648 &weights(row, 0) + n_points,
649 0.0) -
650 1.0) < 1e-10,
651 ExcMessage("The weights for each of the points should sum to "
652 "1!"));
653
654 constexpr std::size_t n_lanes =
655 std::min<std::size_t>(VectorizedArray<double>::size(), 4);
656 using VectorizedArrayType = VectorizedArray<double, n_lanes>;
657 const std::size_t n_regular_cols = (n_points / n_lanes) * n_lanes;
658 for (unsigned int row = 0; row < weights.size(0); row += n_lanes)
659 {
660 std::array<unsigned int, n_lanes> offsets;
661 // ensure to not access out of bounds, possibly duplicating some
662 // entries
663 for (std::size_t i = 0; i < n_lanes; ++i)
664 offsets[i] =
665 std::min<unsigned int>((row + i) * n_points,
666 (weights.size(0) - 1) * n_points);
668 for (std::size_t col = 0; col < n_regular_cols; col += n_lanes)
669 {
670 std::array<VectorizedArrayType, n_lanes> vectorized_weights;
672 &weights(0, 0) + col,
673 offsets.data(),
674 vectorized_weights.data());
675 for (std::size_t i = 0; i < n_lanes; ++i)
676 point += vectorized_weights[i] * surrounding_points[col + i];
677 }
678 for (std::size_t col = n_regular_cols; col < n_points; ++col)
679 {
680 VectorizedArrayType vectorized_weights;
681 vectorized_weights.gather(&weights(0, 0) + col, offsets.data());
682 point += vectorized_weights * surrounding_points[col];
683 }
684 for (unsigned int r = row;
685 r < std::min<unsigned int>(weights.size(0), row + n_lanes);
686 ++r)
687 {
688 // unpack and project to manifold
689 for (unsigned int d = 0; d < spacedim; ++d)
690 new_points[r][d] = point[d][r - row];
691 new_points[r] =
692 project_to_manifold(surrounding_points, new_points[r]);
693 }
694 }
695 }
696 else
697 for (unsigned int row = 0; row < weights.size(0); ++row)
698 new_points[row] =
699 get_new_point(surrounding_points,
700 ArrayView<const double>(&weights(row, 0), n_points));
701}
702
703
704
705template <int dim, int spacedim>
708 const ArrayView<const Point<spacedim>> & /*vertices*/,
709 const Point<spacedim> &candidate) const
710{
711 return candidate;
712}
713
714
715
716template <int dim, int spacedim>
719{
720 return periodicity;
721}
722
723
724
725template <int dim, int spacedim>
728 const Point<spacedim> &x2) const
729{
730 Tensor<1, spacedim> direction = x2 - x1;
731
732 // see if we have to take into account periodicity. if so, we need
733 // to make sure that if a distance in one coordinate direction
734 // is larger than half of the box length, then go the other way
735 // around (i.e., via the periodic box)
736 for (unsigned int d = 0; d < spacedim; ++d)
737 if (periodicity[d] > tolerance)
738 {
739 if (direction[d] < -periodicity[d] / 2)
740 direction[d] += periodicity[d];
741 else if (direction[d] > periodicity[d] / 2)
742 direction[d] -= periodicity[d];
743 }
744
745 return direction;
746}
747
748
749
750template <>
751void
755{
756 Assert(false, ExcImpossibleInDim(1));
757}
758
759
760
761template <>
762void
766{
768}
769
770
771
772template <>
773void
777{
779}
780
781
782
783template <>
784void
787 Manifold<2, 2>::FaceVertexNormals &face_vertex_normals) const
788{
789 const Tensor<1, 2> tangent = face->vertex(1) - face->vertex(0);
790 // We're in 2d. Faces are edges:
791 for (const unsigned int vertex : ReferenceCells::Line.vertex_indices())
792 // compute normals from tangent
793 face_vertex_normals[vertex] = Tensor<1, 2>({tangent[1], -tangent[0]});
794}
795
796
797
798template <>
799void
801 const Triangulation<2, 3>::face_iterator & /*face*/,
802 Manifold<2, 3>::FaceVertexNormals & /*face_vertex_normals*/) const
803{
805}
806
807
808
809template <>
810void
813 Manifold<3, 3>::FaceVertexNormals &face_vertex_normals) const
814{
815 const unsigned int vertices_per_face = GeometryInfo<3>::vertices_per_face;
816
817 static const unsigned int neighboring_vertices[4][2] = {{1, 2},
818 {3, 0},
819 {0, 3},
820 {2, 1}};
821 for (unsigned int vertex = 0; vertex < vertices_per_face; ++vertex)
822 {
823 // first define the two tangent vectors at the vertex by using the
824 // two lines radiating away from this vertex
825 const Tensor<1, 3> tangents[2] = {
826 face->vertex(neighboring_vertices[vertex][0]) - face->vertex(vertex),
827 face->vertex(neighboring_vertices[vertex][1]) - face->vertex(vertex)};
828
829 // then compute the normal by taking the cross product. since the
830 // normal is not required to be normalized, no problem here
831 face_vertex_normals[vertex] = cross_product_3d(tangents[0], tangents[1]);
832 }
833}
834
835
836
837template <>
840 const Point<1> &) const
841{
843 return {};
844}
845
846
847
848template <>
851 const Point<2> &) const
852{
854 return {};
855}
856
857
858
859template <>
862 const Point<3> &) const
863{
865 return {};
866}
867
868
869
870template <>
874 const Point<2> &p) const
875{
876 // In 2d, a face is just a straight line and
877 // we can use the 'standard' implementation.
878 return Manifold<2, 2>::normal_vector(face, p);
879}
880
881
882
883template <int dim, int spacedim>
887 const Point<spacedim> &p) const
888{
889 // I don't think the implementation below will work when dim!=spacedim;
890 // in fact, I believe that we don't even have enough information here,
891 // because we would need to know not only about the tangent vectors
892 // of the face, but also of the cell, to compute the normal vector.
893 // Someone will have to think about this some more.
894 Assert(dim == spacedim, ExcNotImplemented());
895
896 // in order to find out what the normal vector is, we first need to
897 // find the reference coordinates of the point p on the given face,
898 // or at least the reference coordinates of the closest point on the
899 // face
900 //
901 // in other words, we need to find a point xi so that f(xi)=||F(xi)-p||^2->min
902 // where F(xi) is the mapping. this algorithm is implemented in
903 // MappingQ1<dim,spacedim>::transform_real_to_unit_cell but only for cells,
904 // while we need it for faces here. it's also implemented in somewhat
905 // more generality there using the machinery of the MappingQ1 class
906 // while we really only need it for a specific case here
907 //
908 // in any case, the iteration we use here is a Gauss-Newton's iteration with
909 // xi^{n+1} = xi^n - H(xi^n)^{-1} J(xi^n)
910 // where
911 // J(xi) = (grad F(xi))^T (F(xi)-p)
912 // and
913 // H(xi) = [grad F(xi)]^T [grad F(xi)]
914 // In all this,
915 // F(xi) = sum_v vertex[v] phi_v(xi)
916 // We get the shape functions phi_v from an object of type FE_Q<dim-1>(1)
917
918 // We start at the center of the cell. If the face is a line or
919 // square, then the center is at 0.5 or (0.5,0.5). If the face is
920 // a triangle, then we start at the point (1/3,1/3).
921 const unsigned int facedim = dim - 1;
922
924
925 const auto face_reference_cell = face->reference_cell();
926
927 if ((dim <= 2) ||
928 (face_reference_cell == ReferenceCells::get_hypercube<facedim>()))
930 for (unsigned int i = 0; i < facedim; ++i)
931 xi[i] = 1. / 2;
932 }
933 else
934 {
935 for (unsigned int i = 0; i < facedim; ++i)
936 xi[i] = 1. / 3;
937 }
938
939 const double eps = 1e-12;
940 Tensor<1, spacedim> grad_F[facedim];
941 unsigned int iteration = 0;
942 while (true)
943 {
945 for (const unsigned int v : face->vertex_indices())
946 F +=
947 face->vertex(v) * face_reference_cell.d_linear_shape_function(xi, v);
948
949 for (unsigned int i = 0; i < facedim; ++i)
950 {
951 grad_F[i] = 0;
952 for (const unsigned int v : face->vertex_indices())
953 grad_F[i] +=
954 face->vertex(v) *
955 face_reference_cell.d_linear_shape_function_gradient(xi, v)[i];
956 }
957
959 for (unsigned int i = 0; i < facedim; ++i)
960 for (unsigned int j = 0; j < spacedim; ++j)
961 J[i] += grad_F[i][j] * (F - p)[j];
962
964 for (unsigned int i = 0; i < facedim; ++i)
965 for (unsigned int j = 0; j < facedim; ++j)
966 for (unsigned int k = 0; k < spacedim; ++k)
967 H[i][j] += grad_F[i][k] * grad_F[j][k];
968
969 const Tensor<1, facedim> delta_xi = -invert(H) * J;
970 xi += delta_xi;
971 ++iteration;
972
973 AssertThrow(iteration < 10,
975 "The Newton iteration to find the reference point "
976 "did not converge in 10 iterations. Do you have a "
977 "deformed cell? (See the glossary for a definition "
978 "of what a deformed cell is. You may want to output "
979 "the vertices of your cell."));
980
981 // It turns out that the check in reference coordinates with an absolute
982 // tolerance can cause a convergence failure of the Newton method as
983 // seen in tests/manifold/flat_manifold_09.cc. To work around this, also
984 // use a convergence check in world coordinates. This check has to be
985 // relative to the size of the face of course. Here we decided to use
986 // diameter because it works for non-planar faces and is cheap to
987 // compute:
988 const double normalized_delta_world = (F - p).norm() / face->diameter();
989
990 if (delta_xi.norm() < eps || normalized_delta_world < eps)
991 break;
992 }
993
994 // so now we have the reference coordinates xi of the point p.
995 // we then have to compute the normal vector, which we can do
996 // by taking the (normalize) alternating product of all the tangent
997 // vectors given by grad_F
998 return internal::normalized_alternating_product(grad_F);
999}
1000#endif
1001
1002/* -------------------------- ChartManifold --------------------- */
1003template <int dim, int spacedim, int chartdim>
1005 const Tensor<1, chartdim> &periodicity)
1006 : sub_manifold(periodicity)
1007{}
1008
1009
1010
1011template <int dim, int spacedim, int chartdim>
1015 const Point<spacedim> &p2,
1016 const double w) const
1017{
1018 const std::array<Point<spacedim>, 2> points{{p1, p2}};
1019 const std::array<double, 2> weights{{1. - w, w}};
1020 return get_new_point(make_array_view(points), make_array_view(weights));
1021}
1022
1023
1024
1025template <int dim, int spacedim, int chartdim>
1028 const ArrayView<const Point<spacedim>> &surrounding_points,
1029 const ArrayView<const double> &weights) const
1030{
1031 const std::size_t n_points = surrounding_points.size();
1032
1033 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1034
1035 for (unsigned int i = 0; i < n_points; ++i)
1036 chart_points[i] = pull_back(surrounding_points[i]);
1037
1038 const Point<chartdim> p_chart =
1039 sub_manifold.get_new_point(chart_points, weights);
1040
1041 return push_forward(p_chart);
1042}
1043
1044
1045
1046template <int dim, int spacedim, int chartdim>
1047void
1049 const ArrayView<const Point<spacedim>> &surrounding_points,
1050 const Table<2, double> &weights,
1051 ArrayView<Point<spacedim>> new_points) const
1052{
1053 Assert(weights.size(0) > 0, ExcEmptyObject());
1054 AssertDimension(surrounding_points.size(), weights.size(1));
1055
1056 const std::size_t n_points = surrounding_points.size();
1057
1058 boost::container::small_vector<Point<chartdim>, 200> chart_points(n_points);
1059 for (std::size_t i = 0; i < n_points; ++i)
1060 chart_points[i] = pull_back(surrounding_points[i]);
1061
1062 boost::container::small_vector<Point<chartdim>, 200> new_points_on_chart(
1063 weights.size(0));
1064 sub_manifold.get_new_points(chart_points, weights, new_points_on_chart);
1065
1066 for (std::size_t row = 0; row < weights.size(0); ++row)
1067 new_points[row] = push_forward(new_points_on_chart[row]);
1068}
1069
1070
1071
1072template <int dim, int spacedim, int chartdim>
1075 const Point<chartdim> &) const
1076{
1077 // function must be implemented in a derived class to be usable,
1078 // as discussed in this function's documentation
1081}
1082
1083
1084
1085template <int dim, int spacedim, int chartdim>
1088 const Point<spacedim> &x1,
1089 const Point<spacedim> &x2) const
1090{
1092 push_forward_gradient(pull_back(x1));
1093
1094 // ensure that the chart is not singular by asserting that its
1095 // derivative has a positive determinant. we need to make this
1096 // comparison relative to the size of the derivative. since the
1097 // determinant is the product of chartdim factors, take the
1098 // chartdim-th root of it in comparing against the size of the
1099 // derivative
1100 Assert(std::pow(std::abs(F_prime.determinant()), 1. / chartdim) >=
1101 1e-12 * F_prime.norm(),
1102 ExcMessage(
1103 "The derivative of a chart function must not be singular."));
1104
1105 const Tensor<1, chartdim> delta =
1106 sub_manifold.get_tangent_vector(pull_back(x1), pull_back(x2));
1107
1108 Tensor<1, spacedim> result;
1109 for (unsigned int i = 0; i < spacedim; ++i)
1110 result[i] += F_prime[i] * delta;
1111
1112 return result;
1113}
1114
1115
1116
1117template <int dim, int spacedim, int chartdim>
1118const Tensor<1, chartdim> &
1120{
1121 return sub_manifold.get_periodicity();
1122}
1123
1124// explicit instantiations
1125#include "manifold.inst"
1126
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
ChartManifold(const Tensor< 1, chartdim > &periodicity=Tensor< 1, chartdim >())
Definition manifold.cc:1004
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:1048
virtual Tensor< 1, spacedim > get_tangent_vector(const Point< spacedim > &x1, const Point< spacedim > &x2) const override
Definition manifold.cc:1087
virtual DerivativeForm< 1, chartdim, spacedim > push_forward_gradient(const Point< chartdim > &chart_point) const
Definition manifold.cc:1074
const Tensor< 1, chartdim > & get_periodicity() const
Definition manifold.cc:1119
virtual Point< spacedim > get_intermediate_point(const Point< spacedim > &p1, const Point< spacedim > &p2, const double w) const override
Definition manifold.cc:1013
virtual Point< spacedim > get_new_point(const ArrayView< const Point< spacedim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition manifold.cc:1027
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:306
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:111
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) 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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 3 > vertices[4]
unsigned int vertex_indices[2]
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:1691
typename IteratorSelector::quad_iterator quad_iterator
Definition tria.h:1667
typename IteratorSelector::line_iterator line_iterator
Definition tria.h:1643
#define DEAL_II_NOT_IMPLEMENTED()
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
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:191
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:220
::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)