Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-2838-gd85d4b70e9 2025-03-13 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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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
19
20#include <deal.II/fe/fe_dgq.h>
22#include <deal.II/fe/fe_q.h>
24#include <deal.II/fe/fe_tools.h>
25#include <deal.II/fe/mapping.h>
26
28
30
31namespace
32{
33 // TODO: replace this with QProjector once QProjector learns how to project
34 // quadrature points onto separate faces of simplices
35 template <int dim>
37 face_midpoint(const unsigned int face_no)
38 {
39 const auto reference_cell = ReferenceCells::get_simplex<dim>();
40 const auto face_reference_cell =
41 reference_cell.face_reference_cell(face_no);
43 Point<dim> midpoint;
44 for (const auto face_vertex_no : face_reference_cell.vertex_indices())
45 {
46 const auto vertex_no = reference_cell.face_to_cell_vertices(
47 face_no, face_vertex_no, numbers::default_geometric_orientation);
48
49 midpoint += reference_cell.template vertex<dim>(vertex_no);
50 }
51
52 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
53 return midpoint;
54 }
60 std::vector<unsigned int>
61 get_dpo_vector_fe_p(const unsigned int dim, const unsigned int degree)
62 {
63 switch (dim)
64 {
65 case 1:
66 switch (degree)
67 {
68 case 1:
69 return {1, 0};
70 case 2:
71 return {1, 1};
72 case 3:
73 return {1, 2};
74 default:
76 }
77 case 2:
78 switch (degree)
79 {
80 case 1:
81 return {1, 0, 0};
82 case 2:
83 return {1, 1, 0};
84 case 3:
85 return {1, 2, 1};
86 default:
88 }
89 case 3:
90 switch (degree)
91 {
92 case 1:
93 return {1, 0, 0, 0};
94 case 2:
95 return {1, 1, 0, 0};
96 case 3:
97 return {1, 2, 1, 0};
98 default:
100 }
102
104 return {};
105 }
106
108
113 template <int dim>
114 std::vector<Point<dim>>
115 unit_support_points_fe_p(const unsigned int degree)
116 {
117 Assert(dim != 0, ExcInternalError());
118 Assert(degree <= 3, ExcNotImplemented());
119 std::vector<Point<dim>> unit_points;
120 const auto reference_cell = ReferenceCells::get_simplex<dim>();
121
122 // Piecewise constants are a special case: use a support point at the
123 // centroid and only the centroid
124 if (degree == 0)
125 {
126 unit_points.emplace_back(reference_cell.template barycenter<dim>());
127 return unit_points;
128 }
129
130 // otherwise write everything as linear combinations of vertices
131 const auto dpo = get_dpo_vector_fe_p(dim, degree);
132 Assert(dpo.size() == dim + 1, ExcInternalError());
133 Assert(dpo[0] == 1, ExcNotImplemented());
134 // vertices:
135 for (const unsigned int d : reference_cell.vertex_indices())
136 unit_points.push_back(reference_cell.template vertex<dim>(d));
137 // lines:
138 for (const unsigned int l : reference_cell.line_indices())
139 {
140 const Point<dim> p0 =
141 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
142 const Point<dim> p1 =
143 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
144 for (unsigned int p = 0; p < dpo[1]; ++p)
145 unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
146 (double(p + 1) / (dpo[1] + 1)) * p1);
147 }
148 // quads:
149 if (dim > 1 && dpo[2] > 0)
150 {
151 Assert(dpo[2] == 1, ExcNotImplemented());
152 if (dim == 2)
153 unit_points.push_back(reference_cell.template barycenter<dim>());
154 if (dim == 3)
155 for (const unsigned int f : reference_cell.face_indices())
156 unit_points.push_back(face_midpoint<dim>(f));
157 }
158
159 return unit_points;
160 }
161
162 template <>
163 std::vector<Point<0>>
164 unit_support_points_fe_p(const unsigned int /*degree*/)
165 {
166 return {Point<0>()};
167 }
168
173 template <int dim>
174 std::vector<std::vector<Point<dim - 1>>>
175 unit_face_support_points_fe_p(
176 const unsigned int degree,
177 typename FiniteElementData<dim>::Conformity conformity)
178 {
179 // Discontinuous elements don't have face support points
181 return {};
182
183 // this concept doesn't exist in 1d so just return an empty vector
184 if (dim == 1)
185 return {};
186
187 std::vector<std::vector<Point<dim - 1>>> unit_face_points;
188
189 // all faces have the same support points
190 for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
191 {
192 (void)face_n;
193 unit_face_points.emplace_back(
194 unit_support_points_fe_p<dim - 1>(degree));
195 }
196
197 return unit_face_points;
198 }
199
205 template <int dim>
207 constraints_fe_p(const unsigned int /*degree*/)
208 {
209 // no constraints in 1d
210 // constraints in 3d not implemented yet
211 return FullMatrix<double>();
212 }
213
214 template <>
216 constraints_fe_p<2>(const unsigned int degree)
217 {
218 constexpr int dim = 2;
219
220 // the following implements the 2d case
221 // (the 3d case is not implemented yet)
222 //
223 // consult FE_Q_Base::Implementation::initialize_constraints()
224 // for more information
225
226 std::vector<Point<dim - 1>> constraint_points;
227 // midpoint
228 constraint_points.emplace_back(0.5);
229 // subface 0
230 for (unsigned int i = 1; i < degree; ++i)
231 constraint_points.push_back(
233 Point<dim - 1>(i / double(degree)), 0));
234 // subface 1
235 for (unsigned int i = 1; i < degree; ++i)
236 constraint_points.push_back(
238 Point<dim - 1>(i / double(degree)), 1));
239
240 // Now construct relation between destination (child) and source (mother)
241 // dofs.
242
243 const unsigned int n_dofs_constrained = constraint_points.size();
244 unsigned int n_dofs_per_face = degree + 1;
245 FullMatrix<double> interface_constraints(n_dofs_constrained,
246 n_dofs_per_face);
247
248 const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
249
250 for (unsigned int i = 0; i < n_dofs_constrained; ++i)
251 for (unsigned int j = 0; j < n_dofs_per_face; ++j)
252 {
253 interface_constraints(i, j) =
254 poly.compute_value(j, constraint_points[i]);
255
256 // if the value is small up to round-off, then simply set it to zero
257 // to avoid unwanted fill-in of the constraint matrices (which would
258 // then increase the number of other DoFs a constrained DoF would
259 // couple to)
260 if (std::fabs(interface_constraints(i, j)) < 1e-13)
261 interface_constraints(i, j) = 0;
262 }
263 return interface_constraints;
264 }
265
266
267
272 std::vector<unsigned int>
273 get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
274 {
275 // First treat the case of piecewise constant elements:
276 if (degree == 0)
277 {
278 std::vector<unsigned int> dpo(dim + 1, 0U);
279 dpo[dim] = 1;
280 return dpo;
281 }
282 else
283 {
284 // This element has the same degrees of freedom as the continuous one,
285 // but they are all counted for the interior of the cell because
286 // it is continuous. Rather than hard-code how many DoFs the element
287 // has, we just get the numbers from the continuous case and add them
288 // up
289 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
290
291 switch (dim)
292 {
293 case 1:
294 return {0U,
295 ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
296 continuous_dpo[dim]};
297
298 case 2:
299 return {0U,
300 0U,
302 continuous_dpo[0] +
303 ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
304 continuous_dpo[dim]};
305
306 case 3:
307 return {
308 0U,
309 0U,
310 0U,
311 ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
312 ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
313 ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
314 continuous_dpo[dim]};
315 }
316
318 return {};
319 }
320 }
321} // namespace
322
323
324
325template <int dim, int spacedim>
327 const BarycentricPolynomials<dim> polynomials,
328 const FiniteElementData<dim> &fe_data,
329 const bool prolongation_is_additive,
330 const std::vector<Point<dim>> &unit_support_points,
331 const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
332 const FullMatrix<double> &interface_constraints)
333 : ::FE_Poly<dim, spacedim>(
334 polynomials,
335 fe_data,
336 std::vector<bool>(fe_data.dofs_per_cell, prolongation_is_additive),
337 std::vector<ComponentMask>(fe_data.dofs_per_cell,
338 ComponentMask(std::vector<bool>(1, true))))
339{
342 this->interface_constraints = interface_constraints;
343}
344
345
346
347template <int dim, int spacedim>
348std::pair<Table<2, bool>, std::vector<unsigned int>>
350{
351 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
352 constant_modes.fill(true);
353 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
354 constant_modes, std::vector<unsigned int>(1, 0));
355}
356
357
358
359template <int dim, int spacedim>
360const FullMatrix<double> &
362 const unsigned int child,
363 const RefinementCase<dim> &refinement_case) const
364{
365 if (dim == 3)
366 Assert(RefinementCase<dim>(refinement_case) ==
368 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
369 RefinementCase<dim>(refinement_case) ==
371 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
372 RefinementCase<dim>(refinement_case) ==
374 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
376 else
377 Assert(refinement_case ==
380 AssertDimension(dim, spacedim);
381
382 // initialization upon first request
383 if (this->prolongation[refinement_case - 1][child].n() == 0)
384 {
385 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
386
387 // if matrix got updated while waiting for the lock
388 if (this->prolongation[refinement_case - 1][child].n() ==
389 this->n_dofs_per_cell())
390 return this->prolongation[refinement_case - 1][child];
391
392 // now do the work. need to get a non-const version of data in order to
393 // be able to modify them inside a const function
394 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
395
396 if (dim == 2)
397 {
398 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
400 isotropic_matrices.back().resize(
401 this->reference_cell().n_children(
402 RefinementCase<dim>(refinement_case)),
403 FullMatrix<double>(this->n_dofs_per_cell(),
404 this->n_dofs_per_cell()));
405
406 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
407
408 this_nonconst.prolongation[refinement_case - 1] =
409 std::move(isotropic_matrices.back());
410 }
411 else if (dim == 3)
412 {
413 std::vector<std::vector<FullMatrix<double>>> matrices(
414 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49),
415 std::vector<FullMatrix<double>>(
416 this->reference_cell().n_children(
417 RefinementCase<dim>(refinement_case)),
418 FullMatrix<double>(this->n_dofs_per_cell(),
419 this->n_dofs_per_cell())));
420 FETools::compute_embedding_matrices(*this, matrices, true);
421 for (unsigned int refinement_direction = static_cast<unsigned int>(
423 refinement_direction <=
424 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49);
425 refinement_direction++)
426 this_nonconst.prolongation[refinement_direction - 1] =
427 std::move(matrices[refinement_direction - 1]);
428 }
429 else
431 }
432
433 // finally return the matrix
434 return this->prolongation[refinement_case - 1][child];
435}
436
437
438
439template <int dim, int spacedim>
440const FullMatrix<double> &
442 const unsigned int child,
443 const RefinementCase<dim> &refinement_case) const
444{
445 if (dim == 3)
446 Assert(RefinementCase<dim>(refinement_case) ==
448 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
449 RefinementCase<dim>(refinement_case) ==
451 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
452 RefinementCase<dim>(refinement_case) ==
454 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
456 else
459 AssertDimension(dim, spacedim);
460
461 // initialization upon first request
462 if (this->restriction[refinement_case - 1][child].n() == 0)
463 {
464 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
465
466 // if matrix got updated while waiting for the lock
467 if (this->restriction[refinement_case - 1][child].n() ==
468 this->n_dofs_per_cell())
469 return this->restriction[refinement_case - 1][child];
470
471 // get the restriction matrix
472 // Refine a unit cell. As the parent cell is a unit
473 // cell, the reference cell of the children equals the parent, i.e. they
474 // have the support points at the same locations. So we just have to check
475 // if a support point of the parent is one of the interpolation points of
476 // the child. If this is not the case we find the interpolation of the
477 // point.
478
479 const double eps = 1e-12;
480 FullMatrix<double> restriction_mat(this->n_dofs_per_cell(),
481 this->n_dofs_per_cell());
482
483 // first get all support points on the reference cell
484 const std::vector<Point<dim>> unit_support_points =
485 this->get_unit_support_points();
486
487 // now create children on the reference cell
489 GridGenerator::reference_cell(tria, this->reference_cell());
490 tria.begin_active()->set_refine_flag(
492 if (dim == 3)
493 tria.begin_active()->set_refine_choice(refinement_case);
495
496 const auto &child_cell = tria.begin(0)->child(child);
497
498 // iterate over all support points and transform them to the unit cell of
499 // the child
500 for (unsigned int i = 0; i < unit_support_points.size(); i++)
501 {
502 std::vector<Point<dim>> transformed_point(1);
503 const std::vector<Point<spacedim>> unit_support_point = {
504 dim == 2 ? Point<spacedim>(unit_support_points[i][0],
505 unit_support_points[i][1]) :
506 Point<spacedim>(unit_support_points[i][0],
507 unit_support_points[i][1],
508 unit_support_points[i][2])};
509 this->reference_cell()
510 .template get_default_linear_mapping<dim, spacedim>()
511 .transform_points_real_to_unit_cell(
512 child_cell,
513 make_array_view(unit_support_point),
514 make_array_view(transformed_point));
515
516 // if point is inside the unit cell iterate over all shape functions
517 if (this->reference_cell().contains_point(transformed_point[0], eps))
518 for (unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
519 restriction_mat[i][j] =
520 this->shape_value(j, transformed_point[0]);
521 }
522 if constexpr (running_in_debug_mode())
523 {
524 for (unsigned int i = 0; i < this->n_dofs_per_cell(); i++)
525 {
526 double sum = 0.;
527
528 for (unsigned int j = 0; j < this->n_dofs_per_cell(); j++)
529 sum += restriction_mat[i][j];
530
531 Assert(std::fabs(sum - 1) < eps || std::fabs(sum) < eps,
533 "The entries in a row of the local "
534 "restriction matrix do not add to zero or one. "
535 "This typically indicates that the "
536 "polynomial interpolation is "
537 "ill-conditioned such that round-off "
538 "prevents the sum to be one."));
539 }
540 }
541
542 // Remove small entries from the matrix
543 for (unsigned int i = 0; i < restriction_mat.m(); ++i)
544 for (unsigned int j = 0; j < restriction_mat.n(); ++j)
545 {
546 if (std::fabs(restriction_mat(i, j)) < eps)
547 restriction_mat(i, j) = 0.;
548 if (std::fabs(restriction_mat(i, j) - 1) < eps)
549 restriction_mat(i, j) = 1.;
550 }
551
552 const_cast<FullMatrix<double> &>(
553 this->restriction[refinement_case - 1][child]) =
554 std::move(restriction_mat);
555 }
556
557 // finally return the matrix
558 return this->restriction[refinement_case - 1][child];
559}
560
561
562
563template <int dim, int spacedim>
564void
566 const FiniteElement<dim, spacedim> &source_fe,
567 FullMatrix<double> &interpolation_matrix,
568 const unsigned int face_no) const
569{
570 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
571 ExcDimensionMismatch(interpolation_matrix.m(),
572 source_fe.n_dofs_per_face(face_no)));
573
574 // see if source is a P or Q element
575 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
576 nullptr) ||
577 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
578 {
579 const Quadrature<dim - 1> quad_face_support(
580 source_fe.get_unit_face_support_points(face_no));
581
582 const double eps = 2e-13 * this->degree * (dim - 1);
583
584 const std::vector<Point<dim>> face_quadrature_points =
585 QProjector<dim>::project_to_face(this->reference_cell(),
586 quad_face_support,
587 face_no,
589 .get_points();
590
591 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
592 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
593 {
594 double matrix_entry =
595 this->shape_value(this->face_to_cell_index(j, 0),
596 face_quadrature_points[i]);
597
598 // Correct the interpolated value. I.e. if it is close to 1 or
599 // 0, make it exactly 1 or 0. Unfortunately, this is required to
600 // avoid problems with higher order elements.
601 if (std::fabs(matrix_entry - 1.0) < eps)
602 matrix_entry = 1.0;
603 if (std::fabs(matrix_entry) < eps)
604 matrix_entry = 0.0;
605
606 interpolation_matrix(i, j) = matrix_entry;
607 }
608
609 if constexpr (running_in_debug_mode())
610 {
611 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
612 {
613 double sum = 0.;
614
615 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
616 sum += interpolation_matrix(j, i);
617
618 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
619 }
620 }
621 }
622 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
623 {
624 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
625 }
626 else
628 false,
629 (typename FiniteElement<dim,
630 spacedim>::ExcInterpolationNotImplemented()));
631}
632
633
634
635template <int dim, int spacedim>
636void
638 const FiniteElement<dim, spacedim> &source_fe,
639 const unsigned int subface,
640 FullMatrix<double> &interpolation_matrix,
641 const unsigned int face_no) const
642{
643 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
644 ExcDimensionMismatch(interpolation_matrix.m(),
645 source_fe.n_dofs_per_face(face_no)));
646
647 // see if source is a P or Q element
648 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
649 nullptr) ||
650 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
651 {
652 const Quadrature<dim - 1> quad_face_support(
653 source_fe.get_unit_face_support_points(face_no));
654
655 const double eps = 2e-13 * this->degree * (dim - 1);
656
657 const Quadrature<dim> subface_quadrature =
659 this->reference_cell(),
660 quad_face_support,
661 face_no,
662 subface,
665
666 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
667 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
668 {
669 double matrix_entry =
670 this->shape_value(this->face_to_cell_index(j, 0),
671 subface_quadrature.point(i));
672
673 // Correct the interpolated value. I.e. if it is close to 1 or
674 // 0, make it exactly 1 or 0. Unfortunately, this is required to
675 // avoid problems with higher order elements.
676 if (std::fabs(matrix_entry - 1.0) < eps)
677 matrix_entry = 1.0;
678 if (std::fabs(matrix_entry) < eps)
679 matrix_entry = 0.0;
680
681 interpolation_matrix(i, j) = matrix_entry;
682 }
683
684 if constexpr (running_in_debug_mode())
685 {
686 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
687 {
688 double sum = 0.;
689
690 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
691 sum += interpolation_matrix(j, i);
692
693 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
694 }
695 }
696 }
697 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
698 {
699 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
700 }
701 else
703 false,
704 (typename FiniteElement<dim,
705 spacedim>::ExcInterpolationNotImplemented()));
706}
707
708
709
710template <int dim, int spacedim>
711bool
716
717
718
719template <int dim, int spacedim>
720void
723 const std::vector<Vector<double>> &support_point_values,
724 std::vector<double> &nodal_values) const
725{
726 AssertDimension(support_point_values.size(),
727 this->get_unit_support_points().size());
728 AssertDimension(support_point_values.size(), nodal_values.size());
729 AssertDimension(this->dofs_per_cell, nodal_values.size());
730
731 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
732 {
733 AssertDimension(support_point_values[i].size(), 1);
734
735 nodal_values[i] = support_point_values[i](0);
736 }
737}
738
739
740
741template <int dim, int spacedim>
743 : FE_SimplexPoly<dim, spacedim>(
744 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
745 FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
746 ReferenceCells::get_simplex<dim>(),
747 1,
748 degree,
749 FiniteElementData<dim>::H1),
750 false,
751 unit_support_points_fe_p<dim>(degree),
752 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
753 constraints_fe_p<dim>(degree))
754{
755 if (degree > 2)
756 for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
758 this->n_dofs_per_line() - 1 - i - i;
759 // We do not support multiple DoFs per quad yet
761}
762
763
764
765template <int dim, int spacedim>
766std::unique_ptr<FiniteElement<dim, spacedim>>
768{
769 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
770}
771
772
773
774template <int dim, int spacedim>
775std::string
777{
778 std::ostringstream namebuf;
779 namebuf << "FE_SimplexP<" << Utilities::dim_string(dim, spacedim) << ">("
780 << this->degree << ")";
781
782 return namebuf.str();
783}
784
785
786
787template <int dim, int spacedim>
790 const FiniteElement<dim, spacedim> &fe_other,
791 const unsigned int codim) const
792{
793 Assert(codim <= dim, ExcImpossibleInDim(dim));
794
795 // vertex/line/face domination
796 // (if fe_other is derived from FE_SimplexDGP)
797 // ------------------------------------
798 if (codim > 0)
799 if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
800 nullptr)
801 // there are no requirements between continuous and discontinuous
802 // elements
804
805 // vertex/line/face domination
806 // (if fe_other is not derived from FE_SimplexDGP)
807 // & cell domination
808 // ----------------------------------------
809 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
810 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
811 {
812 if (this->degree < fe_p_other->degree)
814 else if (this->degree == fe_p_other->degree)
816 else
818 }
819 else if (const FE_Q<dim, spacedim> *fe_q_other =
820 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
821 {
822 if (this->degree < fe_q_other->degree)
824 else if (this->degree == fe_q_other->degree)
826 else
828 }
829 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
830 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
831 {
832 if (fe_nothing->is_dominating())
834 else
835 // the FE_Nothing has no degrees of freedom and it is typically used
836 // in a context where we don't require any continuity along the
837 // interface
839 }
840
843}
844
845
846
847template <int dim, int spacedim>
848std::vector<std::pair<unsigned int, unsigned int>>
850 const FiniteElement<dim, spacedim> &fe_other) const
851{
852 AssertDimension(dim, 2);
853
854 if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
855 {
856 // there should be exactly one single DoF of each FE at a vertex, and
857 // they should have identical value
858 return {{0U, 0U}};
859 }
860 else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
861 {
862 // there should be exactly one single DoF of each FE at a vertex, and
863 // they should have identical value
864 return {{0U, 0U}};
865 }
866 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
867 {
868 // the FE_Nothing has no degrees of freedom, so there are no
869 // equivalencies to be recorded
870 return {};
871 }
872 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
873 {
874 // if the other element has no DoFs on faces at all,
875 // then it would be impossible to enforce any kind of
876 // continuity even if we knew exactly what kind of element
877 // we have -- simply because the other element declares
878 // that it is discontinuous because it has no DoFs on
879 // its faces. in that case, just state that we have no
880 // constraints to declare
881 return {};
882 }
883 else
884 {
886 return {};
887 }
888}
889
890
891
892template <int dim, int spacedim>
893std::vector<std::pair<unsigned int, unsigned int>>
895 const FiniteElement<dim, spacedim> &fe_other) const
896{
897 AssertDimension(dim, 2);
898
899 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
900 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
901 {
902 // dofs are located along lines, so two dofs are identical if they are
903 // located at identical positions.
904 // Therefore, read the points in unit_support_points for the
905 // first coordinate direction. For FE_SimplexP, they are currently
906 // hard-coded and we iterate over points on the first line which begin
907 // after the 3 vertex points in the complete list of unit support points
908
909 std::vector<std::pair<unsigned int, unsigned int>> identities;
910
911 for (unsigned int i = 0; i < this->degree - 1; ++i)
912 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
913 if (std::fabs(this->unit_support_points[i + 3][0] -
914 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
915 identities.emplace_back(i, j);
916 else
917 {
918 // If nodes are not located in the same place, we have to
919 // interpolate. This is then not handled through the
920 // current function, but via interpolation matrices that
921 // result in constraints, rather than identities. Since
922 // that happens in a different function, there is nothing
923 // for us to do here.
924 }
925
926 return identities;
927 }
928 else if (const FE_Q<dim, spacedim> *fe_q_other =
929 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
930 {
931 // dofs are located along lines, so two dofs are identical if they are
932 // located at identical positions. if we had only equidistant points, we
933 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
934 // might have other support points (e.g. Gauss-Lobatto
935 // points). Therefore, read the points in unit_support_points for the
936 // first coordinate direction. For FE_Q, we take the lexicographic
937 // ordering of the line support points in the first direction (i.e.,
938 // x-direction), which we access between index 1 and p-1 (index 0 and p
939 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
940 // iterate over points on the first line which begin after the 3 vertex
941 // points in the complete list of unit support points
942
943 const std::vector<unsigned int> &index_map_inverse_q_other =
944 fe_q_other->get_poly_space_numbering_inverse();
945
946 std::vector<std::pair<unsigned int, unsigned int>> identities;
947
948 for (unsigned int i = 0; i < this->degree - 1; ++i)
949 for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
950 if (std::fabs(this->unit_support_points[i + 3][0] -
951 fe_q_other->get_unit_support_points()
952 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
953 identities.emplace_back(i, j);
954 else
955 {
956 // If nodes are not located in the same place, we have to
957 // interpolate. This will then also
958 // capture the case where the FE_Q has a different polynomial
959 // degree than the current element. In either case, the resulting
960 // constraints are computed elsewhere, rather than via the
961 // identities this function returns: Since
962 // that happens in a different function, there is nothing
963 // for us to do here.
964 }
965
966 return identities;
967 }
968 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
969 {
970 // The FE_Nothing has no degrees of freedom, so there are no
971 // equivalencies to be recorded. (If the FE_Nothing is dominating,
972 // then this will also leads to constraints, but we are not concerned
973 // with this here.)
974 return {};
975 }
976 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
977 {
978 // if the other element has no elements on faces at all,
979 // then it would be impossible to enforce any kind of
980 // continuity even if we knew exactly what kind of element
981 // we have -- simply because the other element declares
982 // that it is discontinuous because it has no DoFs on
983 // its faces. in that case, just state that we have no
984 // constraints to declare
985 return {};
986 }
987 else
988 {
990 return {};
991 }
992}
993
994
995
996template <int dim, int spacedim>
998 : FE_SimplexPoly<dim, spacedim>(
999 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
1000 FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
1001 ReferenceCells::get_simplex<dim>(),
1002 1,
1003 degree,
1004 FiniteElementData<dim>::L2),
1005 true,
1006 unit_support_points_fe_p<dim>(degree),
1007 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::L2),
1008 constraints_fe_p<dim>(degree))
1009{}
1010
1011
1012
1013template <int dim, int spacedim>
1014std::unique_ptr<FiniteElement<dim, spacedim>>
1016{
1017 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
1018}
1019
1020
1021
1022template <int dim, int spacedim>
1023std::string
1025{
1026 std::ostringstream namebuf;
1027 namebuf << "FE_SimplexDGP<" << Utilities::dim_string(dim, spacedim) << ">("
1028 << this->degree << ")";
1029
1030 return namebuf.str();
1031}
1032
1033
1034template <int dim, int spacedim>
1037 const FiniteElement<dim, spacedim> &fe_other,
1038 const unsigned int codim) const
1039{
1040 Assert(codim <= dim, ExcImpossibleInDim(dim));
1041
1042 // vertex/line/face domination
1043 // ---------------------------
1044 if (codim > 0)
1045 // this is a discontinuous element, so by definition there will
1046 // be no constraints wherever this element comes together with
1047 // any other kind of element
1049
1050 // cell domination
1051 // ---------------
1052 if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
1053 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
1054 {
1055 if (this->degree < fe_dgp_other->degree)
1057 else if (this->degree == fe_dgp_other->degree)
1059 else
1061 }
1062 else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
1063 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
1064 {
1065 if (this->degree < fe_dgq_other->degree)
1067 else if (this->degree == fe_dgq_other->degree)
1069 else
1071 }
1072 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
1073 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
1074 {
1075 if (fe_nothing->is_dominating())
1077 else
1078 // the FE_Nothing has no degrees of freedom and it is typically used
1079 // in a context where we don't require any continuity along the
1080 // interface
1082 }
1083
1086}
1087
1088
1089
1090template <int dim, int spacedim>
1091std::vector<std::pair<unsigned int, unsigned int>>
1093 const FiniteElement<dim, spacedim> &fe_other) const
1094{
1095 (void)fe_other;
1096
1097 return {};
1098}
1099
1100
1101
1102template <int dim, int spacedim>
1103std::vector<std::pair<unsigned int, unsigned int>>
1105 const FiniteElement<dim, spacedim> &fe_other) const
1106{
1107 (void)fe_other;
1108
1109 return {};
1110}
1111
1112
1113
1114template <int dim, int spacedim>
1115const FullMatrix<double> &
1117 const unsigned int child,
1118 const RefinementCase<dim> &refinement_case) const
1119{
1120 if (dim == 3)
1121 Assert(RefinementCase<dim>(refinement_case) ==
1123 static_cast<char>(IsotropicRefinementChoice::cut_tet_68)) ||
1124 RefinementCase<dim>(refinement_case) ==
1126 static_cast<char>(IsotropicRefinementChoice::cut_tet_57)) ||
1127 RefinementCase<dim>(refinement_case) ==
1129 static_cast<char>(IsotropicRefinementChoice::cut_tet_49)),
1131 else
1134 AssertDimension(dim, spacedim);
1135
1136 // initialization upon first request
1137 if (this->restriction[refinement_case - 1][child].n() == 0)
1138 {
1139 std::lock_guard<std::mutex> lock(this->restriction_matrix_mutex);
1140
1141 // if matrix got updated while waiting for the lock
1142 if (this->restriction[refinement_case - 1][child].n() ==
1143 this->n_dofs_per_cell())
1144 return this->restriction[refinement_case - 1][child];
1145
1146 // now do the work. need to get a non-const version of data in order to
1147 // be able to modify them inside a const function
1148 auto &this_nonconst = const_cast<FE_SimplexDGP<dim, spacedim> &>(*this);
1149
1150 if (dim == 2)
1151 {
1152 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
1154 isotropic_matrices.back().resize(
1155 this->reference_cell().n_children(
1156 RefinementCase<dim>(refinement_case)),
1157 FullMatrix<double>(this->n_dofs_per_cell(),
1158 this->n_dofs_per_cell()));
1159
1160 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
1161
1162 this_nonconst.restriction[refinement_case - 1] =
1163 std::move(isotropic_matrices.back());
1164 }
1165 else if (dim == 3)
1166 {
1167 std::vector<std::vector<FullMatrix<double>>> matrices(
1168 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49),
1169 std::vector<FullMatrix<double>>(
1170 this->reference_cell().n_children(
1171 RefinementCase<dim>(refinement_case)),
1172 FullMatrix<double>(this->n_dofs_per_cell(),
1173 this->n_dofs_per_cell())));
1174 FETools::compute_projection_matrices(*this, matrices, true);
1175 for (unsigned int refinement_direction = static_cast<unsigned int>(
1177 refinement_direction <=
1178 static_cast<unsigned int>(IsotropicRefinementChoice::cut_tet_49);
1179 refinement_direction++)
1180 this_nonconst.restriction[refinement_direction - 1] =
1181 std::move(matrices[refinement_direction - 1]);
1182 }
1183 else
1185 }
1186
1187 // finally return the matrix
1188 return this->restriction[refinement_case - 1][child];
1189}
1190
1191// explicit instantiations
1192#include "fe/fe_simplex_p.inst"
1193
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
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:2568
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:2610
std::vector< Point< dim > > unit_support_points
Definition fe.h:2561
FullMatrix< double > interface_constraints
Definition fe.h:2549
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:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#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)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#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 const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:555
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:346
STL namespace.