Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
reference_cell.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
21
29
32#include <deal.II/grid/tria.h>
33
34#include <iostream>
35#include <memory>
36
38
39namespace
40{
41 namespace VTKCellType
42 {
43 // Define VTK constants for linear, quadratic and
44 // high-order Lagrange geometrices
45 enum : unsigned int
46 {
47 VTK_VERTEX = 1,
48 // Linear cells
49 VTK_LINE = 3,
50 VTK_TRIANGLE = 5,
51 VTK_QUAD = 9,
52 VTK_TETRA = 10,
53 VTK_HEXAHEDRON = 12,
54 VTK_WEDGE = 13,
55 VTK_PYRAMID = 14,
56 // Quadratic cells
57 VTK_QUADRATIC_EDGE = 21,
58 VTK_QUADRATIC_TRIANGLE = 22,
59 VTK_QUADRATIC_QUAD = 23,
60 VTK_QUADRATIC_TETRA = 24,
61 VTK_QUADRATIC_HEXAHEDRON = 25,
62 VTK_QUADRATIC_WEDGE = 26,
63 VTK_QUADRATIC_PYRAMID = 27,
64 // Lagrange cells
65 VTK_LAGRANGE_CURVE = 68,
66 VTK_LAGRANGE_TRIANGLE = 69,
67 VTK_LAGRANGE_QUADRILATERAL = 70,
68 VTK_LAGRANGE_TETRAHEDRON = 71,
69 VTK_LAGRANGE_HEXAHEDRON = 72,
70 VTK_LAGRANGE_WEDGE = 73,
71 VTK_LAGRANGE_PYRAMID = 74,
72 // Invalid code
74 };
75
76 } // namespace VTKCellType
77
78} // namespace
79
81
84
87
88std::string
90{
91 switch (this->kind)
92 {
94 return "Vertex";
96 return "Line";
98 return "Tri";
100 return "Quad";
102 return "Tet";
104 return "Pyramid";
106 return "Wedge";
108 return "Hex";
110 return "Invalid";
111 default:
112 Assert(false, ExcNotImplemented());
113 }
114
115 return "Invalid";
116}
117
118
119
120template <int dim, int spacedim>
121std::unique_ptr<Mapping<dim, spacedim>>
122ReferenceCell::get_default_mapping(const unsigned int degree) const
123{
125
126 if (is_hyper_cube())
127 return std::make_unique<MappingQ<dim, spacedim>>(degree);
128 else if (is_simplex())
129 return std::make_unique<MappingFE<dim, spacedim>>(
131 else if (*this == ReferenceCells::Pyramid)
132 return std::make_unique<MappingFE<dim, spacedim>>(
134 else if (*this == ReferenceCells::Wedge)
135 return std::make_unique<MappingFE<dim, spacedim>>(
137 else
138 {
139 Assert(false, ExcNotImplemented());
140 }
141
142 return std::make_unique<MappingQ<dim, spacedim>>(degree);
143}
144
145
146
147template <int dim, int spacedim>
150{
152
153 if (is_hyper_cube())
154 {
156 }
157 else if (is_simplex())
158 {
159 static const MappingFE<dim, spacedim> mapping(
161 return mapping;
162 }
163 else if (*this == ReferenceCells::Pyramid)
164 {
165 static const MappingFE<dim, spacedim> mapping(
167 return mapping;
168 }
169 else if (*this == ReferenceCells::Wedge)
170 {
171 static const MappingFE<dim, spacedim> mapping(
173 return mapping;
174 }
175 else
176 {
177 Assert(false, ExcNotImplemented());
178 }
179
180 return StaticMappingQ1<dim, spacedim>::mapping; // never reached
181}
182
183
184
185template <int dim>
187ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1d) const
188{
190
191 if (is_hyper_cube())
192 return QGauss<dim>(n_points_1d);
193 else if (is_simplex())
194 return QGaussSimplex<dim>(n_points_1d);
195 else if (*this == ReferenceCells::Pyramid)
196 return QGaussPyramid<dim>(n_points_1d);
197 else if (*this == ReferenceCells::Wedge)
198 return QGaussWedge<dim>(n_points_1d);
199 else
200 Assert(false, ExcNotImplemented());
201
202 return Quadrature<dim>(); // never reached
203}
204
205
206
207template <int dim>
208const Quadrature<dim> &
210{
212
213 // A function that is used to fill a quadrature object of the
214 // desired type the first time we encounter a particular
215 // reference cell
216 const auto create_quadrature = [](const ReferenceCell &reference_cell) {
217 std::vector<Point<dim>> vertices(reference_cell.n_vertices());
218 for (const unsigned int v : reference_cell.vertex_indices())
219 vertices[v] = reference_cell.vertex<dim>(v);
220
222 };
223
224 if (is_hyper_cube())
225 {
226 static const Quadrature<dim> quadrature = create_quadrature(*this);
227 return quadrature;
228 }
229 else if (is_simplex())
230 {
231 static const Quadrature<dim> quadrature = create_quadrature(*this);
232 return quadrature;
233 }
234 else if (*this == ReferenceCells::Pyramid)
235 {
236 static const Quadrature<dim> quadrature = create_quadrature(*this);
237 return quadrature;
238 }
239 else if (*this == ReferenceCells::Wedge)
240 {
241 static const Quadrature<dim> quadrature = create_quadrature(*this);
242 return quadrature;
243 }
244 else
245 Assert(false, ExcNotImplemented());
246
247 static const Quadrature<dim> dummy;
248 return dummy; // never reached
249}
250
251
252
253unsigned int
254ReferenceCell::exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
255{
256 AssertIndexRange(vertex_n, n_vertices());
257
258 switch (this->kind)
259 {
262 return vertex_n;
264 {
265 constexpr std::array<unsigned int, 4> exodus_to_deal{{0, 1, 3, 2}};
266 return exodus_to_deal[vertex_n];
267 }
269 return vertex_n;
271 {
272 constexpr std::array<unsigned int, 8> exodus_to_deal{
273 {0, 1, 3, 2, 4, 5, 7, 6}};
274 return exodus_to_deal[vertex_n];
275 }
277 {
278 constexpr std::array<unsigned int, 6> exodus_to_deal{
279 {2, 1, 0, 5, 4, 3}};
280 return exodus_to_deal[vertex_n];
281 }
283 {
284 constexpr std::array<unsigned int, 5> exodus_to_deal{{0, 1, 3, 2, 4}};
285 return exodus_to_deal[vertex_n];
286 }
287 default:
288 Assert(false, ExcNotImplemented());
289 }
290
292}
293
294
295
296unsigned int
297ReferenceCell::exodusii_face_to_deal_face(const unsigned int face_n) const
298{
299 AssertIndexRange(face_n, n_faces());
300
301 switch (this->kind)
302 {
304 return 0;
307 return face_n;
309 {
310 constexpr std::array<unsigned int, 4> exodus_to_deal{{2, 1, 3, 0}};
311 return exodus_to_deal[face_n];
312 }
314 {
315 constexpr std::array<unsigned int, 4> exodus_to_deal{{1, 3, 2, 0}};
316 return exodus_to_deal[face_n];
317 }
319 {
320 constexpr std::array<unsigned int, 6> exodus_to_deal{
321 {2, 1, 3, 0, 4, 5}};
322 return exodus_to_deal[face_n];
323 }
325 {
326 constexpr std::array<unsigned int, 6> exodus_to_deal{{3, 4, 2, 0, 1}};
327 return exodus_to_deal[face_n];
328 }
330 {
331 constexpr std::array<unsigned int, 5> exodus_to_deal{{3, 2, 4, 1, 0}};
332 return exodus_to_deal[face_n];
333 }
334 default:
335 Assert(false, ExcNotImplemented());
336 }
337
339}
340
341
342
343unsigned int
344ReferenceCell::unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
345{
346 AssertIndexRange(vertex_n, n_vertices());
347 // Information on this file format isn't easy to find - the documents here
348 //
349 // https://www.ceas3.uc.edu/sdrluff/
350 //
351 // Don't actually explain anything about the sections we care about (2412) in
352 // any detail. For node numbering I worked backwards from what is actually in
353 // our test files (since that's supposed to work), which all use some
354 // non-standard clockwise numbering scheme which starts at the bottom right
355 // vertex.
356 if (*this == ReferenceCells::Line)
357 {
358 return vertex_n;
359 }
360 else if (*this == ReferenceCells::Quadrilateral)
361 {
362 constexpr std::array<unsigned int, 4> unv_to_deal{{1, 0, 2, 3}};
363 return unv_to_deal[vertex_n];
364 }
365 else if (*this == ReferenceCells::Hexahedron)
366 {
367 constexpr std::array<unsigned int, 8> unv_to_deal{
368 {6, 7, 5, 4, 2, 3, 1, 0}};
369 return unv_to_deal[vertex_n];
370 }
371
372 Assert(false, ExcNotImplemented());
373
375}
376
377
378
379unsigned int
381{
382 switch (this->kind)
383 {
385 return VTKCellType::VTK_VERTEX;
387 return VTKCellType::VTK_LINE;
389 return VTKCellType::VTK_TRIANGLE;
391 return VTKCellType::VTK_QUAD;
393 return VTKCellType::VTK_TETRA;
395 return VTKCellType::VTK_PYRAMID;
397 return VTKCellType::VTK_WEDGE;
399 return VTKCellType::VTK_HEXAHEDRON;
401 return VTKCellType::VTK_INVALID;
402 default:
403 Assert(false, ExcNotImplemented());
404 }
405
406 return VTKCellType::VTK_INVALID;
407}
408
409
410
411unsigned int
413{
414 switch (this->kind)
415 {
417 return VTKCellType::VTK_VERTEX;
419 return VTKCellType::VTK_QUADRATIC_EDGE;
421 return VTKCellType::VTK_QUADRATIC_TRIANGLE;
423 return VTKCellType::VTK_QUADRATIC_QUAD;
425 return VTKCellType::VTK_QUADRATIC_TETRA;
427 return VTKCellType::VTK_QUADRATIC_PYRAMID;
429 return VTKCellType::VTK_QUADRATIC_WEDGE;
431 return VTKCellType::VTK_QUADRATIC_HEXAHEDRON;
433 return VTKCellType::VTK_INVALID;
434 default:
435 Assert(false, ExcNotImplemented());
436 }
437
438 return VTKCellType::VTK_INVALID;
439}
440
441
442
443unsigned int
445{
446 switch (this->kind)
447 {
449 return VTKCellType::VTK_VERTEX;
451 return VTKCellType::VTK_LAGRANGE_CURVE;
453 return VTKCellType::VTK_LAGRANGE_TRIANGLE;
455 return VTKCellType::VTK_LAGRANGE_QUADRILATERAL;
457 return VTKCellType::VTK_LAGRANGE_TETRAHEDRON;
459 return VTKCellType::VTK_LAGRANGE_PYRAMID;
461 return VTKCellType::VTK_LAGRANGE_WEDGE;
463 return VTKCellType::VTK_LAGRANGE_HEXAHEDRON;
465 return VTKCellType::VTK_INVALID;
466 default:
467 Assert(false, ExcNotImplemented());
468 }
469
470 return VTKCellType::VTK_INVALID;
471}
472
473
474
475template <>
476unsigned int
477ReferenceCell::vtk_lexicographic_to_node_index<0>(
478 const std::array<unsigned, 0> &,
479 const std::array<unsigned, 0> &,
480 const bool) const
481{
482 Assert(false, ExcNotImplemented());
483 return 0;
484}
485
486
487
488template <>
489unsigned int
490ReferenceCell::vtk_lexicographic_to_node_index<1>(
491 const std::array<unsigned, 1> &,
492 const std::array<unsigned, 1> &,
493 const bool) const
494{
495 Assert(false, ExcNotImplemented());
496 return 0;
497}
498
499
500
505template <>
506unsigned int
507ReferenceCell::vtk_lexicographic_to_node_index<2>(
508 const std::array<unsigned, 2> &node_indices,
509 const std::array<unsigned, 2> &nodes_per_direction,
510 const bool) const
511{
513
514 const unsigned int i = node_indices[0];
515 const unsigned int j = node_indices[1];
516
517 const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
518 const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
519 // How many boundaries do we lie on at once?
520 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0);
521
522 if (nbdy == 2) // Vertex DOF
523 { // ijk is a corner node. Return the proper index (somewhere in [0,3]):
524 return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0));
525 }
526
527 int offset = 4;
528 if (nbdy == 1) // Edge DOF
529 {
530 if (!ibdy)
531 { // On i axis
532 return (i - 1) +
533 (j != 0u ?
534 nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 :
535 0) +
536 offset;
537 }
538
539 if (!jbdy)
540 { // On j axis
541 return (j - 1) +
542 (i != 0u ? nodes_per_direction[0] - 1 :
543 2 * (nodes_per_direction[0] - 1) +
544 nodes_per_direction[1] - 1) +
545 offset;
546 }
547 }
548
549 offset += 2 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1);
550 // nbdy == 0: Face DOF
551 return offset + (i - 1) + (nodes_per_direction[0] - 1) * ((j - 1));
552}
553
554
555
566template <>
567unsigned int
568ReferenceCell::vtk_lexicographic_to_node_index<3>(
569 const std::array<unsigned, 3> &node_indices,
570 const std::array<unsigned, 3> &nodes_per_direction,
571 const bool legacy_format) const
572{
574
575 const unsigned int i = node_indices[0];
576 const unsigned int j = node_indices[1];
577 const unsigned int k = node_indices[2];
578
579 const bool ibdy = (i == 0 || i == nodes_per_direction[0]);
580 const bool jbdy = (j == 0 || j == nodes_per_direction[1]);
581 const bool kbdy = (k == 0 || k == nodes_per_direction[2]);
582 // How many boundaries do we lie on at once?
583 const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0);
584
585 if (nbdy == 3) // Vertex DOF
586 { // ijk is a corner node. Return the proper index (somewhere in [0,7]):
587 return (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
588 (k != 0u ? 4 : 0);
589 }
590
591 int offset = 8;
592 if (nbdy == 2) // Edge DOF
593 {
594 if (!ibdy)
595 { // On i axis
596 return (i - 1) +
597 (j != 0u ?
598 nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 :
599 0) +
600 (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
601 nodes_per_direction[1] - 1) :
602 0) +
603 offset;
604 }
605 if (!jbdy)
606 { // On j axis
607 return (j - 1) +
608 (i != 0u ? nodes_per_direction[0] - 1 :
609 2 * (nodes_per_direction[0] - 1) +
610 nodes_per_direction[1] - 1) +
611 (k != 0u ? 2 * (nodes_per_direction[0] - 1 +
612 nodes_per_direction[1] - 1) :
613 0) +
614 offset;
615 }
616 // !kbdy, On k axis
617 offset +=
618 4 * (nodes_per_direction[0] - 1) + 4 * (nodes_per_direction[1] - 1);
619 if (legacy_format)
620 return (k - 1) +
621 (nodes_per_direction[2] - 1) *
622 (i != 0u ? (j != 0u ? 3 : 1) : (j != 0u ? 2 : 0)) +
623 offset;
624 else
625 return (k - 1) +
626 (nodes_per_direction[2] - 1) *
627 (i != 0u ? (j != 0u ? 2 : 1) : (j != 0u ? 3 : 0)) +
628 offset;
629 }
630
631 offset += 4 * (nodes_per_direction[0] - 1 + nodes_per_direction[1] - 1 +
632 nodes_per_direction[2] - 1);
633 if (nbdy == 1) // Face DOF
634 {
635 if (ibdy) // On i-normal face
636 {
637 return (j - 1) + ((nodes_per_direction[1] - 1) * (k - 1)) +
638 (i != 0u ? (nodes_per_direction[1] - 1) *
639 (nodes_per_direction[2] - 1) :
640 0) +
641 offset;
642 }
643 offset += 2 * (nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1);
644 if (jbdy) // On j-normal face
645 {
646 return (i - 1) + ((nodes_per_direction[0] - 1) * (k - 1)) +
647 (j != 0u ? (nodes_per_direction[2] - 1) *
648 (nodes_per_direction[0] - 1) :
649 0) +
650 offset;
651 }
652 offset += 2 * (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1);
653 // kbdy, On k-normal face
654 return (i - 1) + ((nodes_per_direction[0] - 1) * (j - 1)) +
655 (k != 0u ?
656 (nodes_per_direction[0] - 1) * (nodes_per_direction[1] - 1) :
657 0) +
658 offset;
659 }
660
661 // nbdy == 0: Body DOF
662 offset += 2 * ((nodes_per_direction[1] - 1) * (nodes_per_direction[2] - 1) +
663 (nodes_per_direction[2] - 1) * (nodes_per_direction[0] - 1) +
664 (nodes_per_direction[0] - 1) * (nodes_per_direction[1] - 1));
665 return offset + (i - 1) +
666 (nodes_per_direction[0] - 1) *
667 ((j - 1) + (nodes_per_direction[1] - 1) * ((k - 1)));
668}
669
670
671
672unsigned int
673ReferenceCell::vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
674{
675 AssertIndexRange(vertex_index, n_vertices());
676
677 // For some of the following, deal.II uses the same ordering as VTK
678 // and in that case, we only need to return 'vertex_index' (i.e.,
679 // use the identity mapping). For some others, we need to translate.
680 //
681 // For the ordering, see the VTK manual (for example at
682 // http://www.princeton.edu/~efeibush/viscourse/vtk.pdf, page 9).
683 switch (this->kind)
684 {
686 return vertex_index;
688 return vertex_index;
690 return vertex_index;
692 {
693 static constexpr std::array<unsigned int, 4> index_translation_table =
694 {{0, 1, 3, 2}};
695 return index_translation_table[vertex_index];
696 }
698 return vertex_index;
700 {
701 static constexpr std::array<unsigned int, 5> index_translation_table =
702 {{0, 1, 3, 2, 4}};
703 return index_translation_table[vertex_index];
704 }
706 return vertex_index;
708 {
709 static constexpr std::array<unsigned int, 8> index_translation_table =
710 {{0, 1, 3, 2, 4, 5, 7, 6}};
711 return index_translation_table[vertex_index];
712 }
714 {
715 Assert(false, ExcNotImplemented());
717 }
718 default:
719 Assert(false, ExcNotImplemented());
720 }
721
723}
724
725
726
727unsigned int
729{
730 /*
731 From the GMSH documentation:
732
733 elm-type
734 defines the geometrical type of the n-th element:
735
736 1
737 Line (2 nodes).
738
739 2
740 Triangle (3 nodes).
741
742 3
743 Quadrangle (4 nodes).
744
745 4
746 Tetrahedron (4 nodes).
747
748 5
749 Hexahedron (8 nodes).
750
751 6
752 Prism (6 nodes).
753
754 7
755 Pyramid (5 nodes).
756
757 8
758 Second order line (3 nodes: 2 associated with the vertices and 1 with the
759 edge).
760
761 9
762 Second order triangle (6 nodes: 3 associated with the vertices and 3 with
763 the edges).
764
765 10 Second order quadrangle (9 nodes: 4 associated with the
766 vertices, 4 with the edges and 1 with the face).
767
768 11 Second order tetrahedron (10 nodes: 4 associated with the vertices and 6
769 with the edges).
770
771 12 Second order hexahedron (27 nodes: 8 associated with the vertices, 12
772 with the edges, 6 with the faces and 1 with the volume).
773
774 13 Second order prism (18 nodes: 6 associated with the vertices, 9 with the
775 edges and 3 with the quadrangular faces).
776
777 14 Second order pyramid (14 nodes: 5 associated with the vertices, 8 with
778 the edges and 1 with the quadrangular face).
779
780 15 Point (1 node).
781 */
782
783 switch (this->kind)
784 {
786 return 15;
788 return 1;
790 return 2;
792 return 3;
794 return 4;
796 return 7;
798 return 6;
800 return 5;
802 default:
803 Assert(false, ExcNotImplemented());
804 }
805
807}
808
809
810
811namespace
812{
813 // Compute the nearest point to @p on the line segment and the square of its
814 // distance to @p.
815 template <int dim>
816 std::pair<Point<dim>, double>
817 project_to_line(const Point<dim> &x0,
818 const Point<dim> &x1,
819 const Point<dim> &p)
820 {
821 Assert(x0 != x1, ExcInternalError());
822 // t is the convex combination coefficient (x = (1 - t) * x0 + t * x1)
823 // defining the position of the closest point on the line (not line segment)
824 // to p passing through x0 and x1. This formula is equivalent to the
825 // standard 'project a vector onto another vector', where each vector is
826 // shifted to start at x0.
827 const double t = ((x1 - x0) * (p - x0)) / ((x1 - x0).norm_square());
828
829 // Only consider points between x0 and x1
830 if (t <= 0)
831 return std::make_pair(x0, x0.distance_square(p));
832 else if (t <= 1)
833 {
834 const auto p2 = x0 + t * (x1 - x0);
835 return std::make_pair(p2, p2.distance_square(p));
836 }
837 else
838 return std::make_pair(x1, x1.distance_square(p));
839 }
840
841 // template base-case
842 template <int dim>
843 std::pair<Point<dim>, double>
844 project_to_quad(const std::array<Point<dim>, 3> & /*vertices*/,
845 const Point<dim> & /*p*/,
846 const ReferenceCell /*reference_cell*/)
847 {
848 Assert(false, ExcInternalError());
849 return std::make_pair(Point<dim>(),
850 std::numeric_limits<double>::signaling_NaN());
851 }
852
870 template <>
871 std::pair<Point<3>, double>
872 project_to_quad(const std::array<Point<3>, 3> &vertices,
873 const Point<3> & p,
874 const ReferenceCell face_reference_cell)
875 {
876 Assert(face_reference_cell == ReferenceCells::Triangle ||
877 face_reference_cell == ReferenceCells::Quadrilateral,
879
880 // Make the problem slightly easier by shifting everything to avoid a point
881 // at the origin (this way we can invert the matrix of vertices). Use 2.0 so
882 // that the bottom left vertex of a Pyramid is now at x = 1.
883 std::array<Point<3>, 3> shifted_vertices = vertices;
884 const Tensor<1, 3> shift{{2.0, 2.0, 2.0}};
885 for (Point<3> &shifted_vertex : shifted_vertices)
886 shifted_vertex += shift;
887 const Point<3> shifted_p = p + shift;
888
889 // As we are projecting onto a face of a reference cell, the vectors
890 // describing its local coordinate system should be orthogonal. We don't
891 // know which of the three vectors computed from p are mutually orthogonal
892 // for triangles so that case requires an extra check.
893 Tensor<1, 3> e0;
894 Tensor<1, 3> e1;
895 const Point<3> vertex = shifted_vertices[0];
896 // Triangles are difficult because of two cases:
897 // 1. the top face of a Tetrahedron, which does not have a right angle
898 // 2. wedges and pyramids, whose faces do not lie on the reference simplex
899 //
900 // Deal with both by creating a locally orthogonal (but not necessarily
901 // orthonormal) coordinate system and testing if the projected point is in
902 // the triangle by expressing it as a convex combination of the vertices.
903 if (face_reference_cell == ReferenceCells::Triangle)
904 {
905 e0 = shifted_vertices[1] - shifted_vertices[0];
906 e1 = shifted_vertices[2] - shifted_vertices[0];
907 e1 -= (e0 * e1) * e0 / (e0.norm_square());
908 }
909 else
910 {
911 e0 = shifted_vertices[1] - shifted_vertices[0];
912 e1 = shifted_vertices[2] - shifted_vertices[0];
913 }
914 Assert(std::abs(e0 * e1) <= 1e-14, ExcInternalError());
915 // the quadrilaterals on pyramids and wedges don't necessarily have edge
916 // lengths of 1 so we cannot skip the denominator
917 const double c0 = e0 * (shifted_p - vertex) / e0.norm_square();
918 const double c1 = e1 * (shifted_p - vertex) / e1.norm_square();
919 const Point<3> projected_shifted_p = vertex + c0 * e0 + c1 * e1;
920
921 bool in_quad = false;
922 if (face_reference_cell == ReferenceCells::Triangle)
923 {
924 Tensor<2, 3> shifted_vertex_matrix;
925 for (unsigned int i = 0; i < 3; ++i)
926 shifted_vertex_matrix[i] = shifted_vertices[i];
927 const Tensor<1, 3> combination_coordinates =
928 invert(transpose(shifted_vertex_matrix)) * projected_shifted_p;
929 bool is_convex_combination = true;
930 for (unsigned int i = 0; i < 3; ++i)
931 is_convex_combination = is_convex_combination &&
932 (0.0 <= combination_coordinates[i]) &&
933 (combination_coordinates[i] <= 1.0);
934 in_quad = is_convex_combination;
935 }
936 else
937 in_quad = (0.0 <= c0 && c0 <= 1.0 && 0.0 <= c1 && c1 <= 1.0);
938
939 if (in_quad)
940 return std::make_pair(projected_shifted_p - shift,
941 shifted_p.distance_square(projected_shifted_p));
942 else
943 return std::make_pair(Point<3>(), std::numeric_limits<double>::max());
944 }
945} // namespace
946
947
948
949template <int dim>
952{
954
955 // Handle simple cases first:
956 if (dim == 0)
957 return p;
958 if (contains_point(p, 0.0))
959 return p;
960 if (dim == 1)
961 return project_to_line(vertex<dim>(0), vertex<dim>(1), p).first;
962
963 // Find the closest vertex so that we only need to check adjacent faces and
964 // lines.
965 Point<dim> result;
966 unsigned int closest_vertex_no = 0;
967 double closest_vertex_distance_square = vertex<dim>(0).distance_square(p);
968 for (unsigned int i = 1; i < n_vertices(); ++i)
969 {
970 const double new_vertex_distance_square =
971 vertex<dim>(i).distance_square(p);
972 if (new_vertex_distance_square < closest_vertex_distance_square)
973 {
974 closest_vertex_no = i;
975 closest_vertex_distance_square = new_vertex_distance_square;
976 }
977 }
978
979 double min_distance_square = std::numeric_limits<double>::max();
980 if (dim == 2)
981 {
982 for (const unsigned int face_no :
983 faces_for_given_vertex(closest_vertex_no))
984 {
985 const Point<dim> v0 = vertex<dim>(line_to_cell_vertices(face_no, 0));
986 const Point<dim> v1 = vertex<dim>(line_to_cell_vertices(face_no, 1));
987
988 auto pair = project_to_line(v0, v1, p);
989 if (pair.second < min_distance_square)
990 {
991 result = pair.first;
992 min_distance_square = pair.second;
993 }
994 }
995 }
996 else
997 {
998 // Check faces and then lines.
999 //
1000 // For reference cells with sloped faces (i.e., all 3D shapes except
1001 // Hexahedra), we might be able to do a valid normal projection to a face
1002 // with a different slope which is on the 'other side' of the reference
1003 // cell. To catch that case we have to unconditionally check lines after
1004 // checking faces.
1005 //
1006 // For pyramids the closest vertex might not be on the closest face: for
1007 // example, the origin is closest to vertex 4 which is not on the bottom
1008 // plane. Get around that by just checking all faces for pyramids.
1009 const std::array<unsigned int, 5> all_pyramid_faces{{0, 1, 2, 3, 4}};
1010 const auto &faces = *this == ReferenceCells::Pyramid ?
1011 ArrayView<const unsigned int>(all_pyramid_faces) :
1012 faces_for_given_vertex(closest_vertex_no);
1013 for (const unsigned int face_no : faces)
1014 {
1015 auto face_cell = face_reference_cell(face_no);
1016 // We only need the first three points since for quads the last point
1017 // is redundant
1018 std::array<Point<dim>, 3> vertices;
1019 for (unsigned int vertex_no = 0; vertex_no < 3; ++vertex_no)
1020 vertices[vertex_no] = vertex<dim>(face_to_cell_vertices(
1021 face_no, vertex_no, default_combined_face_orientation()));
1022
1023 auto pair = project_to_quad(vertices, p, face_cell);
1024 if (pair.second < min_distance_square)
1025 {
1026 result = pair.first;
1027 min_distance_square = pair.second;
1028 }
1029 }
1030
1031 for (const unsigned int face_no :
1032 faces_for_given_vertex(closest_vertex_no))
1033 {
1034 auto face_cell = face_reference_cell(face_no);
1035 for (const unsigned int face_line_no : face_cell.line_indices())
1036 {
1037 const auto cell_line_no =
1038 face_to_cell_lines(face_no,
1039 face_line_no,
1041 const auto v0 =
1042 vertex<dim>(line_to_cell_vertices(cell_line_no, 0));
1043 const auto v1 =
1044 vertex<dim>(line_to_cell_vertices(cell_line_no, 1));
1045 auto pair = project_to_line(v0, v1, p);
1046 if (pair.second < min_distance_square)
1047 {
1048 result = pair.first;
1049 min_distance_square = pair.second;
1050 }
1051 }
1052 }
1053 }
1054
1055 Assert(min_distance_square < std::numeric_limits<double>::max(),
1057
1058 // If necessary, slightly adjust the computed point so that it is closer to
1059 // being on the surface of the reference cell. Due to roundoff it is difficult
1060 // to place points on sloped surfaces (e.g., for Pyramids) so this check isn't
1061 // perfect but does improve the accuracy of the projected point.
1062 if (!contains_point(result, 0.0))
1063 {
1064 constexpr unsigned int x_index = 0;
1065 constexpr unsigned int y_index = (dim >= 2 ? 1 : 0);
1066 constexpr unsigned int z_index = (dim >= 3 ? 2 : 0);
1067 switch (this->kind)
1068 {
1070 Assert(false, ExcInternalError());
1071 break;
1072 // the bounds for each dimension of a hypercube are mutually
1073 // independent:
1077 for (unsigned int d = 0; d < dim; ++d)
1078 result[d] = std_cxx17::clamp(result[d], 0.0, 1.0);
1079 // simplices can use the standard definition of a simplex:
1080 break;
1082 result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
1083 result[y_index] =
1084 std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1085 break;
1087 result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
1088 result[y_index] =
1089 std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1090 result[z_index] =
1091 std_cxx17::clamp(result[z_index],
1092 0.0,
1093 1.0 - result[x_index] - result[y_index]);
1094 break;
1095 // wedges and pyramids are more ad-hoc:
1097 result[x_index] = std_cxx17::clamp(result[x_index], 0.0, 1.0);
1098 result[y_index] =
1099 std_cxx17::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1100 result[z_index] = std_cxx17::clamp(result[z_index], 0.0, 1.0);
1101 break;
1103 {
1104 result[x_index] = std_cxx17::clamp(result[x_index], -1.0, 1.0);
1105 result[y_index] = std_cxx17::clamp(result[y_index], -1.0, 1.0);
1106 // It suffices to transform everything to the first quadrant to
1107 // adjust z:
1108 const auto x_abs = std::abs(result[x_index]);
1109 const auto y_abs = std::abs(result[y_index]);
1110
1111 if (y_abs <= x_abs)
1112 result[z_index] =
1113 std_cxx17::clamp(result[z_index], 0.0, 1.0 - x_abs);
1114 else
1115 result[z_index] =
1116 std_cxx17::clamp(result[z_index], 0.0, 1.0 - y_abs);
1117 }
1118 break;
1119 default:
1120 Assert(false, ExcNotImplemented());
1121 }
1122 }
1123 // We should be within 4 * eps of the cell by this point. The roundoff error
1124 // comes from, e.g., computing (1 - x) + x when moving points onto the top of
1125 // a Pyramid.
1126 Assert(contains_point(result, 4.0 * std::numeric_limits<double>::epsilon()),
1128
1129 return result;
1130}
1131
1132
1133
1134std::ostream &
1135operator<<(std::ostream &out, const ReferenceCell &reference_cell)
1136{
1137 AssertThrow(out.fail() == false, ExcIO());
1138
1139 // Output as an integer to avoid outputting it as a character with
1140 // potentially non-printing value:
1141 out << static_cast<unsigned int>(reference_cell.kind);
1142 return out;
1143}
1144
1145
1146
1147std::istream &
1148operator>>(std::istream &in, ReferenceCell &reference_cell)
1149{
1150 AssertThrow(in.fail() == false, ExcIO());
1151
1152 // Read the information as an integer and convert it to the correct type
1153 unsigned int value;
1154 in >> value;
1155 reference_cell.kind = static_cast<decltype(reference_cell.kind)>(value);
1156
1157 // Ensure that the object we read is valid
1158 Assert(
1159 (reference_cell == ReferenceCells::Vertex) ||
1160 (reference_cell == ReferenceCells::Line) ||
1161 (reference_cell == ReferenceCells::Triangle) ||
1162 (reference_cell == ReferenceCells::Quadrilateral) ||
1163 (reference_cell == ReferenceCells::Tetrahedron) ||
1164 (reference_cell == ReferenceCells::Hexahedron) ||
1165 (reference_cell == ReferenceCells::Wedge) ||
1166 (reference_cell == ReferenceCells::Pyramid) ||
1167 (reference_cell == ReferenceCells::Invalid),
1168 ExcMessage(
1169 "The reference cell kind just read does not correspond to one of the "
1170 "valid choices. There must be an error."));
1171
1172 return in;
1173}
1174
1175// explicitly instantiate dimension 0 quadrature in addition to the standard
1176// dimensions
1177template Quadrature<0>
1178ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const;
1179template const Quadrature<0> &
1181
1182#include "reference_cell.inst"
1183
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
ArrayView< const unsigned int > faces_for_given_vertex(const unsigned int vertex_index) const
unsigned int n_vertices() const
Quadrature< dim > get_gauss_type_quadrature(const unsigned n_points_1d) const
bool is_hyper_cube() const
unsigned int vtk_linear_type() const
unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const unsigned char face_orientation) const
static constexpr unsigned char default_combined_face_orientation()
const Quadrature< dim > & get_nodal_type_quadrature() const
static constexpr ndarray< unsigned int, 2, 2 > line_vertex_permutations
bool contains_point(const Point< dim > &p, const double tolerance=0) const
unsigned int gmsh_element_type() const
static constexpr ndarray< unsigned int, 8, 4 > quadrilateral_vertex_permutations
unsigned int n_faces() const
unsigned int line_to_cell_vertices(const unsigned int line, const unsigned int vertex) const
std::uint8_t kind
unsigned int exodusii_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int unv_vertex_to_deal_vertex(const unsigned int vertex_n) const
unsigned int vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
std::unique_ptr< Mapping< dim, spacedim > > get_default_mapping(const unsigned int degree) const
unsigned int vtk_quadratic_type() const
unsigned int vtk_lagrange_type() const
unsigned int get_dimension() const
ReferenceCell face_reference_cell(const unsigned int face_no) const
unsigned int exodusii_face_to_deal_face(const unsigned int face_n) const
const Mapping< dim, spacedim > & get_default_linear_mapping() const
std::string to_string() const
bool is_simplex() const
unsigned int face_to_cell_lines(const unsigned int face, const unsigned int line, const unsigned char face_orientation) const
Point< dim > closest_point(const Point< dim > &p) const
static constexpr ndarray< unsigned int, 6, 3 > triangle_vertex_permutations
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > vertices[4]
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
const unsigned int v0
const unsigned int v1
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Quadrilateral
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Invalid
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Hexahedron
constexpr const ReferenceCell Vertex
constexpr const ReferenceCell Line
static const unsigned int invalid_unsigned_int
Definition types.h:213
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:108
std::istream & operator>>(std::istream &in, ReferenceCell &reference_cell)
std::ostream & operator<<(std::ostream &out, const ReferenceCell &reference_cell)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)