Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00:02+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\}}\)
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 
20 
24 #include <deal.II/fe/fe_wedge_p.h>
25 #include <deal.II/fe/mapping_fe.h>
26 #include <deal.II/fe/mapping_q.h>
27 #include <deal.II/fe/mapping_q1.h>
28 
31 #include <deal.II/grid/tria.h>
32 
33 #include <algorithm>
34 #include <iostream>
35 #include <memory>
36 
38 
39 namespace
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
73  VTK_INVALID = numbers::invalid_unsigned_int
74  };
75 
76  } // namespace VTKCellType
77 
78 } // namespace
79 
81 
84 
87 
88 std::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 
120 template <int dim, int spacedim>
121 std::unique_ptr<Mapping<dim, spacedim>>
122 ReferenceCell::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>>(
136  FE_WedgeP<dim, spacedim>(degree));
137  else
138  {
139  Assert(false, ExcNotImplemented());
140  }
141 
142  return std::make_unique<MappingQ<dim, spacedim>>(degree);
143 }
144 
145 
146 
147 template <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 
185 template <int dim>
187 ReferenceCell::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 
207 template <int dim>
208 const 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 
221  return Quadrature<dim>(vertices);
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 
253 unsigned int
254 ReferenceCell::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 
296 unsigned int
297 ReferenceCell::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 
343 unsigned int
344 ReferenceCell::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 
379 unsigned 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 
411 unsigned 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 
443 unsigned 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 
475 template <>
476 unsigned int
477 ReferenceCell::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 
488 template <>
489 unsigned int
490 ReferenceCell::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 
505 template <>
506 unsigned int
507 ReferenceCell::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 
566 template <>
567 unsigned int
568 ReferenceCell::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 
672 unsigned int
673 ReferenceCell::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 
727 unsigned 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 
811 namespace
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 
949 template <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(),
1056  ExcInternalError());
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:
1074  case ReferenceCells::Line:
1077  for (unsigned int d = 0; d < dim; ++d)
1078  result[d] = std::clamp(result[d], 0.0, 1.0);
1079  // simplices can use the standard definition of a simplex:
1080  break;
1082  result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1083  result[y_index] =
1084  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1085  break;
1087  result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1088  result[y_index] =
1089  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1090  result[z_index] =
1091  std::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:
1096  case ReferenceCells::Wedge:
1097  result[x_index] = std::clamp(result[x_index], 0.0, 1.0);
1098  result[y_index] =
1099  std::clamp(result[y_index], 0.0, 1.0 - result[x_index]);
1100  result[z_index] = std::clamp(result[z_index], 0.0, 1.0);
1101  break;
1103  {
1104  result[x_index] = std::clamp(result[x_index], -1.0, 1.0);
1105  result[y_index] = std::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] = std::clamp(result[z_index], 0.0, 1.0 - x_abs);
1113  else
1114  result[z_index] = std::clamp(result[z_index], 0.0, 1.0 - y_abs);
1115  }
1116  break;
1117  default:
1118  Assert(false, ExcNotImplemented());
1119  }
1120  }
1121  // We should be within 4 * eps of the cell by this point. The roundoff error
1122  // comes from, e.g., computing (1 - x) + x when moving points onto the top of
1123  // a Pyramid.
1125  ExcInternalError());
1126 
1127  return result;
1128 }
1129 
1130 
1131 
1132 std::ostream &
1133 operator<<(std::ostream &out, const ReferenceCell &reference_cell)
1134 {
1135  AssertThrow(out.fail() == false, ExcIO());
1136 
1137  // Output as an integer to avoid outputting it as a character with
1138  // potentially non-printing value:
1139  out << static_cast<unsigned int>(reference_cell.kind);
1140  return out;
1141 }
1142 
1143 
1144 
1145 std::istream &
1147 {
1148  AssertThrow(in.fail() == false, ExcIO());
1149 
1150  // Read the information as an integer and convert it to the correct type
1151  unsigned int value;
1152  in >> value;
1153  reference_cell.kind = static_cast<decltype(reference_cell.kind)>(value);
1154 
1155  // Ensure that the object we read is valid
1156  Assert(
1166  ExcMessage(
1167  "The reference cell kind just read does not correspond to one of the "
1168  "valid choices. There must be an error."));
1169 
1170  return in;
1171 }
1172 
1173 // explicitly instantiate dimension 0 quadrature in addition to the standard
1174 // dimensions
1175 template Quadrature<0>
1176 ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const;
1177 template const Quadrature<0> &
1179 
1180 #include "reference_cell.inst"
1181 
OutputOperator< VectorType > & operator<<(OutputOperator< VectorType > &out, unsigned int step)
Definition: operator.h:165
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: point.h:112
numbers::NumberTraits< Number >::real_type distance_square(const Point< dim, Number > &p) const
std::istream & operator>>(std::istream &in, Point< dim, Number > &p)
Definition: point.h:699
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 DEAL_II_HOST SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &t)
Definition: tensor.h:516
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
const unsigned int v0
Definition: grid_tools.cc:1063
const unsigned int v1
Definition: grid_tools.cc:1063
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void shift(const Tensor< 1, spacedim > &shift_vector, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2132
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition: ndarray.h:108