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