Reference documentation for deal.II version 9.3.3
\(\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\}}\)
fe_simplex_p.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2020 - 2021 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 std::vector<unsigned int> dpo(dim + 1, 0U);
39
40 if (degree == 1)
41 {
42 // one dof at each vertex
43 dpo[0] = 1;
44 }
45 else if (degree == 2)
46 {
47 // one dof at each vertex and in the middle of each line
48 dpo[0] = 1;
49 dpo[1] = 1;
50 dpo[2] = 0;
51 }
52 else
53 {
54 Assert(false, ExcNotImplemented());
55 }
56
57 return dpo;
58 }
59
64 template <int dim>
65 std::vector<Point<dim>>
66 unit_support_points_fe_poly(const unsigned int degree)
67 {
68 std::vector<Point<dim>> unit_points;
70 // Piecewise constants are a special case: use a support point at the
71 // centroid and only the centroid
72 if (degree == 0)
73 {
74 Point<dim> centroid;
75 std::fill(centroid.begin_raw(),
76 centroid.end_raw(),
77 1.0 / double(dim + 1));
78 unit_points.emplace_back(centroid);
79 return unit_points;
80 }
81
82 if (dim == 1)
83 {
84 // We don't really have dim = 1 support for simplex elements yet, but
85 // its convenient for populating the face array
86 Assert(degree <= 2, ExcNotImplemented());
87 if (degree >= 1)
88 {
89 unit_points.emplace_back(0.0);
90 unit_points.emplace_back(1.0);
91
92 if (degree == 2)
93 unit_points.emplace_back(0.5);
94 }
95 }
96 else if (dim == 2)
97 {
98 Assert(degree <= 2, ExcNotImplemented());
99 if (degree >= 1)
100 {
101 unit_points.emplace_back(0.0, 0.0);
102 unit_points.emplace_back(1.0, 0.0);
103 unit_points.emplace_back(0.0, 1.0);
104
105 if (degree == 2)
106 {
107 unit_points.emplace_back(0.5, 0.0);
108 unit_points.emplace_back(0.5, 0.5);
109 unit_points.emplace_back(0.0, 0.5);
110 }
111 }
112 }
113 else if (dim == 3)
114 {
115 Assert(degree <= 2, ExcNotImplemented());
116 if (degree >= 1)
117 {
118 unit_points.emplace_back(0.0, 0.0, 0.0);
119 unit_points.emplace_back(1.0, 0.0, 0.0);
120 unit_points.emplace_back(0.0, 1.0, 0.0);
121 unit_points.emplace_back(0.0, 0.0, 1.0);
122
123 if (degree == 2)
124 {
125 unit_points.emplace_back(0.5, 0.0, 0.0);
126 unit_points.emplace_back(0.5, 0.5, 0.0);
127 unit_points.emplace_back(0.0, 0.5, 0.0);
128 unit_points.emplace_back(0.0, 0.0, 0.5);
129 unit_points.emplace_back(0.5, 0.0, 0.5);
130 unit_points.emplace_back(0.0, 0.5, 0.5);
131 }
132 }
133 }
134 else
135 {
136 Assert(false, ExcNotImplemented());
137 }
138
139 return unit_points;
140 }
141
146 template <int dim>
147 std::vector<std::vector<Point<dim - 1>>>
148 unit_face_support_points_fe_poly(const unsigned int degree)
149 {
150 // this concept doesn't exist in 1D so just return an empty vector
151 if (dim == 1)
152 return {};
153
154 std::vector<std::vector<Point<dim - 1>>> unit_face_points;
155
156 // all faces have the same support points
157 for (auto face_n :
159 .face_indices())
160 {
161 (void)face_n;
162 unit_face_points.emplace_back(
163 unit_support_points_fe_poly<dim - 1>(degree));
164 }
165
166 return unit_face_points;
167 }
168
174 template <int dim>
176 constraints_fe_poly(const unsigned int /*degree*/)
177 {
178 // no constraints in 1d
179 // constraints in 3d not implemented yet
180 return FullMatrix<double>();
181 }
182
183 template <>
185 constraints_fe_poly<2>(const unsigned int degree)
186 {
187 const unsigned int dim = 2;
188
189 Assert(degree <= 2, ExcNotImplemented());
190
191 // the following implements the 2d case
192 // (the 3d case is not implemented yet)
193 //
194 // consult FE_Q_Base::Implementation::initialize_constraints()
195 // for more information
196
197 std::vector<Point<dim - 1>> constraint_points;
198 // midpoint
199 constraint_points.emplace_back(0.5);
200 if (degree == 2)
201 {
202 // midpoint on subface 0
203 constraint_points.emplace_back(0.25);
204 // midpoint on subface 1
205 constraint_points.emplace_back(0.75);
206 }
207
208 // Now construct relation between destination (child) and source (mother)
209 // dofs.
210
211 const unsigned int n_dofs_constrained = constraint_points.size();
212 unsigned int n_dofs_per_face = degree + 1;
213 FullMatrix<double> interface_constraints(n_dofs_constrained,
214 n_dofs_per_face);
215
216 const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
217
218 for (unsigned int i = 0; i < n_dofs_constrained; ++i)
219 for (unsigned int j = 0; j < n_dofs_per_face; ++j)
220 {
221 interface_constraints(i, j) =
222 poly.compute_value(j, constraint_points[i]);
223
224 // if the value is small up to round-off, then simply set it to zero
225 // to avoid unwanted fill-in of the constraint matrices (which would
226 // then increase the number of other DoFs a constrained DoF would
227 // couple to)
228 if (std::fabs(interface_constraints(i, j)) < 1e-13)
229 interface_constraints(i, j) = 0;
230 }
231 return interface_constraints;
232 }
233
238 std::vector<unsigned int>
239 get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
240 {
241 std::vector<unsigned int> dpo(dim + 1, 0U);
242
243 // all dofs are internal
244 if (dim == 2 && degree == 1)
245 dpo[dim] = 3;
246 else if (dim == 2 && degree == 2)
247 dpo[dim] = 6;
248 else if (dim == 3 && degree == 1)
249 dpo[dim] = 4;
250 else if (dim == 3 && degree == 2)
251 dpo[dim] = 10;
252 else
253 {
254 Assert(false, ExcNotImplemented());
255 }
256
257 return dpo;
258 }
259} // namespace
260
261
262
263template <int dim, int spacedim>
265 const unsigned int degree,
266 const std::vector<unsigned int> & dpo_vector,
267 const typename FiniteElementData<dim>::Conformity conformity)
268 : ::FE_Poly<dim, spacedim>(
269 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
270 FiniteElementData<dim>(dpo_vector,
271 dim == 2 ? ReferenceCells::Triangle :
273 1,
274 degree,
275 conformity),
276 std::vector<bool>(FiniteElementData<dim>(dpo_vector,
277 dim == 2 ?
280 1,
281 degree)
282 .dofs_per_cell,
283 true),
284 std::vector<ComponentMask>(
285 FiniteElementData<dim>(dpo_vector,
286 dim == 2 ? ReferenceCells::Triangle :
288 1,
289 degree)
290 .dofs_per_cell,
291 std::vector<bool>(1, true)))
292{
293 this->unit_support_points = unit_support_points_fe_poly<dim>(degree);
294 // Discontinuous elements don't have face support points
297 unit_face_support_points_fe_poly<dim>(degree);
298 this->interface_constraints = constraints_fe_poly<dim>(degree);
299}
300
301
302
303template <int dim, int spacedim>
304std::pair<Table<2, bool>, std::vector<unsigned int>>
306{
307 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
308 constant_modes.fill(true);
309 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
310 constant_modes, std::vector<unsigned int>(1, 0));
311}
312
313
314
315template <int dim, int spacedim>
316const FullMatrix<double> &
318 const unsigned int child,
319 const RefinementCase<dim> &refinement_case) const
320{
323 AssertDimension(dim, spacedim);
324
325 // initialization upon first request
326 if (this->prolongation[refinement_case - 1][child].n() == 0)
327 {
328 std::lock_guard<std::mutex> lock(this->mutex);
329
330 // if matrix got updated while waiting for the lock
331 if (this->prolongation[refinement_case - 1][child].n() ==
332 this->n_dofs_per_cell())
333 return this->prolongation[refinement_case - 1][child];
334
335 // now do the work. need to get a non-const version of data in order to
336 // be able to modify them inside a const function
337 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
338
339 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
341 isotropic_matrices.back().resize(
343 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
344
345 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
346
347 this_nonconst.prolongation[refinement_case - 1].swap(
348 isotropic_matrices.back());
349 }
350
351 // finally return the matrix
352 return this->prolongation[refinement_case - 1][child];
353}
354
355
356
357template <int dim, int spacedim>
358const FullMatrix<double> &
360 const unsigned int child,
361 const RefinementCase<dim> &refinement_case) const
362{
365 AssertDimension(dim, spacedim);
366
367 // initialization upon first request
368 if (this->restriction[refinement_case - 1][child].n() == 0)
369 {
370 std::lock_guard<std::mutex> lock(this->mutex);
371
372 // if matrix got updated while waiting for the lock
373 if (this->restriction[refinement_case - 1][child].n() ==
374 this->n_dofs_per_cell())
375 return this->restriction[refinement_case - 1][child];
376
377 // now do the work. need to get a non-const version of data in order to
378 // be able to modify them inside a const function
379 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
380
381 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
383 isotropic_matrices.back().resize(
385 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
386
387 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
388
389 this_nonconst.restriction[refinement_case - 1].swap(
390 isotropic_matrices.back());
391 }
392
393 // finally return the matrix
394 return this->restriction[refinement_case - 1][child];
395}
396
397
398
399template <int dim, int spacedim>
400void
402 const FiniteElement<dim, spacedim> &source_fe,
403 FullMatrix<double> & interpolation_matrix,
404 const unsigned int face_no) const
405{
406 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
407 ExcDimensionMismatch(interpolation_matrix.m(),
408 source_fe.n_dofs_per_face(face_no)));
409
410 // see if source is a P or Q element
411 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
412 nullptr) ||
413 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
414 {
415 const Quadrature<dim - 1> quad_face_support(
416 source_fe.get_unit_face_support_points(face_no));
417
418 const double eps = 2e-13 * this->degree * (dim - 1);
419
420 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
422 quad_face_support,
423 face_no,
424 face_quadrature_points);
425
426 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
427 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
428 {
429 double matrix_entry =
430 this->shape_value(this->face_to_cell_index(j, 0),
431 face_quadrature_points[i]);
432
433 // Correct the interpolated value. I.e. if it is close to 1 or
434 // 0, make it exactly 1 or 0. Unfortunately, this is required to
435 // avoid problems with higher order elements.
436 if (std::fabs(matrix_entry - 1.0) < eps)
437 matrix_entry = 1.0;
438 if (std::fabs(matrix_entry) < eps)
439 matrix_entry = 0.0;
440
441 interpolation_matrix(i, j) = matrix_entry;
442 }
443
444#ifdef DEBUG
445 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
446 {
447 double sum = 0.;
448
449 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
450 sum += interpolation_matrix(j, i);
451
453 }
454#endif
455 }
456 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
457 {
458 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
459 }
460 else
462 false,
463 (typename FiniteElement<dim,
464 spacedim>::ExcInterpolationNotImplemented()));
465}
466
467
468
469template <int dim, int spacedim>
470void
472 const FiniteElement<dim, spacedim> &source_fe,
473 const unsigned int subface,
474 FullMatrix<double> & interpolation_matrix,
475 const unsigned int face_no) const
476{
477 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
478 ExcDimensionMismatch(interpolation_matrix.m(),
479 source_fe.n_dofs_per_face(face_no)));
480
481 // see if source is a P or Q element
482 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
483 nullptr) ||
484 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
485 {
486 const Quadrature<dim - 1> quad_face_support(
487 source_fe.get_unit_face_support_points(face_no));
488
489 const double eps = 2e-13 * this->degree * (dim - 1);
490
491 std::vector<Point<dim>> subface_quadrature_points(
492 quad_face_support.size());
494 quad_face_support,
495 face_no,
496 subface,
497 subface_quadrature_points);
498
499 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
500 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
501 {
502 double matrix_entry =
503 this->shape_value(this->face_to_cell_index(j, 0),
504 subface_quadrature_points[i]);
505
506 // Correct the interpolated value. I.e. if it is close to 1 or
507 // 0, make it exactly 1 or 0. Unfortunately, this is required to
508 // avoid problems with higher order elements.
509 if (std::fabs(matrix_entry - 1.0) < eps)
510 matrix_entry = 1.0;
511 if (std::fabs(matrix_entry) < eps)
512 matrix_entry = 0.0;
513
514 interpolation_matrix(i, j) = matrix_entry;
515 }
516
517#ifdef DEBUG
518 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
519 {
520 double sum = 0.;
521
522 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
523 sum += interpolation_matrix(j, i);
524
526 }
527#endif
528 }
529 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
530 {
531 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
532 }
533 else
535 false,
536 (typename FiniteElement<dim,
537 spacedim>::ExcInterpolationNotImplemented()));
538}
539
540
541
542template <int dim, int spacedim>
543bool
545{
546 return true;
547}
548
549
550
551template <int dim, int spacedim>
552void
555 const std::vector<Vector<double>> &support_point_values,
556 std::vector<double> & nodal_values) const
557{
558 AssertDimension(support_point_values.size(),
559 this->get_unit_support_points().size());
560 AssertDimension(support_point_values.size(), nodal_values.size());
561 AssertDimension(this->dofs_per_cell, nodal_values.size());
562
563 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
564 {
565 AssertDimension(support_point_values[i].size(), 1);
566
567 nodal_values[i] = support_point_values[i](0);
568 }
569}
570
571
572
573template <int dim, int spacedim>
575 : FE_SimplexPoly<dim, spacedim>(degree,
576 get_dpo_vector_fe_p(dim, degree),
577 FiniteElementData<dim>::H1)
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 elements 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 Assert(this->degree <= 2, ExcNotImplemented());
715
716 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
717 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
718 {
719 // dofs are located along lines, so two dofs are identical if they are
720 // located at identical positions.
721 // Therefore, read the points in unit_support_points for the
722 // first coordinate direction. For FE_SimplexP, they are currently
723 // hard-coded and we iterate over points on the first line which begin
724 // after the 3 vertex points in the complete list of unit support points
725
726 Assert(fe_p_other->degree <= 2, ExcNotImplemented());
727
728 std::vector<std::pair<unsigned int, unsigned int>> identities;
729
730 for (unsigned int i = 0; i < this->degree - 1; ++i)
731 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
732 if (std::fabs(this->unit_support_points[i + 3][0] -
733 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
734 identities.emplace_back(i, j);
735
736 return identities;
737 }
738 else if (const FE_Q<dim, spacedim> *fe_q_other =
739 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
740 {
741 // dofs are located along lines, so two dofs are identical if they are
742 // located at identical positions. if we had only equidistant points, we
743 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
744 // might have other support points (e.g. Gauss-Lobatto
745 // points). Therefore, read the points in unit_support_points for the
746 // first coordinate direction. For FE_Q, we take the lexicographic
747 // ordering of the line support points in the first direction (i.e.,
748 // x-direction), which we access between index 1 and p-1 (index 0 and p
749 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
750 // iterate over points on the first line which begin after the 3 vertex
751 // points in the complete list of unit support points
752
753 const std::vector<unsigned int> &index_map_inverse_q_other =
754 fe_q_other->get_poly_space_numbering_inverse();
755
756 std::vector<std::pair<unsigned int, unsigned int>> identities;
757
758 for (unsigned int i = 0; i < this->degree - 1; ++i)
759 for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
760 if (std::fabs(this->unit_support_points[i + 3][0] -
761 fe_q_other->get_unit_support_points()
762 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
763 identities.emplace_back(i, j);
764
765 return identities;
766 }
767 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
768 {
769 // the FE_Nothing has no degrees of freedom, so there are no
770 // equivalencies to be recorded
771 return {};
772 }
773 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
774 {
775 // if the other element has no elements on faces at all,
776 // then it would be impossible to enforce any kind of
777 // continuity even if we knew exactly what kind of element
778 // we have -- simply because the other element declares
779 // that it is discontinuous because it has no DoFs on
780 // its faces. in that case, just state that we have no
781 // constraints to declare
782 return {};
783 }
784 else
785 {
786 Assert(false, ExcNotImplemented());
787 return {};
788 }
789}
790
791
792
793template <int dim, int spacedim>
795 : FE_SimplexPoly<dim, spacedim>(degree,
796 get_dpo_vector_fe_dgp(dim, degree),
797 FiniteElementData<dim>::L2)
798{}
799
800
801
802template <int dim, int spacedim>
803std::unique_ptr<FiniteElement<dim, spacedim>>
805{
806 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
807}
808
809
810
811template <int dim, int spacedim>
812std::string
814{
815 std::ostringstream namebuf;
816 namebuf << "FE_SimplexDGP<" << dim << ">(" << this->degree << ")";
817
818 return namebuf.str();
819}
820
821
822template <int dim, int spacedim>
825 const FiniteElement<dim, spacedim> &fe_other,
826 const unsigned int codim) const
827{
828 Assert(codim <= dim, ExcImpossibleInDim(dim));
829
830 // vertex/line/face domination
831 // ---------------------------
832 if (codim > 0)
833 // this is a discontinuous element, so by definition there will
834 // be no constraints wherever this element comes together with
835 // any other kind of element
837
838 // cell domination
839 // ---------------
840 if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
841 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
842 {
843 if (this->degree < fe_dgp_other->degree)
845 else if (this->degree == fe_dgp_other->degree)
847 else
849 }
850 else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
851 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
852 {
853 if (this->degree < fe_dgq_other->degree)
855 else if (this->degree == fe_dgq_other->degree)
857 else
859 }
860 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
861 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
862 {
863 if (fe_nothing->is_dominating())
865 else
866 // the FE_Nothing has no degrees of freedom and it is typically used
867 // in a context where we don't require any continuity along the
868 // interface
870 }
871
872 Assert(false, ExcNotImplemented());
874}
875
876
877
878template <int dim, int spacedim>
879std::vector<std::pair<unsigned int, unsigned int>>
881 const FiniteElement<dim, spacedim> &fe_other) const
882{
883 (void)fe_other;
884
885 return {};
886}
887
888
889
890template <int dim, int spacedim>
891std::vector<std::pair<unsigned int, unsigned int>>
893 const FiniteElement<dim, spacedim> &fe_other) const
894{
895 (void)fe_other;
896
897 return {};
898}
899
900// explicit instantiations
901#include "fe_simplex_p.inst"
902
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition: fe_dgq.h:111
Definition: fe_q.h:549
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
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
FE_SimplexPoly(const unsigned int degree, const std::vector< unsigned int > &dpo_vector, const typename FiniteElementData< dim >::Conformity conformity)
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_base.h:435
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:2449
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:2442
FullMatrix< double > interface_constraints
Definition: fe.h:2430
size_type m() const
Definition: point.h:111
static void project_to_subface(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 SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim > > &q_points)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Number * begin_raw()
Number * end_raw()
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
Expression fabs(const Expression &x)
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)
static const char U
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 > > &line_support_points, const std::vector< unsigned int > &renumbering)
STL namespace.