deal.II version GIT relicensing-3307-g81a1c05a67 2025-05-12 20:50: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
fe_simplex_p.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
15#include <deal.II/base/config.h>
16
20#include <deal.II/base/types.h>
21
22#include <deal.II/fe/fe_dgq.h>
24#include <deal.II/fe/fe_q.h>
26#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/mapping.h>
28
31
33
34namespace
35{
36 // TODO: replace this with QProjector once QProjector learns how to project
37 // quadrature points onto separate faces of simplices
38 template <int dim>
40 face_midpoint(const unsigned int face_no)
41 {
42 const auto reference_cell = ReferenceCells::get_simplex<dim>();
43 const auto face_reference_cell =
44 reference_cell.face_reference_cell(face_no);
45
46 Point<dim> midpoint;
47 for (const auto face_vertex_no : face_reference_cell.vertex_indices())
48 {
49 const auto vertex_no = reference_cell.face_to_cell_vertices(
50 face_no, face_vertex_no, numbers::default_geometric_orientation);
51
52 midpoint += reference_cell.template vertex<dim>(vertex_no);
53 }
54
55 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
56 return midpoint;
57 }
58
63 std::vector<unsigned int>
64 get_dpo_vector_fe_p(const unsigned int dim, const unsigned int degree)
65 {
66 switch (dim)
67 {
68 case 1:
69 switch (degree)
70 {
71 case 1:
72 return {1, 0};
73 case 2:
74 return {1, 1};
75 case 3:
76 return {1, 2};
77 default:
79 }
80 case 2:
81 switch (degree)
82 {
83 case 1:
84 return {1, 0, 0};
85 case 2:
86 return {1, 1, 0};
87 case 3:
88 return {1, 2, 1};
89 default:
91 }
92 case 3:
93 switch (degree)
94 {
95 case 1:
96 return {1, 0, 0, 0};
97 case 2:
98 return {1, 1, 0, 0};
99 case 3:
100 return {1, 2, 1, 0};
101 default:
103 }
104 }
105
107 return {};
109
110
111
116 template <int dim>
117 std::vector<Point<dim>>
118 unit_support_points_fe_p(const unsigned int degree)
119 {
120 Assert(dim != 0, ExcInternalError());
121 Assert(degree <= 3, ExcNotImplemented());
122 std::vector<Point<dim>> unit_points;
123 const auto reference_cell = ReferenceCells::get_simplex<dim>();
125 // Piecewise constants are a special case: use a support point at the
126 // centroid and only the centroid
127 if (degree == 0)
128 {
129 unit_points.emplace_back(reference_cell.template barycenter<dim>());
130 return unit_points;
131 }
132
133 // otherwise write everything as linear combinations of vertices
134 const auto dpo = get_dpo_vector_fe_p(dim, degree);
135 Assert(dpo.size() == dim + 1, ExcInternalError());
136 Assert(dpo[0] == 1, ExcNotImplemented());
137 // vertices:
138 for (const unsigned int d : reference_cell.vertex_indices())
139 unit_points.push_back(reference_cell.template vertex<dim>(d));
140 // lines:
141 for (const unsigned int l : reference_cell.line_indices())
142 {
143 const Point<dim> p0 =
144 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
145 const Point<dim> p1 =
146 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
147 for (unsigned int p = 0; p < dpo[1]; ++p)
148 unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
149 (double(p + 1) / (dpo[1] + 1)) * p1);
150 }
151 // quads:
152 if (dim > 1 && dpo[2] > 0)
153 {
154 Assert(dpo[2] == 1, ExcNotImplemented());
155 if (dim == 2)
156 unit_points.push_back(reference_cell.template barycenter<dim>());
157 if (dim == 3)
158 for (const unsigned int f : reference_cell.face_indices())
159 unit_points.push_back(face_midpoint<dim>(f));
160 }
161
162 return unit_points;
163 }
164
165 template <>
166 std::vector<Point<0>>
167 unit_support_points_fe_p(const unsigned int /*degree*/)
168 {
169 return {Point<0>()};
170 }
171
176 template <int dim>
177 std::vector<std::vector<Point<dim - 1>>>
178 unit_face_support_points_fe_p(
179 const unsigned int degree,
180 typename FiniteElementData<dim>::Conformity conformity)
181 {
182 // Discontinuous elements don't have face support points
184 return {};
185
186 // this concept doesn't exist in 1d so just return an empty vector
187 if (dim == 1)
188 return {};
189
190 std::vector<std::vector<Point<dim - 1>>> unit_face_points;
191
192 // all faces have the same support points
193 for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
194 {
195 (void)face_n;
196 unit_face_points.emplace_back(
197 unit_support_points_fe_p<dim - 1>(degree));
198 }
199
200 return unit_face_points;
201 }
202
208 template <int dim>
210 constraints_fe_p(const unsigned int /*degree*/)
211 {
212 // no constraints in 1d
213 // constraints in 3d not implemented yet
214 return FullMatrix<double>();
215 }
216
217 template <>
219 constraints_fe_p<2>(const unsigned int degree)
220 {
221 constexpr int dim = 2;
222
223 // the following implements the 2d case
224 // (the 3d case is not implemented yet)
225 //
226 // consult FE_Q_Base::Implementation::initialize_constraints()
227 // for more information
228
229 std::vector<Point<dim - 1>> constraint_points;
230 // midpoint
231 constraint_points.emplace_back(0.5);
232 // subface 0
233 for (unsigned int i = 1; i < degree; ++i)
234 constraint_points.push_back(
236 Point<dim - 1>(i / double(degree)), 0));
237 // subface 1
238 for (unsigned int i = 1; i < degree; ++i)
239 constraint_points.push_back(
241 Point<dim - 1>(i / double(degree)), 1));
242
243 // Now construct relation between destination (child) and source (mother)
244 // dofs.
245
246 const unsigned int n_dofs_constrained = constraint_points.size();
247 unsigned int n_dofs_per_face = degree + 1;
248 FullMatrix<double> interface_constraints(n_dofs_constrained,
249 n_dofs_per_face);
250
251 const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
252
253 for (unsigned int i = 0; i < n_dofs_constrained; ++i)
254 for (unsigned int j = 0; j < n_dofs_per_face; ++j)
255 {
256 interface_constraints(i, j) =
257 poly.compute_value(j, constraint_points[i]);
258
259 // if the value is small up to round-off, then simply set it to zero
260 // to avoid unwanted fill-in of the constraint matrices (which would
261 // then increase the number of other DoFs a constrained DoF would
262 // couple to)
263 if (std::fabs(interface_constraints(i, j)) < 1e-13)
264 interface_constraints(i, j) = 0;
265 }
266 return interface_constraints;
267 }
268
269
270
275 std::vector<unsigned int>
276 get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
277 {
278 // First treat the case of piecewise constant elements:
279 if (degree == 0)
280 {
281 std::vector<unsigned int> dpo(dim + 1, 0U);
282 dpo[dim] = 1;
283 return dpo;
284 }
285 else
286 {
287 // This element has the same degrees of freedom as the continuous one,
288 // but they are all counted for the interior of the cell because
289 // it is continuous. Rather than hard-code how many DoFs the element
290 // has, we just get the numbers from the continuous case and add them
291 // up
292 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
293
294 switch (dim)
295 {
296 case 1:
297 return {0U,
298 ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
299 continuous_dpo[dim]};
300
301 case 2:
302 return {0U,
303 0U,
305 continuous_dpo[0] +
306 ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
307 continuous_dpo[dim]};
308
309 case 3:
310 return {
311 0U,
312 0U,
313 0U,
314 ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
315 ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
316 ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
317 continuous_dpo[dim]};
318 }
319
321 return {};
322 }
323 }
324} // namespace
325
326
327
328template <int dim, int spacedim>
330 const BarycentricPolynomials<dim> polynomials,
331 const FiniteElementData<dim> &fe_data,
332 const bool prolongation_is_additive,
333 const std::vector<Point<dim>> &unit_support_points,
334 const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
335 const FullMatrix<double> &interface_constraints)
336 : ::FE_Poly<dim, spacedim>(
337 polynomials,
338 fe_data,
339 std::vector<bool>(fe_data.dofs_per_cell, prolongation_is_additive),
340 std::vector<ComponentMask>(fe_data.dofs_per_cell,
341 ComponentMask(std::vector<bool>(1, true))))
342{
345 this->interface_constraints = interface_constraints;
346}
347
348
349
350template <int dim, int spacedim>
351std::pair<Table<2, bool>, std::vector<unsigned int>>
353{
354 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
355 constant_modes.fill(true);
356 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
357 constant_modes, std::vector<unsigned int>(1, 0));
358}
359
360
361
362template <int dim, int spacedim>
363const FullMatrix<double> &
365 const unsigned int child,
366 const RefinementCase<dim> &refinement_case) const
367{
368 if (dim == 3)
369 Assert(RefinementCase<dim>(refinement_case) ==
371 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
372 RefinementCase<dim>(refinement_case) ==
374 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
375 RefinementCase<dim>(refinement_case) ==
377 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
379 else
380 Assert(refinement_case ==
383 AssertDimension(dim, spacedim);
384
385 // initialization upon first request
386 if (this->prolongation[refinement_case - 1][child].n() == 0)
387 {
388 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
389
390 // if matrix got updated while waiting for the lock
391 if (this->prolongation[refinement_case - 1][child].n() ==
392 this->n_dofs_per_cell())
393 return this->prolongation[refinement_case - 1][child];
394
395 // now do the work. need to get a non-const version of data in order to
396 // be able to modify them inside a const function
397 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
398
399 if (dim == 2)
400 {
401 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
403 isotropic_matrices.back().resize(
404 this->reference_cell().n_children(
405 RefinementCase<dim>(refinement_case)),
406 FullMatrix<double>(this->n_dofs_per_cell(),
407 this->n_dofs_per_cell()));
408
409 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
410
411 this_nonconst.prolongation[refinement_case - 1] =
412 std::move(isotropic_matrices.back());
413 }
414 else if (dim == 3)
415 {
416 std::vector<std::vector<FullMatrix<double>>> matrices(
417 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49),
418 std::vector<FullMatrix<double>>(
419 this->reference_cell().n_children(
420 RefinementCase<dim>(refinement_case)),
421 FullMatrix<double>(this->n_dofs_per_cell(),
422 this->n_dofs_per_cell())));
423 FETools::compute_embedding_matrices(*this, matrices, true);
424 for (unsigned int refinement_direction = static_cast<unsigned int>(
426 refinement_direction <=
427 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49);
428 refinement_direction++)
429 this_nonconst.prolongation[refinement_direction - 1] =
430 std::move(matrices[refinement_direction - 1]);
431 }
432 else
434 }
435
436 // finally return the matrix
437 return this->prolongation[refinement_case - 1][child];
438}
439
440
441
442template <int dim, int spacedim>
443unsigned int
445 const unsigned int face_dof_index,
446 const unsigned int face,
447 const types::geometric_orientation combined_orientation) const
448{
449 // Function largely lifted from FE_Q_Base::face_to_cell_index()
450 AssertIndexRange(face_dof_index, this->n_dofs_per_face(face));
451 AssertIndexRange(face, this->reference_cell().n_faces());
452
453 // TODO: once the default orientation is switched to 0 then we can remove this
454 // special case for 1D.
455 if (dim == 1)
456 Assert(combined_orientation == numbers::default_geometric_orientation,
457 ExcMessage("In 1D, all faces must have the default orientation."));
458 else
459 AssertIndexRange(combined_orientation,
460 this->reference_cell().n_face_orientations(face));
461
462 // we need to distinguish between DoFs on vertices, lines and in 3d quads.
463 // do so in a sequence of if-else statements
464 if (face_dof_index < this->get_first_face_line_index(face))
465 // DoF is on a vertex
466 {
467 // get the number of the vertex on the face that corresponds to this DoF,
468 // along with the number of the DoF on this vertex
469 const unsigned int face_vertex =
470 face_dof_index / this->n_dofs_per_vertex();
471 const unsigned int dof_index_on_vertex =
472 face_dof_index % this->n_dofs_per_vertex();
473
474 // then get the number of this vertex on the cell and translate
475 // this to a DoF number on the cell
476 return this->reference_cell().face_to_cell_vertices(
477 face, face_vertex, combined_orientation) *
478 this->n_dofs_per_vertex() +
479 dof_index_on_vertex;
480 }
481 else if (face_dof_index < this->get_first_face_quad_index(face))
482 // DoF is on a line
483 {
484 // do the same kind of translation as before. we need to only consider
485 // DoFs on the lines, i.e., ignoring those on the vertices
486 const unsigned int index =
487 face_dof_index - this->get_first_face_line_index(face);
488
489 const unsigned int face_line = index / this->n_dofs_per_line();
490 const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
491
492 // we now also need to adjust the line index for the case of
493 // face orientation, face flips, etc
494 unsigned int adjusted_dof_index_on_line = 0;
495 switch (dim)
496 {
497 case 1:
499 break;
500
501 case 2:
502 if (combined_orientation == numbers::default_geometric_orientation)
503 adjusted_dof_index_on_line = dof_index_on_line;
504 else
505 adjusted_dof_index_on_line =
506 this->n_dofs_per_line() - dof_index_on_line - 1;
507 break;
508
509 case 3:
510 // in 3d, things are difficult. someone will have to think
511 // about how this code here should look like, by drawing a bunch
512 // of pictures of how all the faces can look like with the various
513 // flips and rotations.
514 //
515 // that said, we can implement a couple other situations easily:
516 // if the face orientation is
517 // numbers::default_combined_face_orientation then things are
518 // simple. Likewise if the face orientation is
519 // internal::combined_face_orientation(false,false,false) then we
520 // just know we need to reverse the DoF order on each line.
521 //
522 // this at least allows for tetrahedral meshes where periodic face
523 // pairs are each in one of the two above configurations. Doing so
524 // may turn out to be always possible or it may not, but at least
525 // this is less restrictive.
526 Assert((this->n_dofs_per_line() <= 1) ||
527 combined_orientation ==
529 combined_orientation ==
530 internal::combined_face_orientation(false, false, false),
532
533 if (combined_orientation == numbers::default_geometric_orientation)
534 {
535 adjusted_dof_index_on_line = dof_index_on_line;
536 }
537 else
538 {
539 adjusted_dof_index_on_line =
540 this->n_dofs_per_line() - dof_index_on_line - 1;
541 }
542 break;
543
544 default:
546 }
547
548 return (this->get_first_line_index() +
549 this->reference_cell().face_to_cell_lines(face,
550 face_line,
551 combined_orientation) *
552 this->n_dofs_per_line() +
553 adjusted_dof_index_on_line);
554 }
555 else
556 // DoF is on a quad
557 {
558 Assert(dim >= 3, ExcInternalError());
559
560 // ignore vertex and line dofs
561 const unsigned int index =
562 face_dof_index - this->get_first_face_quad_index(face);
563
564 // the same is true here as above for the 3d case -- someone will
565 // just have to draw a bunch of pictures. in the meantime,
566 // we can implement the degree <= 3 case in which it is simple
567 Assert((this->n_dofs_per_quad(face) <= 1) ||
568 combined_orientation == numbers::default_geometric_orientation,
570 return (this->get_first_quad_index(face) + index);
571 }
572}
573
574
575
576template <int dim, int spacedim>
577const FullMatrix<double> &
579 const unsigned int child,
580 const RefinementCase<dim> &refinement_case) const
581{
582 if (dim == 3)
583 Assert(RefinementCase<dim>(refinement_case) ==
585 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
586 RefinementCase<dim>(refinement_case) ==
588 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
589 RefinementCase<dim>(refinement_case) ==
591 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
593 else
596 AssertDimension(dim, spacedim);
597
598 // initialization upon first request
599 if (this->restriction[refinement_case - 1][child].n() == 0)
600 {
601 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
602
603 // if matrix got updated while waiting for the lock
604 if (this->restriction[refinement_case - 1][child].n() ==
605 this->n_dofs_per_cell())
606 return this->restriction[refinement_case - 1][child];
607
608 // get the restriction matrix
609 // Refine a unit cell. As the parent cell is a unit
610 // cell, the reference cell of the children equals the parent, i.e. they
611 // have the support points at the same locations. So we just have to check
612 // if a support point of the parent is one of the interpolation points of
613 // the child. If this is not the case we find the interpolation of the
614 // point.
615
616 const double eps = 1e-12;
617 FullMatrix<double> restriction_mat(this->n_dofs_per_cell(),
618 this->n_dofs_per_cell());
619
620 // first get all support points on the reference cell
621 const std::vector<Point<dim>> unit_support_points =
622 this->get_unit_support_points();
623
624 // now create children on the reference cell
626 GridGenerator::reference_cell(tria, this->reference_cell());
627 tria.begin_active()->set_refine_flag(
629 if (dim == 3)
630 tria.begin_active()->set_refine_choice(refinement_case);
632
633 const auto &child_cell = tria.begin(0)->child(child);
634
635 // iterate over all support points and transform them to the unit cell of
636 // the child
637 for (unsigned int i = 0; i < unit_support_points.size(); i++)
638 {
639 std::vector<Point<dim>> transformed_point(1);
640 const std::vector<Point<spacedim>> unit_support_point = {
641 dim == 2 ? Point<spacedim>(unit_support_points[i][0],
642 unit_support_points[i][1]) :
643 Point<spacedim>(unit_support_points[i][0],
644 unit_support_points[i][1],
645 unit_support_points[i][2])};
646 this->reference_cell()
647 .template get_default_linear_mapping<dim, spacedim>()
648 .transform_points_real_to_unit_cell(
649 child_cell,
650 make_array_view(unit_support_point),
651 make_array_view(transformed_point));
652
653 // if point is inside the unit cell iterate over all shape functions
654 if (this->reference_cell().contains_point(transformed_point[0], eps))
655 for (unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
656 restriction_mat[i][j] =
657 this->shape_value(j, transformed_point[0]);
658 }
659 if constexpr (running_in_debug_mode())
660 {
661 for (unsigned int i = 0; i < this->n_dofs_per_cell(); i++)
662 {
663 double sum = 0.;
664
665 for (unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
666 sum += restriction_mat[i][j];
667
668 Assert(std::fabs(sum - 1) < eps || std::fabs(sum) < eps,
670 "The entries in a row of the local "
671 "restriction matrix do not add to zero or one. "
672 "This typically indicates that the "
673 "polynomial interpolation is "
674 "ill-conditioned such that round-off "
675 "prevents the sum to be one."));
676 }
677 }
678
679 // Remove small entries from the matrix
680 for (unsigned int i = 0; i < restriction_mat.m(); ++i)
681 for (unsigned int j = 0; j < restriction_mat.n(); ++j)
682 {
683 if (std::fabs(restriction_mat(i, j)) < eps)
684 restriction_mat(i, j) = 0.;
685 if (std::fabs(restriction_mat(i, j) - 1) < eps)
686 restriction_mat(i, j) = 1.;
687 }
688
689 const_cast<FullMatrix<double> &>(
690 this->restriction[refinement_case - 1][child]) =
691 std::move(restriction_mat);
692 }
693
694 // finally return the matrix
695 return this->restriction[refinement_case - 1][child];
696}
697
698
699
700template <int dim, int spacedim>
701void
703 const FiniteElement<dim, spacedim> &source_fe,
704 FullMatrix<double> &interpolation_matrix,
705 const unsigned int face_no) const
706{
707 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
708 ExcDimensionMismatch(interpolation_matrix.m(),
709 source_fe.n_dofs_per_face(face_no)));
710
711 // see if source is a P or Q element
712 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
713 nullptr) ||
714 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
715 {
716 const Quadrature<dim - 1> quad_face_support(
717 source_fe.get_unit_face_support_points(face_no));
718
719 const double eps = 2e-13 * this->degree * (dim - 1);
720
721 const std::vector<Point<dim>> face_quadrature_points =
722 QProjector<dim>::project_to_face(this->reference_cell(),
723 quad_face_support,
724 face_no,
726 .get_points();
727
728 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
729 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
730 {
731 double matrix_entry =
732 this->shape_value(this->face_to_cell_index(j, 0),
733 face_quadrature_points[i]);
734
735 // Correct the interpolated value. I.e. if it is close to 1 or
736 // 0, make it exactly 1 or 0. Unfortunately, this is required to
737 // avoid problems with higher order elements.
738 if (std::fabs(matrix_entry - 1.0) < eps)
739 matrix_entry = 1.0;
740 if (std::fabs(matrix_entry) < eps)
741 matrix_entry = 0.0;
742
743 interpolation_matrix(i, j) = matrix_entry;
744 }
745
746 if constexpr (running_in_debug_mode())
747 {
748 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
749 {
750 double sum = 0.;
751
752 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
753 sum += interpolation_matrix(j, i);
754
755 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
756 }
757 }
758 }
759 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
760 {
761 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
762 }
763 else
765 false,
766 (typename FiniteElement<dim,
767 spacedim>::ExcInterpolationNotImplemented()));
768}
769
770
771
772template <int dim, int spacedim>
773void
775 const FiniteElement<dim, spacedim> &source_fe,
776 const unsigned int subface,
777 FullMatrix<double> &interpolation_matrix,
778 const unsigned int face_no) const
779{
780 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
781 ExcDimensionMismatch(interpolation_matrix.m(),
782 source_fe.n_dofs_per_face(face_no)));
783
784 // see if source is a P or Q element
785 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
786 nullptr) ||
787 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
788 {
789 const Quadrature<dim - 1> quad_face_support(
790 source_fe.get_unit_face_support_points(face_no));
791
792 const double eps = 2e-13 * this->degree * (dim - 1);
793
794 const Quadrature<dim> subface_quadrature =
796 this->reference_cell(),
797 quad_face_support,
798 face_no,
799 subface,
802
803 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
804 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
805 {
806 double matrix_entry =
807 this->shape_value(this->face_to_cell_index(j, 0),
808 subface_quadrature.point(i));
809
810 // Correct the interpolated value. I.e. if it is close to 1 or
811 // 0, make it exactly 1 or 0. Unfortunately, this is required to
812 // avoid problems with higher order elements.
813 if (std::fabs(matrix_entry - 1.0) < eps)
814 matrix_entry = 1.0;
815 if (std::fabs(matrix_entry) < eps)
816 matrix_entry = 0.0;
817
818 interpolation_matrix(i, j) = matrix_entry;
819 }
820
821 if constexpr (running_in_debug_mode())
822 {
823 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
824 {
825 double sum = 0.;
826
827 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
828 sum += interpolation_matrix(j, i);
829
830 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
831 }
832 }
833 }
834 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
835 {
836 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
837 }
838 else
840 false,
841 (typename FiniteElement<dim,
842 spacedim>::ExcInterpolationNotImplemented()));
843}
844
845
846
847template <int dim, int spacedim>
848bool
853
854
855
856template <int dim, int spacedim>
857void
860 const std::vector<Vector<double>> &support_point_values,
861 std::vector<double> &nodal_values) const
862{
863 AssertDimension(support_point_values.size(),
864 this->get_unit_support_points().size());
865 AssertDimension(support_point_values.size(), nodal_values.size());
866 AssertDimension(this->dofs_per_cell, nodal_values.size());
867
868 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
869 {
870 AssertDimension(support_point_values[i].size(), 1);
871
872 nodal_values[i] = support_point_values[i](0);
873 }
874}
875
876
877
878template <int dim, int spacedim>
880 : FE_SimplexPoly<dim, spacedim>(
881 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
882 FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
883 ReferenceCells::get_simplex<dim>(),
884 1,
885 degree,
886 FiniteElementData<dim>::H1),
887 false,
888 unit_support_points_fe_p<dim>(degree),
889 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
890 constraints_fe_p<dim>(degree))
891{
892 if (degree > 2)
893 for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
895 this->n_dofs_per_line() - 1 - i - i;
896 // We do not support multiple DoFs per quad yet
898}
899
900
901
902template <int dim, int spacedim>
903std::unique_ptr<FiniteElement<dim, spacedim>>
905{
906 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
907}
908
909
910
911template <int dim, int spacedim>
912std::string
914{
915 std::ostringstream namebuf;
916 namebuf << "FE_SimplexP<" << Utilities::dim_string(dim, spacedim) << ">("
917 << this->degree << ")";
918
919 return namebuf.str();
920}
921
922
923
924template <int dim, int spacedim>
927 const FiniteElement<dim, spacedim> &fe_other,
928 const unsigned int codim) const
929{
930 Assert(codim <= dim, ExcImpossibleInDim(dim));
931
932 // vertex/line/face domination
933 // (if fe_other is derived from FE_SimplexDGP)
934 // ------------------------------------
935 if (codim > 0)
936 if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
937 nullptr)
938 // there are no requirements between continuous and discontinuous
939 // elements
941
942 // vertex/line/face domination
943 // (if fe_other is not derived from FE_SimplexDGP)
944 // & cell domination
945 // ----------------------------------------
946 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
947 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
948 {
949 if (this->degree < fe_p_other->degree)
951 else if (this->degree == fe_p_other->degree)
953 else
955 }
956 else if (const FE_Q<dim, spacedim> *fe_q_other =
957 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
958 {
959 if (this->degree < fe_q_other->degree)
961 else if (this->degree == fe_q_other->degree)
963 else
965 }
966 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
967 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
968 {
969 if (fe_nothing->is_dominating())
971 else
972 // the FE_Nothing has no degrees of freedom and it is typically used
973 // in a context where we don't require any continuity along the
974 // interface
976 }
977
980}
981
982
983
984template <int dim, int spacedim>
985std::vector<std::pair<unsigned int, unsigned int>>
987 const FiniteElement<dim, spacedim> &fe_other) const
988{
989 AssertDimension(dim, 2);
990
991 if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
992 {
993 // there should be exactly one single DoF of each FE at a vertex, and
994 // they should have identical value
995 return {{0U, 0U}};
996 }
997 else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
998 {
999 // there should be exactly one single DoF of each FE at a vertex, and
1000 // they should have identical value
1001 return {{0U, 0U}};
1002 }
1003 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
1004 {
1005 // the FE_Nothing has no degrees of freedom, so there are no
1006 // equivalencies to be recorded
1007 return {};
1008 }
1009 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
1010 {
1011 // if the other element has no DoFs on faces at all,
1012 // then it would be impossible to enforce any kind of
1013 // continuity even if we knew exactly what kind of element
1014 // we have -- simply because the other element declares
1015 // that it is discontinuous because it has no DoFs on
1016 // its faces. in that case, just state that we have no
1017 // constraints to declare
1018 return {};
1019 }
1020 else
1021 {
1023 return {};
1024 }
1025}
1026
1027
1028
1029template <int dim, int spacedim>
1030std::vector<std::pair<unsigned int, unsigned int>>
1032 const FiniteElement<dim, spacedim> &fe_other) const
1033{
1034 AssertDimension(dim, 2);
1035
1036 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
1037 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
1038 {
1039 // dofs are located along lines, so two dofs are identical if they are
1040 // located at identical positions.
1041 // Therefore, read the points in unit_support_points for the
1042 // first coordinate direction. For FE_SimplexP, they are currently
1043 // hard-coded and we iterate over points on the first line which begin
1044 // after the 3 vertex points in the complete list of unit support points
1045
1046 std::vector<std::pair<unsigned int, unsigned int>> identities;
1047
1048 for (unsigned int i = 0; i < this->degree - 1; ++i)
1049 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
1050 if (std::fabs(this->unit_support_points[i + 3][0] -
1051 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
1052 identities.emplace_back(i, j);
1053 else
1054 {
1055 // If nodes are not located in the same place, we have to
1056 // interpolate. This is then not handled through the
1057 // current function, but via interpolation matrices that
1058 // result in constraints, rather than identities. Since
1059 // that happens in a different function, there is nothing
1060 // for us to do here.
1061 }
1062
1063 return identities;
1064 }
1065 else if (const FE_Q<dim, spacedim> *fe_q_other =
1066 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
1067 {
1068 // dofs are located along lines, so two dofs are identical if they are
1069 // located at identical positions. if we had only equidistant points, we
1070 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
1071 // might have other support points (e.g. Gauss-Lobatto
1072 // points). Therefore, read the points in unit_support_points for the
1073 // first coordinate direction. For FE_Q, we take the lexicographic
1074 // ordering of the line support points in the first direction (i.e.,
1075 // x-direction), which we access between index 1 and p-1 (index 0 and p
1076 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
1077 // iterate over points on the first line which begin after the 3 vertex
1078 // points in the complete list of unit support points
1079
1080 const std::vector<unsigned int> &index_map_inverse_q_other =
1081 fe_q_other->get_poly_space_numbering_inverse();
1082
1083 std::vector<std::pair<unsigned int, unsigned int>> identities;
1084
1085 for (unsigned int i = 0; i < this->degree - 1; ++i)
1086 for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
1087 if (std::fabs(this->unit_support_points[i + 3][0] -
1088 fe_q_other->get_unit_support_points()
1089 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
1090 identities.emplace_back(i, j);
1091 else
1092 {
1093 // If nodes are not located in the same place, we have to
1094 // interpolate. This will then also
1095 // capture the case where the FE_Q has a different polynomial
1096 // degree than the current element. In either case, the resulting
1097 // constraints are computed elsewhere, rather than via the
1098 // identities this function returns: Since
1099 // that happens in a different function, there is nothing
1100 // for us to do here.
1101 }
1102
1103 return identities;
1104 }
1105 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
1106 {
1107 // The FE_Nothing has no degrees of freedom, so there are no
1108 // equivalencies to be recorded. (If the FE_Nothing is dominating,
1109 // then this will also leads to constraints, but we are not concerned
1110 // with this here.)
1111 return {};
1112 }
1113 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
1114 {
1115 // if the other element has no elements on faces at all,
1116 // then it would be impossible to enforce any kind of
1117 // continuity even if we knew exactly what kind of element
1118 // we have -- simply because the other element declares
1119 // that it is discontinuous because it has no DoFs on
1120 // its faces. in that case, just state that we have no
1121 // constraints to declare
1122 return {};
1123 }
1124 else
1125 {
1127 return {};
1128 }
1129}
1130
1131
1132
1133template <int dim, int spacedim>
1135 : FE_SimplexPoly<dim, spacedim>(
1136 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
1137 FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
1138 ReferenceCells::get_simplex<dim>(),
1139 1,
1140 degree,
1141 FiniteElementData<dim>::L2),
1142 true,
1143 unit_support_points_fe_p<dim>(degree),
1144 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::L2),
1145 constraints_fe_p<dim>(degree))
1146{}
1147
1148
1149
1150template <int dim, int spacedim>
1151std::unique_ptr<FiniteElement<dim, spacedim>>
1153{
1154 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
1155}
1156
1157
1158
1159template <int dim, int spacedim>
1160std::string
1162{
1163 std::ostringstream namebuf;
1164 namebuf << "FE_SimplexDGP<" << Utilities::dim_string(dim, spacedim) << ">("
1165 << this->degree << ")";
1166
1167 return namebuf.str();
1168}
1169
1170
1171template <int dim, int spacedim>
1174 const FiniteElement<dim, spacedim> &fe_other,
1175 const unsigned int codim) const
1176{
1177 Assert(codim <= dim, ExcImpossibleInDim(dim));
1178
1179 // vertex/line/face domination
1180 // ---------------------------
1181 if (codim > 0)
1182 // this is a discontinuous element, so by definition there will
1183 // be no constraints wherever this element comes together with
1184 // any other kind of element
1186
1187 // cell domination
1188 // ---------------
1189 if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
1190 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
1191 {
1192 if (this->degree < fe_dgp_other->degree)
1194 else if (this->degree == fe_dgp_other->degree)
1196 else
1198 }
1199 else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
1200 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
1201 {
1202 if (this->degree < fe_dgq_other->degree)
1204 else if (this->degree == fe_dgq_other->degree)
1206 else
1208 }
1209 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
1210 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
1211 {
1212 if (fe_nothing->is_dominating())
1214 else
1215 // the FE_Nothing has no degrees of freedom and it is typically used
1216 // in a context where we don't require any continuity along the
1217 // interface
1219 }
1220
1223}
1224
1225
1226
1227template <int dim, int spacedim>
1228std::vector<std::pair<unsigned int, unsigned int>>
1230 const FiniteElement<dim, spacedim> &fe_other) const
1231{
1232 (void)fe_other;
1233
1234 return {};
1235}
1236
1237
1238
1239template <int dim, int spacedim>
1240std::vector<std::pair<unsigned int, unsigned int>>
1242 const FiniteElement<dim, spacedim> &fe_other) const
1243{
1244 (void)fe_other;
1245
1246 return {};
1247}
1248
1249
1250
1251template <int dim, int spacedim>
1252const FullMatrix<double> &
1254 const unsigned int child,
1255 const RefinementCase<dim> &refinement_case) const
1256{
1257 if (dim == 3)
1258 Assert(RefinementCase<dim>(refinement_case) ==
1260 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
1261 RefinementCase<dim>(refinement_case) ==
1263 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
1264 RefinementCase<dim>(refinement_case) ==
1266 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
1268 else
1271 AssertDimension(dim, spacedim);
1272
1273 // initialization upon first request
1274 if (this->restriction[refinement_case - 1][child].n() == 0)
1275 {
1276 std::lock_guard<std::mutex> lock(this->restriction_matrix_mutex);
1277
1278 // if matrix got updated while waiting for the lock
1279 if (this->restriction[refinement_case - 1][child].n() ==
1280 this->n_dofs_per_cell())
1281 return this->restriction[refinement_case - 1][child];
1282
1283 // now do the work. need to get a non-const version of data in order to
1284 // be able to modify them inside a const function
1285 auto &this_nonconst = const_cast<FE_SimplexDGP<dim, spacedim> &>(*this);
1286
1287 if (dim == 2)
1288 {
1289 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1291 isotropic_matrices.back().resize(
1292 this->reference_cell().n_children(
1293 RefinementCase<dim>(refinement_case)),
1294 FullMatrix<double>(this->n_dofs_per_cell(),
1295 this->n_dofs_per_cell()));
1296
1297 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
1298
1299 this_nonconst.restriction[refinement_case - 1] =
1300 std::move(isotropic_matrices.back());
1301 }
1302 else if (dim == 3)
1303 {
1304 std::vector<std::vector<FullMatrix<double>>> matrices(
1305 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49),
1306 std::vector<FullMatrix<double>>(
1307 this->reference_cell().n_children(
1308 RefinementCase<dim>(refinement_case)),
1309 FullMatrix<double>(this->n_dofs_per_cell(),
1310 this->n_dofs_per_cell())));
1311 FETools::compute_projection_matrices(*this, matrices, true);
1312 for (unsigned int refinement_direction = static_cast<unsigned int>(
1314 refinement_direction <=
1315 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49);
1316 refinement_direction++)
1317 this_nonconst.restriction[refinement_direction - 1] =
1318 std::move(matrices[refinement_direction - 1]);
1319 }
1320 else
1322 }
1323
1324 // finally return the matrix
1325 return this->restriction[refinement_case - 1][child];
1326}
1327
1328// explicit instantiations
1329#include "fe/fe_simplex_p.inst"
1330
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition fe_q.h:554
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::string get_name() const override
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_SimplexDGP(const unsigned int degree)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_SimplexP(const unsigned int degree)
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
std::string get_name() const override
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_SimplexPoly(const BarycentricPolynomials< dim > polynomials, const FiniteElementData< dim > &fe_data, const bool prolongation_is_additive, const std::vector< Point< dim > > &unit_support_points, const std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points, const FullMatrix< double > &interface_constraints)
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &x_source_fe, const unsigned int subface, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
bool hp_constraints_are_implemented() const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const types::geometric_orientation combined_orientation=numbers::default_geometric_orientation) const override
void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source_fe, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition fe.h:2569
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
std::vector< int > adjust_line_dof_index_for_line_orientation_table
Definition fe.h:2611
std::vector< Point< dim > > unit_support_points
Definition fe.h:2562
FullMatrix< double > interface_constraints
Definition fe.h:2550
size_type n() const
size_type m() const
Definition point.h:113
static void project_to_subface(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim > > &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
unsigned int n_vertices() const
unsigned int n_faces() const
unsigned int n_lines() const
cell_iterator begin(const unsigned int level=0) const
virtual void execute_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int vertex_indices[2]
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::size_t size
Definition mpi.cc:745
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
constexpr char U
constexpr ReferenceCell Triangle
constexpr ReferenceCell Tetrahedron
constexpr const ReferenceCell & get_simplex()
constexpr ReferenceCell Line
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:551
types::geometric_orientation combined_face_orientation(const bool face_orientation, const bool face_rotation, const bool face_flip)
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
STL namespace.
unsigned char geometric_orientation
Definition types.h:40