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