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