Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_simplex_p.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
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
27
28namespace
29{
30 // TODO: replace this with QProjector once QProjector learns how to project
31 // quadrature points onto separate faces of simplices
32 template <int dim>
34 face_midpoint(const unsigned int face_no)
35 {
36 const auto reference_cell = ReferenceCells::get_simplex<dim>();
37 const auto face_reference_cell =
38 reference_cell.face_reference_cell(face_no);
39
40 Point<dim> midpoint;
41 for (const auto face_vertex_no : face_reference_cell.vertex_indices())
42 {
43 const auto vertex_no = reference_cell.face_to_cell_vertices(
44 face_no,
45 face_vertex_no,
47
48 midpoint += reference_cell.template vertex<dim>(vertex_no);
49 }
50
51 midpoint /= reference_cell.face_reference_cell(0).n_vertices();
52 return midpoint;
53 }
59 std::vector<unsigned int>
60 get_dpo_vector_fe_p(const unsigned int dim, const unsigned int degree)
61 {
62 switch (dim)
63 {
64 case 1:
65 switch (degree)
66 {
67 case 1:
68 return {1, 0};
69 case 2:
70 return {1, 1};
71 case 3:
72 return {1, 2};
73 default:
75 }
76 case 2:
77 switch (degree)
78 {
79 case 1:
80 return {1, 0, 0};
81 case 2:
82 return {1, 1, 0};
83 case 3:
84 return {1, 2, 1};
85 default:
87 }
88 case 3:
89 switch (degree)
90 {
91 case 1:
92 return {1, 0, 0, 0};
93 case 2:
94 return {1, 1, 0, 0};
95 case 3:
96 return {1, 2, 1, 0};
97 default:
99 }
101
103 return {};
104 }
105
107
112 template <int dim>
113 std::vector<Point<dim>>
114 unit_support_points_fe_p(const unsigned int degree)
115 {
116 Assert(dim != 0, ExcInternalError());
117 Assert(degree <= 3, ExcNotImplemented());
118 std::vector<Point<dim>> unit_points;
119 const auto reference_cell = ReferenceCells::get_simplex<dim>();
120
121 // Piecewise constants are a special case: use a support point at the
122 // centroid and only the centroid
123 if (degree == 0)
124 {
125 unit_points.emplace_back(reference_cell.template barycenter<dim>());
126 return unit_points;
127 }
128
129 // otherwise write everything as linear combinations of vertices
130 const auto dpo = get_dpo_vector_fe_p(dim, degree);
131 Assert(dpo.size() == dim + 1, ExcInternalError());
132 Assert(dpo[0] == 1, ExcNotImplemented());
133 // vertices:
134 for (const unsigned int d : reference_cell.vertex_indices())
135 unit_points.push_back(reference_cell.template vertex<dim>(d));
136 // lines:
137 for (const unsigned int l : reference_cell.line_indices())
138 {
139 const Point<dim> p0 =
140 unit_points[reference_cell.line_to_cell_vertices(l, 0)];
141 const Point<dim> p1 =
142 unit_points[reference_cell.line_to_cell_vertices(l, 1)];
143 for (unsigned int p = 0; p < dpo[1]; ++p)
144 unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 +
145 (double(p + 1) / (dpo[1] + 1)) * p1);
146 }
147 // quads:
148 if (dim > 1 && dpo[2] > 0)
149 {
150 Assert(dpo[2] == 1, ExcNotImplemented());
151 if (dim == 2)
152 unit_points.push_back(reference_cell.template barycenter<dim>());
153 if (dim == 3)
154 for (const unsigned int f : reference_cell.face_indices())
155 unit_points.push_back(face_midpoint<dim>(f));
156 }
157
158 return unit_points;
159 }
160
161 template <>
162 std::vector<Point<0>>
163 unit_support_points_fe_p(const unsigned int /*degree*/)
164 {
165 return {Point<0>()};
166 }
167
172 template <int dim>
173 std::vector<std::vector<Point<dim - 1>>>
174 unit_face_support_points_fe_p(
175 const unsigned int degree,
176 typename FiniteElementData<dim>::Conformity conformity)
177 {
178 // Discontinuous elements don't have face support points
180 return {};
181
182 // this concept doesn't exist in 1d so just return an empty vector
183 if (dim == 1)
184 return {};
185
186 std::vector<std::vector<Point<dim - 1>>> unit_face_points;
187
188 // all faces have the same support points
189 for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
190 {
191 (void)face_n;
192 unit_face_points.emplace_back(
193 unit_support_points_fe_p<dim - 1>(degree));
194 }
195
196 return unit_face_points;
197 }
198
204 template <int dim>
206 constraints_fe_p(const unsigned int /*degree*/)
207 {
208 // no constraints in 1d
209 // constraints in 3d not implemented yet
210 return FullMatrix<double>();
211 }
212
213 template <>
215 constraints_fe_p<2>(const unsigned int degree)
216 {
217 constexpr int dim = 2;
218 Assert(degree <= 3, ExcNotImplemented());
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 if (degree == 2)
230 {
231 // midpoint on subface 0
232 constraint_points.emplace_back(0.25);
233 // midpoint on subface 1
234 constraint_points.emplace_back(0.75);
235 }
236
237 // Now construct relation between destination (child) and source (mother)
238 // dofs.
239
240 const unsigned int n_dofs_constrained = constraint_points.size();
241 unsigned int n_dofs_per_face = degree + 1;
242 FullMatrix<double> interface_constraints(n_dofs_constrained,
243 n_dofs_per_face);
244
245 const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
246
247 for (unsigned int i = 0; i < n_dofs_constrained; ++i)
248 for (unsigned int j = 0; j < n_dofs_per_face; ++j)
249 {
250 interface_constraints(i, j) =
251 poly.compute_value(j, constraint_points[i]);
252
253 // if the value is small up to round-off, then simply set it to zero
254 // to avoid unwanted fill-in of the constraint matrices (which would
255 // then increase the number of other DoFs a constrained DoF would
256 // couple to)
257 if (std::fabs(interface_constraints(i, j)) < 1e-13)
258 interface_constraints(i, j) = 0;
259 }
260 return interface_constraints;
261 }
262
263
264
269 std::vector<unsigned int>
270 get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
271 {
272 // First treat the case of piecewise constant elements:
273 if (degree == 0)
274 {
275 std::vector<unsigned int> dpo(dim + 1, 0U);
276 dpo[dim] = 1;
277 return dpo;
278 }
279 else
280 {
281 // This element has the same degrees of freedom as the continuous one,
282 // but they are all counted for the interior of the cell because
283 // it is continuous. Rather than hard-code how many DoFs the element
284 // has, we just get the numbers from the continuous case and add them
285 // up
286 const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
287
288 switch (dim)
289 {
290 case 1:
291 return {0U,
292 ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
293 continuous_dpo[dim]};
294
295 case 2:
296 return {0U,
297 0U,
299 continuous_dpo[0] +
300 ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
301 continuous_dpo[dim]};
302
303 case 3:
304 return {
305 0U,
306 0U,
307 0U,
308 ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
309 ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
310 ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
311 continuous_dpo[dim]};
312 }
313
315 return {};
316 }
317 }
318} // namespace
319
320
321
322template <int dim, int spacedim>
324 const BarycentricPolynomials<dim> polynomials,
325 const FiniteElementData<dim> &fe_data,
326 const std::vector<Point<dim>> &unit_support_points,
327 const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
328 const FullMatrix<double> &interface_constraints)
329 : ::FE_Poly<dim, spacedim>(
330 polynomials,
331 fe_data,
332 std::vector<bool>(fe_data.dofs_per_cell),
333 std::vector<ComponentMask>(fe_data.dofs_per_cell,
334 ComponentMask(std::vector<bool>(1, true))))
335{
338 this->interface_constraints = interface_constraints;
339}
340
341
342
343template <int dim, int spacedim>
344std::pair<Table<2, bool>, std::vector<unsigned int>>
346{
347 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
348 constant_modes.fill(true);
349 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
350 constant_modes, std::vector<unsigned int>(1, 0));
351}
352
353
354
355template <int dim, int spacedim>
356const FullMatrix<double> &
358 const unsigned int child,
359 const RefinementCase<dim> &refinement_case) const
360{
363 AssertDimension(dim, spacedim);
364
365 // initialization upon first request
366 if (this->prolongation[refinement_case - 1][child].n() == 0)
367 {
368 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
369
370 // if matrix got updated while waiting for the lock
371 if (this->prolongation[refinement_case - 1][child].n() ==
372 this->n_dofs_per_cell())
373 return this->prolongation[refinement_case - 1][child];
374
375 // now do the work. need to get a non-const version of data in order to
376 // be able to modify them inside a const function
377 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
378
379 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
381 isotropic_matrices.back().resize(
383 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
384
385 FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
386
387 this_nonconst.prolongation[refinement_case - 1] =
388 std::move(isotropic_matrices.back());
389 }
390
391 // finally return the matrix
392 return this->prolongation[refinement_case - 1][child];
393}
394
395
396
397template <int dim, int spacedim>
398const FullMatrix<double> &
400 const unsigned int child,
401 const RefinementCase<dim> &refinement_case) const
402{
405 AssertDimension(dim, spacedim);
406
407 // initialization upon first request
408 if (this->restriction[refinement_case - 1][child].n() == 0)
409 {
410 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
411
412 // if matrix got updated while waiting for the lock
413 if (this->restriction[refinement_case - 1][child].n() ==
414 this->n_dofs_per_cell())
415 return this->restriction[refinement_case - 1][child];
416
417 // now do the work. need to get a non-const version of data in order to
418 // be able to modify them inside a const function
419 auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
420
421 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
423 isotropic_matrices.back().resize(
425 FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
426
427 FETools::compute_projection_matrices(*this, isotropic_matrices, true);
428
429 this_nonconst.restriction[refinement_case - 1] =
430 std::move(isotropic_matrices.back());
431 }
432
433 // finally return the matrix
434 return this->restriction[refinement_case - 1][child];
435}
436
437
438
439template <int dim, int spacedim>
440void
442 const FiniteElement<dim, spacedim> &source_fe,
443 FullMatrix<double> &interpolation_matrix,
444 const unsigned int face_no) const
445{
446 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
447 ExcDimensionMismatch(interpolation_matrix.m(),
448 source_fe.n_dofs_per_face(face_no)));
449
450 // see if source is a P or Q element
451 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
452 nullptr) ||
453 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
454 {
455 const Quadrature<dim - 1> quad_face_support(
456 source_fe.get_unit_face_support_points(face_no));
457
458 const double eps = 2e-13 * this->degree * (dim - 1);
459
460 std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
461 QProjector<dim>::project_to_face(this->reference_cell(),
462 quad_face_support,
463 face_no,
464 face_quadrature_points);
465
466 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
467 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
468 {
469 double matrix_entry =
470 this->shape_value(this->face_to_cell_index(j, 0),
471 face_quadrature_points[i]);
472
473 // Correct the interpolated value. I.e. if it is close to 1 or
474 // 0, make it exactly 1 or 0. Unfortunately, this is required to
475 // avoid problems with higher order elements.
476 if (std::fabs(matrix_entry - 1.0) < eps)
477 matrix_entry = 1.0;
478 if (std::fabs(matrix_entry) < eps)
479 matrix_entry = 0.0;
480
481 interpolation_matrix(i, j) = matrix_entry;
482 }
483
484#ifdef DEBUG
485 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
486 {
487 double sum = 0.;
488
489 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
490 sum += interpolation_matrix(j, i);
491
492 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
493 }
494#endif
495 }
496 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
497 {
498 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
499 }
500 else
502 false,
503 (typename FiniteElement<dim,
504 spacedim>::ExcInterpolationNotImplemented()));
505}
506
507
508
509template <int dim, int spacedim>
510void
512 const FiniteElement<dim, spacedim> &source_fe,
513 const unsigned int subface,
514 FullMatrix<double> &interpolation_matrix,
515 const unsigned int face_no) const
516{
517 Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
518 ExcDimensionMismatch(interpolation_matrix.m(),
519 source_fe.n_dofs_per_face(face_no)));
520
521 // see if source is a P or Q element
522 if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
523 nullptr) ||
524 (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
525 {
526 const Quadrature<dim - 1> quad_face_support(
527 source_fe.get_unit_face_support_points(face_no));
528
529 const double eps = 2e-13 * this->degree * (dim - 1);
530
531 std::vector<Point<dim>> subface_quadrature_points(
532 quad_face_support.size());
533 QProjector<dim>::project_to_subface(this->reference_cell(),
534 quad_face_support,
535 face_no,
536 subface,
537 subface_quadrature_points);
538
539 for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
540 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
541 {
542 double matrix_entry =
543 this->shape_value(this->face_to_cell_index(j, 0),
544 subface_quadrature_points[i]);
545
546 // Correct the interpolated value. I.e. if it is close to 1 or
547 // 0, make it exactly 1 or 0. Unfortunately, this is required to
548 // avoid problems with higher order elements.
549 if (std::fabs(matrix_entry - 1.0) < eps)
550 matrix_entry = 1.0;
551 if (std::fabs(matrix_entry) < eps)
552 matrix_entry = 0.0;
553
554 interpolation_matrix(i, j) = matrix_entry;
555 }
556
557#ifdef DEBUG
558 for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
559 {
560 double sum = 0.;
561
562 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
563 sum += interpolation_matrix(j, i);
564
565 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
566 }
567#endif
568 }
569 else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
570 {
571 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
572 }
573 else
575 false,
576 (typename FiniteElement<dim,
577 spacedim>::ExcInterpolationNotImplemented()));
578}
579
580
581
582template <int dim, int spacedim>
583bool
588
589
590
591template <int dim, int spacedim>
592void
595 const std::vector<Vector<double>> &support_point_values,
596 std::vector<double> &nodal_values) const
597{
598 AssertDimension(support_point_values.size(),
599 this->get_unit_support_points().size());
600 AssertDimension(support_point_values.size(), nodal_values.size());
601 AssertDimension(this->dofs_per_cell, nodal_values.size());
602
603 for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
604 {
605 AssertDimension(support_point_values[i].size(), 1);
606
607 nodal_values[i] = support_point_values[i](0);
608 }
609}
610
611
612
613template <int dim, int spacedim>
615 : FE_SimplexPoly<dim, spacedim>(
616 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
617 FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
618 ReferenceCells::get_simplex<dim>(),
619 1,
620 degree,
621 FiniteElementData<dim>::H1),
622 unit_support_points_fe_p<dim>(degree),
623 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
624 constraints_fe_p<dim>(degree))
625{
626 if (degree > 2)
627 for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i)
629 this->n_dofs_per_line() - 1 - i - i;
630 // We do not support multiple DoFs per quad yet
632}
633
634
635
636template <int dim, int spacedim>
637std::unique_ptr<FiniteElement<dim, spacedim>>
639{
640 return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
641}
642
643
644
645template <int dim, int spacedim>
646std::string
648{
649 std::ostringstream namebuf;
650 namebuf << "FE_SimplexP<" << Utilities::dim_string(dim, spacedim) << ">("
651 << this->degree << ")";
652
653 return namebuf.str();
654}
655
656
657
658template <int dim, int spacedim>
661 const FiniteElement<dim, spacedim> &fe_other,
662 const unsigned int codim) const
663{
664 Assert(codim <= dim, ExcImpossibleInDim(dim));
665
666 // vertex/line/face domination
667 // (if fe_other is derived from FE_SimplexDGP)
668 // ------------------------------------
669 if (codim > 0)
670 if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
671 nullptr)
672 // there are no requirements between continuous and discontinuous
673 // elements
675
676 // vertex/line/face domination
677 // (if fe_other is not derived from FE_SimplexDGP)
678 // & cell domination
679 // ----------------------------------------
680 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
681 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
682 {
683 if (this->degree < fe_p_other->degree)
685 else if (this->degree == fe_p_other->degree)
687 else
689 }
690 else if (const FE_Q<dim, spacedim> *fe_q_other =
691 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
692 {
693 if (this->degree < fe_q_other->degree)
695 else if (this->degree == fe_q_other->degree)
697 else
699 }
700 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
701 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
702 {
703 if (fe_nothing->is_dominating())
705 else
706 // the FE_Nothing has no degrees of freedom and it is typically used
707 // in a context where we don't require any continuity along the
708 // interface
710 }
711
714}
715
716
717
718template <int dim, int spacedim>
719std::vector<std::pair<unsigned int, unsigned int>>
721 const FiniteElement<dim, spacedim> &fe_other) const
722{
723 AssertDimension(dim, 2);
724
725 if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
726 {
727 // there should be exactly one single DoF of each FE at a vertex, and
728 // they should have identical value
729 return {{0U, 0U}};
730 }
731 else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
732 {
733 // there should be exactly one single DoF of each FE at a vertex, and
734 // they should have identical value
735 return {{0U, 0U}};
736 }
737 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
738 {
739 // the FE_Nothing has no degrees of freedom, so there are no
740 // equivalencies to be recorded
741 return {};
742 }
743 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
744 {
745 // if the other element has no DoFs on faces at all,
746 // then it would be impossible to enforce any kind of
747 // continuity even if we knew exactly what kind of element
748 // we have -- simply because the other element declares
749 // that it is discontinuous because it has no DoFs on
750 // its faces. in that case, just state that we have no
751 // constraints to declare
752 return {};
753 }
754 else
755 {
757 return {};
758 }
759}
760
761
762
763template <int dim, int spacedim>
764std::vector<std::pair<unsigned int, unsigned int>>
766 const FiniteElement<dim, spacedim> &fe_other) const
767{
768 AssertDimension(dim, 2);
769
770 if (const FE_SimplexP<dim, spacedim> *fe_p_other =
771 dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
772 {
773 // dofs are located along lines, so two dofs are identical if they are
774 // located at identical positions.
775 // Therefore, read the points in unit_support_points for the
776 // first coordinate direction. For FE_SimplexP, they are currently
777 // hard-coded and we iterate over points on the first line which begin
778 // after the 3 vertex points in the complete list of unit support points
779
780 std::vector<std::pair<unsigned int, unsigned int>> identities;
781
782 for (unsigned int i = 0; i < this->degree - 1; ++i)
783 for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
784 if (std::fabs(this->unit_support_points[i + 3][0] -
785 fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
786 identities.emplace_back(i, j);
787 else
788 {
789 // If nodes are not located in the same place, we have to
790 // interpolate. This is then not handled through the
791 // current function, but via interpolation matrices that
792 // result in constraints, rather than identities. Since
793 // that happens in a different function, there is nothing
794 // for us to do here.
795 }
796
797 return identities;
798 }
799 else if (const FE_Q<dim, spacedim> *fe_q_other =
800 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
801 {
802 // dofs are located along lines, so two dofs are identical if they are
803 // located at identical positions. if we had only equidistant points, we
804 // could simply check for similarity like (i+1)*q == (j+1)*p, but we
805 // might have other support points (e.g. Gauss-Lobatto
806 // points). Therefore, read the points in unit_support_points for the
807 // first coordinate direction. For FE_Q, we take the lexicographic
808 // ordering of the line support points in the first direction (i.e.,
809 // x-direction), which we access between index 1 and p-1 (index 0 and p
810 // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
811 // iterate over points on the first line which begin after the 3 vertex
812 // points in the complete list of unit support points
813
814 const std::vector<unsigned int> &index_map_inverse_q_other =
815 fe_q_other->get_poly_space_numbering_inverse();
816
817 std::vector<std::pair<unsigned int, unsigned int>> identities;
818
819 for (unsigned int i = 0; i < this->degree - 1; ++i)
820 for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
821 if (std::fabs(this->unit_support_points[i + 3][0] -
822 fe_q_other->get_unit_support_points()
823 [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
824 identities.emplace_back(i, j);
825 else
826 {
827 // If nodes are not located in the same place, we have to
828 // interpolate. This will then also
829 // capture the case where the FE_Q has a different polynomial
830 // degree than the current element. In either case, the resulting
831 // constraints are computed elsewhere, rather than via the
832 // identities this function returns: Since
833 // that happens in a different function, there is nothing
834 // for us to do here.
835 }
836
837 return identities;
838 }
839 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
840 {
841 // The FE_Nothing has no degrees of freedom, so there are no
842 // equivalencies to be recorded. (If the FE_Nothing is dominating,
843 // then this will also leads to constraints, but we are not concerned
844 // with this here.)
845 return {};
846 }
847 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
848 {
849 // if the other element has no elements on faces at all,
850 // then it would be impossible to enforce any kind of
851 // continuity even if we knew exactly what kind of element
852 // we have -- simply because the other element declares
853 // that it is discontinuous because it has no DoFs on
854 // its faces. in that case, just state that we have no
855 // constraints to declare
856 return {};
857 }
858 else
859 {
861 return {};
862 }
863}
864
865
866
867template <int dim, int spacedim>
869 : FE_SimplexPoly<dim, spacedim>(
870 BarycentricPolynomials<dim>::get_fe_p_basis(degree),
871 FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
872 ReferenceCells::get_simplex<dim>(),
873 1,
874 degree,
875 FiniteElementData<dim>::L2),
876 unit_support_points_fe_p<dim>(degree),
877 unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::L2),
878 constraints_fe_p<dim>(degree))
879{}
880
881
882
883template <int dim, int spacedim>
884std::unique_ptr<FiniteElement<dim, spacedim>>
886{
887 return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
888}
889
890
891
892template <int dim, int spacedim>
893std::string
895{
896 std::ostringstream namebuf;
897 namebuf << "FE_SimplexDGP<" << Utilities::dim_string(dim, spacedim) << ">("
898 << this->degree << ")";
899
900 return namebuf.str();
901}
902
903
904template <int dim, int spacedim>
907 const FiniteElement<dim, spacedim> &fe_other,
908 const unsigned int codim) const
909{
910 Assert(codim <= dim, ExcImpossibleInDim(dim));
911
912 // vertex/line/face domination
913 // ---------------------------
914 if (codim > 0)
915 // this is a discontinuous element, so by definition there will
916 // be no constraints wherever this element comes together with
917 // any other kind of element
919
920 // cell domination
921 // ---------------
922 if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
923 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
924 {
925 if (this->degree < fe_dgp_other->degree)
927 else if (this->degree == fe_dgp_other->degree)
929 else
931 }
932 else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
933 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
934 {
935 if (this->degree < fe_dgq_other->degree)
937 else if (this->degree == fe_dgq_other->degree)
939 else
941 }
942 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
943 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
944 {
945 if (fe_nothing->is_dominating())
947 else
948 // the FE_Nothing has no degrees of freedom and it is typically used
949 // in a context where we don't require any continuity along the
950 // interface
952 }
953
956}
957
958
959
960template <int dim, int spacedim>
961std::vector<std::pair<unsigned int, unsigned int>>
963 const FiniteElement<dim, spacedim> &fe_other) const
964{
965 (void)fe_other;
966
967 return {};
968}
969
970
971
972template <int dim, int spacedim>
973std::vector<std::pair<unsigned int, unsigned int>>
975 const FiniteElement<dim, spacedim> &fe_other) const
976{
977 (void)fe_other;
978
979 return {};
980}
981
982// explicit instantiations
983#include "fe_simplex_p.inst"
984
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition fe_q.h:550
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
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:2443
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:2485
std::vector< Point< dim > > unit_support_points
Definition fe.h:2436
FullMatrix< double > interface_constraints
Definition fe.h:2424
size_type m() const
Definition point.h:111
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
static constexpr unsigned char default_combined_face_orientation()
unsigned int n_faces() const
unsigned int n_lines() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
#define DEAL_II_NOT_IMPLEMENTED()
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 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
STL namespace.