Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_simplex_p.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
16#include <deal.II/base/config.h>
17
20
21#include <deal.II/fe/fe_dgq.h>
23#include <deal.II/fe/fe_q.h>
25#include <deal.II/fe/fe_tools.h>
26
28
29namespace
30{
35 std::vector<unsigned int>
36 get_dpo_vector_fe_p(const unsigned int dim, const unsigned int degree)
37 {
38 switch (dim)
39 {
40 case 1:
41 switch (degree)
42 {
43 case 1:
44 return {1, 0};
45 case 2:
46 return {1, 1};
47 default:
48 Assert(false, ExcNotImplemented());
49 }
50 case 2:
51 switch (degree)
52 {
53 case 1:
54 return {1, 0, 0};
55 case 2:
56 return {1, 1, 0};
57 default:
58 Assert(false, ExcNotImplemented());
59 }
60 case 3:
61 switch (degree)
62 {
63 case 1:
64 return {1, 0, 0, 0};
65 case 2:
66 return {1, 1, 0, 0};
67 default:
68 Assert(false, ExcNotImplemented());
69 }
70 }
71
72 Assert(false, ExcNotImplemented());
73 return {};
74 }
75
76
77
82 template <int dim>
83 std::vector<Point<dim>>
84 unit_support_points_fe_p(const unsigned int degree)
85 {
86 std::vector<Point<dim>> unit_points;
87 // If we do dim - 1 we can get here in dim = 0:
88 if (dim == 0)
89 {
90 unit_points.emplace_back();
91 return unit_points;
92 }
93
94 const auto reference_cell = ReferenceCells::get_simplex<dim>();
95 // Piecewise constants are a special case: use a support point at the
96 // centroid and only the centroid
97 if (degree == 0)
98 {
99 unit_points.emplace_back(reference_cell.template barycenter<dim>());
100 return unit_points;
102
103 Assert(degree <= 2, ExcNotImplemented());
104 for (const auto &vertex_no : reference_cell.vertex_indices())
105 unit_points.emplace_back(reference_cell.template vertex<dim>(vertex_no));
106
107 if (degree == 2)
108 for (const auto &line_no : reference_cell.line_indices())
109 {
110 const auto v0 = reference_cell.template vertex<dim>(
111 reference_cell.line_to_cell_vertices(line_no, 0));
112 const auto v1 = reference_cell.template vertex<dim>(
113 reference_cell.line_to_cell_vertices(line_no, 1));
114 unit_points.emplace_back((v0 + v1) / 2.0);
115 }
116
117 return unit_points;
118 }
119
124 template <int dim>
125 std::vector<std::vector<Point<dim - 1>>>
126 unit_face_support_points_fe_p(
127 const unsigned int degree,
128 typename FiniteElementData<dim>::Conformity conformity)
129 {
130 // Discontinuous elements don't have face support points
132 return {};
133
134 // this concept doesn't exist in 1d so just return an empty vector
135 if (dim == 1)
136 return {};
137
138 std::vector<std::vector<Point<dim - 1>>> unit_face_points;
139
140 // all faces have the same support points
141 for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
142 {
143 (void)face_n;
144 unit_face_points.emplace_back(
145 unit_support_points_fe_p<dim - 1>(degree));
146 }
147
148 return unit_face_points;
149 }
150
156 template <int dim>
158 constraints_fe_p(const unsigned int /*degree*/)
159 {
160 // no constraints in 1d
161 // constraints in 3d not implemented yet
162 return FullMatrix<double>();
163 }
164
165 template <>
167 constraints_fe_p<2>(const unsigned int degree)
168 {
169 const unsigned int dim = 2;
170
171 Assert(degree <= 2, ExcNotImplemented());
172
173 // the following implements the 2d case
174 // (the 3d case is not implemented yet)
175 //
176 // consult FE_Q_Base::Implementation::initialize_constraints()
177 // for more information
178
179 std::vector<Point<dim - 1>> constraint_points;
180 // midpoint
181 constraint_points.emplace_back(0.5);
182 if (degree == 2)
183 {
184 // midpoint on subface 0
185 constraint_points.emplace_back(0.25);
186 // midpoint on subface 1
187 constraint_points.emplace_back(0.75);
188 }
189
190 // Now construct relation between destination (child) and source (mother)
191 // dofs.
192
193 const unsigned int n_dofs_constrained = constraint_points.size();
194 unsigned int n_dofs_per_face = degree + 1;
195 FullMatrix<double> interface_constraints(n_dofs_constrained,
196 n_dofs_per_face);
197
198 const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
199
200 for (unsigned int i = 0; i < n_dofs_constrained; ++i)
201 for (unsigned int j = 0; j < n_dofs_per_face; ++j)
202 {
203 interface_constraints(i, j) =
204 poly.compute_value(j, constraint_points[i]);
205
206 // if the value is small up to round-off, then simply set it to zero
207 // to avoid unwanted fill-in of the constraint matrices (which would
208 // then increase the number of other DoFs a constrained DoF would
209 // couple to)
210 if (std::fabs(interface_constraints(i, j)) < 1e-13)
211 interface_constraints(i, j) = 0;
212 }
213 return interface_constraints;
214 }
215
216
217
222 std::vector<unsigned int>
223 get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
224 {
225 // First treat the case of piecewise constant elements:
226 if (degree == 0)
227 {
228 std::vector<unsigned int> dpo(dim + 1, 0U);
229 dpo[dim] = 1;
230 return dpo;
231 }
232 else
233 {
234 // This element has the same degrees of freedom as the continuous one,
235 // but they are all counted for the interior of the cell because
236 // it is continuous. Rather than hard-code how many DoFs the element
237 // has, we just get the numbers from the continuous case and add them
238 // up
239 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
240
241 switch (dim)
242 {
243 case 1:
244 return {0U,
245 ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
246 continuous_dpo[dim]};
247
248 case 2:
249 return {0U,
250 0U,
252 continuous_dpo[0] +
253 ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
254 continuous_dpo[dim]};
255
256 case 3:
257 return {
258 0U,
259 0U,
260 0U,
261 ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
262 ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
263 ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
264 continuous_dpo[dim]};
265 }
266
267 Assert(false, ExcNotImplemented());
268 return {};
269 }
270 }
271} // namespace
272
273
274
275template <int dim, int spacedim>
277 const BarycentricPolynomials<dim> polynomials,
278 const FiniteElementData<dim> & fe_data,
279 const std::vector<Point<dim>> & unit_support_points,
280 const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
281 const FullMatrix<double> & interface_constraints)
282 : ::FE_Poly<dim, spacedim>(
283 polynomials,
284 fe_data,
285 std::vector<bool>(fe_data.dofs_per_cell),
286 std::vector<ComponentMask>(fe_data.dofs_per_cell,
287 std::vector<bool>(1, true)))
288{
291 this->interface_constraints = interface_constraints;
292}
293
294
295
296template <int dim, int spacedim>
297std::pair<Table<2, bool>, std::vector<unsigned int>>
299{
300 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
301 constant_modes.fill(true);
302 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
303 constant_modes, std::vector<unsigned int>(1, 0));
304}
305
306
307
308template <int dim, int spacedim>
309const FullMatrix<double> &
311 const unsigned int child,
312 const RefinementCase<dim> &refinement_case) const
313{
316 AssertDimension(dim, spacedim);
317
318 // initialization upon first request
319 if (this->prolongation[refinement_case - 1][child].n() == 0)
320 {
321 std::lock_guard<std::mutex> lock(this->mutex);
322
323 // if matrix got updated while waiting for the lock
324 if (this->prolongation[refinement_case - 1][child].n() ==
325 this->n_dofs_per_cell())
326 return this->prolongation[refinement_case - 1][child];
327
328 // now do the work. need to get a non-const version of data in order to
329 // be able to modify them inside a const function
330 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
331
332 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
334 isotropic_matrices.back().resize(
336 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
337
338 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
339
340 this_nonconst.prolongation[refinement_case - 1].swap(
341 isotropic_matrices.back());
342 }
343
344 // finally return the matrix
345 return this->prolongation[refinement_case - 1][child];
346}
347
348
349
350template <int dim, int spacedim>
351const FullMatrix<double> &
353 const unsigned int child,
354 const RefinementCase<dim> &refinement_case) const
355{
358 AssertDimension(dim, spacedim);
359
360 // initialization upon first request
361 if (this->restriction[refinement_case - 1][child].n() == 0)
362 {
363 std::lock_guard<std::mutex> lock(this->mutex);
364
365 // if matrix got updated while waiting for the lock
366 if (this->restriction[refinement_case - 1][child].n() ==
367 this->n_dofs_per_cell())
368 return this->restriction[refinement_case - 1][child];
369
370 // now do the work. need to get a non-const version of data in order to
371 // be able to modify them inside a const function
372 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
373
374 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
376 isotropic_matrices.back().resize(
378 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
379
380 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
381
382 this_nonconst.restriction[refinement_case - 1].swap(
383 isotropic_matrices.back());
384 }
385
386 // finally return the matrix
387 return this->restriction[refinement_case - 1][child];
388}
389
390
391
392template <int dim, int spacedim>
393void
395 const FiniteElement<dim, spacedim> &source_fe,
396 FullMatrix<double> & interpolation_matrix,
397 const unsigned int face_no) const
398{
399 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
400 ExcDimensionMismatch(interpolation_matrix.m(),
401 source_fe.n_dofs_per_face(face_no)));
402
403 // see if source is a P or Q element
404 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
405 nullptr) ||
406 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
407 {
408 const Quadrature<dim - 1> quad_face_support(
409 source_fe.get_unit_face_support_points(face_no));
410
411 const double eps = 2e-13 * this->degree * (dim - 1);
412
413 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
414 QProjector<dim>::project_to_face(this->reference_cell(),
415 quad_face_support,
416 face_no,
417 face_quadrature_points);
418
419 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
420 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
421 {
422 double matrix_entry =
423 this->shape_value(this->face_to_cell_index(j, 0),
424 face_quadrature_points[i]);
425
426 // Correct the interpolated value. I.e. if it is close to 1 or
427 // 0, make it exactly 1 or 0. Unfortunately, this is required to
428 // avoid problems with higher order elements.
429 if (std::fabs(matrix_entry - 1.0) < eps)
430 matrix_entry = 1.0;
431 if (std::fabs(matrix_entry) < eps)
432 matrix_entry = 0.0;
433
434 interpolation_matrix(i, j) = matrix_entry;
435 }
436
437#ifdef DEBUG
438 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
439 {
440 double sum = 0.;
441
442 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
443 sum += interpolation_matrix(j, i);
444
445 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
446 }
447#endif
448 }
449 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
450 {
451 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
452 }
453 else
455 false,
456 (typename FiniteElement<dim,
457 spacedim>::ExcInterpolationNotImplemented()));
458}
459
460
461
462template <int dim, int spacedim>
463void
465 const FiniteElement<dim, spacedim> &source_fe,
466 const unsigned int subface,
467 FullMatrix<double> & interpolation_matrix,
468 const unsigned int face_no) const
469{
470 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
471 ExcDimensionMismatch(interpolation_matrix.m(),
472 source_fe.n_dofs_per_face(face_no)));
473
474 // see if source is a P or Q element
475 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
476 nullptr) ||
477 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
478 {
479 const Quadrature<dim - 1> quad_face_support(
480 source_fe.get_unit_face_support_points(face_no));
481
482 const double eps = 2e-13 * this->degree * (dim - 1);
483
484 std::vector<Point<dim>> subface_quadrature_points(
485 quad_face_support.size());
486 QProjector<dim>::project_to_subface(this->reference_cell(),
487 quad_face_support,
488 face_no,
489 subface,
490 subface_quadrature_points);
491
492 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
493 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
494 {
495 double matrix_entry =
496 this->shape_value(this->face_to_cell_index(j, 0),
497 subface_quadrature_points[i]);
498
499 // Correct the interpolated value. I.e. if it is close to 1 or
500 // 0, make it exactly 1 or 0. Unfortunately, this is required to
501 // avoid problems with higher order elements.
502 if (std::fabs(matrix_entry - 1.0) < eps)
503 matrix_entry = 1.0;
504 if (std::fabs(matrix_entry) < eps)
505 matrix_entry = 0.0;
506
507 interpolation_matrix(i, j) = matrix_entry;
508 }
509
510#ifdef DEBUG
511 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
512 {
513 double sum = 0.;
514
515 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
516 sum += interpolation_matrix(j, i);
517
518 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
519 }
520#endif
521 }
522 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
523 {
524 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
525 }
526 else
528 false,
529 (typename FiniteElement<dim,
530 spacedim>::ExcInterpolationNotImplemented()));
531}
532
533
534
535template <int dim, int spacedim>
536bool
538{
539 return true;
540}
541
542
543
544template <int dim, int spacedim>
545void
548 const std::vector<Vector<double>> &support_point_values,
549 std::vector<double> & nodal_values) const
550{
551 AssertDimension(support_point_values.size(),
552 this->get_unit_support_points().size());
553 AssertDimension(support_point_values.size(), nodal_values.size());
554 AssertDimension(this->dofs_per_cell, nodal_values.size());
555
556 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
557 {
558 AssertDimension(support_point_values[i].size(), 1);
559
560 nodal_values[i] = support_point_values[i](0);
561 }
562}
563
564
565
566template <int dim, int spacedim>
568 : FE_SimplexPoly<dim, spacedim>(
569 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
570 FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
571 ReferenceCells::get_simplex<dim>(),
572 1,
573 degree,
574 FiniteElementData<dim>::H1),
575 unit_support_points_fe_p<dim>(degree),
576 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
577 constraints_fe_p<dim>(degree))
578{}
579
580
581
582template <int dim, int spacedim>
583std::unique_ptr<FiniteElement<dim, spacedim>>
585{
586 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
587}
588
589
590
591template <int dim, int spacedim>
592std::string
594{
595 std::ostringstream namebuf;
596 namebuf << "FE_SimplexP<" << dim << ">(" << this->degree << ")";
597
598 return namebuf.str();
599}
600
601
602
603template <int dim, int spacedim>
606 const FiniteElement<dim, spacedim> &fe_other,
607 const unsigned int codim) const
608{
609 Assert(codim <= dim, ExcImpossibleInDim(dim));
610
611 // vertex/line/face domination
612 // (if fe_other is derived from FE_SimplexDGP)
613 // ------------------------------------
614 if (codim > 0)
615 if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
616 nullptr)
617 // there are no requirements between continuous and discontinuous
618 // elements
620
621 // vertex/line/face domination
622 // (if fe_other is not derived from FE_SimplexDGP)
623 // & cell domination
624 // ----------------------------------------
625 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
626 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
627 {
628 if (this->degree < fe_p_other->degree)
630 else if (this->degree == fe_p_other->degree)
632 else
634 }
635 else if (const FE_Q<dim, spacedim> *fe_q_other =
636 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
637 {
638 if (this->degree < fe_q_other->degree)
640 else if (this->degree == fe_q_other->degree)
642 else
644 }
645 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
646 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
647 {
648 if (fe_nothing->is_dominating())
650 else
651 // the FE_Nothing has no degrees of freedom and it is typically used
652 // in a context where we don't require any continuity along the
653 // interface
655 }
656
657 Assert(false, ExcNotImplemented());
659}
660
661
662
663template <int dim, int spacedim>
664std::vector<std::pair<unsigned int, unsigned int>>
666 const FiniteElement<dim, spacedim> &fe_other) const
667{
668 AssertDimension(dim, 2);
669
670 if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
671 {
672 // there should be exactly one single DoF of each FE at a vertex, and
673 // they should have identical value
674 return {{0U, 0U}};
675 }
676 else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
677 {
678 // there should be exactly one single DoF of each FE at a vertex, and
679 // they should have identical value
680 return {{0U, 0U}};
681 }
682 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
683 {
684 // the FE_Nothing has no degrees of freedom, so there are no
685 // equivalencies to be recorded
686 return {};
687 }
688 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
689 {
690 // if the other element has no DoFs on faces at all,
691 // then it would be impossible to enforce any kind of
692 // continuity even if we knew exactly what kind of element
693 // we have -- simply because the other element declares
694 // that it is discontinuous because it has no DoFs on
695 // its faces. in that case, just state that we have no
696 // constraints to declare
697 return {};
698 }
699 else
700 {
701 Assert(false, ExcNotImplemented());
702 return {};
703 }
704}
705
706
707
708template <int dim, int spacedim>
709std::vector<std::pair<unsigned int, unsigned int>>
711 const FiniteElement<dim, spacedim> &fe_other) const
712{
713 AssertDimension(dim, 2);
714
715 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
716 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
717 {
718 // dofs are located along lines, so two dofs are identical if they are
719 // located at identical positions.
720 // Therefore, read the points in unit_support_points for the
721 // first coordinate direction. For FE_SimplexP, they are currently
722 // hard-coded and we iterate over points on the first line which begin
723 // after the 3 vertex points in the complete list of unit support points
724
725 std::vector<std::pair<unsigned int, unsigned int>> identities;
726
727 for (unsigned int i = 0; i < this->degree - 1; ++i)
728 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
729 if (std::fabs(this->unit_support_points[i + 3][0] -
730 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
731 identities.emplace_back(i, j);
732 else
733 {
734 // If nodes are not located in the same place, we have to
735 // interpolate. This is then not handled through the
736 // current function, but via interpolation matrices that
737 // result in constraints, rather than identities. Since
738 // that happens in a different function, there is nothing
739 // for us to do here.
740 }
741
742 return identities;
743 }
744 else if (const FE_Q<dim, spacedim> *fe_q_other =
745 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
746 {
747 // dofs are located along lines, so two dofs are identical if they are
748 // located at identical positions. if we had only equidistant points, we
749 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
750 // might have other support points (e.g. Gauss-Lobatto
751 // points). Therefore, read the points in unit_support_points for the
752 // first coordinate direction. For FE_Q, we take the lexicographic
753 // ordering of the line support points in the first direction (i.e.,
754 // x-direction), which we access between index 1 and p-1 (index 0 and p
755 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
756 // iterate over points on the first line which begin after the 3 vertex
757 // points in the complete list of unit support points
758
759 const std::vector<unsigned int> &index_map_inverse_q_other =
760 fe_q_other->get_poly_space_numbering_inverse();
761
762 std::vector<std::pair<unsigned int, unsigned int>> identities;
763
764 for (unsigned int i = 0; i < this->degree - 1; ++i)
765 for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
766 if (std::fabs(this->unit_support_points[i + 3][0] -
767 fe_q_other->get_unit_support_points()
768 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
769 identities.emplace_back(i, j);
770 else
771 {
772 // If nodes are not located in the same place, we have to
773 // interpolate. This will then also
774 // capture the case where the FE_Q has a different polynomial
775 // degree than the current element. In either case, the resulting
776 // constraints are computed elsewhere, rather than via the
777 // identities this function returns: Since
778 // that happens in a different function, there is nothing
779 // for us to do here.
780 }
781
782 return identities;
783 }
784 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
785 {
786 // The FE_Nothing has no degrees of freedom, so there are no
787 // equivalencies to be recorded. (If the FE_Nothing is dominating,
788 // then this will also leads to constraints, but we are not concerned
789 // with this here.)
790 return {};
791 }
792 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
793 {
794 // if the other element has no elements on faces at all,
795 // then it would be impossible to enforce any kind of
796 // continuity even if we knew exactly what kind of element
797 // we have -- simply because the other element declares
798 // that it is discontinuous because it has no DoFs on
799 // its faces. in that case, just state that we have no
800 // constraints to declare
801 return {};
802 }
803 else
804 {
805 Assert(false, ExcNotImplemented());
806 return {};
807 }
808}
809
810
811
812template <int dim, int spacedim>
814 : FE_SimplexPoly<dim, spacedim>(
815 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
816 FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
817 ReferenceCells::get_simplex<dim>(),
818 1,
819 degree,
820 FiniteElementData<dim>::L2),
821 unit_support_points_fe_p<dim>(degree),
822 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::L2),
823 constraints_fe_p<dim>(degree))
824{}
825
826
827
828template <int dim, int spacedim>
829std::unique_ptr<FiniteElement<dim, spacedim>>
831{
832 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
833}
834
835
836
837template <int dim, int spacedim>
838std::string
840{
841 std::ostringstream namebuf;
842 namebuf << "FE_SimplexDGP<" << dim << ">(" << this->degree << ")";
843
844 return namebuf.str();
845}
846
847
848template <int dim, int spacedim>
851 const FiniteElement<dim, spacedim> &fe_other,
852 const unsigned int codim) const
853{
854 Assert(codim <= dim, ExcImpossibleInDim(dim));
855
856 // vertex/line/face domination
857 // ---------------------------
858 if (codim > 0)
859 // this is a discontinuous element, so by definition there will
860 // be no constraints wherever this element comes together with
861 // any other kind of element
863
864 // cell domination
865 // ---------------
866 if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
867 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
868 {
869 if (this->degree < fe_dgp_other->degree)
871 else if (this->degree == fe_dgp_other->degree)
873 else
875 }
876 else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
877 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
878 {
879 if (this->degree < fe_dgq_other->degree)
881 else if (this->degree == fe_dgq_other->degree)
883 else
885 }
886 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
887 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
888 {
889 if (fe_nothing->is_dominating())
891 else
892 // the FE_Nothing has no degrees of freedom and it is typically used
893 // in a context where we don't require any continuity along the
894 // interface
896 }
897
898 Assert(false, ExcNotImplemented());
900}
901
902
903
904template <int dim, int spacedim>
905std::vector<std::pair<unsigned int, unsigned int>>
907 const FiniteElement<dim, spacedim> &fe_other) const
908{
909 (void)fe_other;
910
911 return {};
912}
913
914
915
916template <int dim, int spacedim>
917std::vector<std::pair<unsigned int, unsigned int>>
919 const FiniteElement<dim, spacedim> &fe_other) const
920{
921 (void)fe_other;
922
923 return {};
924}
925
926// explicit instantiations
927#include "fe_simplex_p.inst"
928
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition fe_q.h:551
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)
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
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
FE_SimplexPoly(const BarycentricPolynomials< dim > polynomials, const FiniteElementData< dim > &fe_data, 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 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
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:2447
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
std::vector< Point< dim > > unit_support_points
Definition fe.h:2440
FullMatrix< double > interface_constraints
Definition fe.h:2428
size_type m() const
Definition point.h:112
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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int vertex_indices[2]
const unsigned int v0
const unsigned int v1
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)
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)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
STL namespace.