Reference documentation for deal.II version 9.6.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_face.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 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
17
18#include <deal.II/fe/fe_face.h>
20#include <deal.II/fe/fe_poly_face.templates.h>
21#include <deal.II/fe/fe_tools.h>
22
24
25#include <memory>
26#include <sstream>
27
29
30
31namespace internal
32{
33 namespace FE_FaceQImplementation
34 {
35 namespace
36 {
37 std::vector<Point<1>>
38 get_QGaussLobatto_points(const unsigned int degree)
39 {
40 if (degree > 0)
41 return QGaussLobatto<1>(degree + 1).get_points();
42 else
43 return std::vector<Point<1>>(1, Point<1>(0.5));
44 }
45 } // namespace
46 } // namespace FE_FaceQImplementation
47} // namespace internal
48
49template <int dim, int spacedim>
50FE_FaceQ<dim, spacedim>::FE_FaceQ(const unsigned int degree)
51 : FE_PolyFace<TensorProductPolynomials<dim - 1>, dim, spacedim>(
53 Polynomials::generate_complete_Lagrange_basis(
54 internal::FE_FaceQImplementation::get_QGaussLobatto_points(degree))),
55 FiniteElementData<dim>(get_dpo_vector(degree),
56 1,
57 degree,
58 FiniteElementData<dim>::L2),
59 std::vector<bool>(1, true))
60{
61 // initialize unit face support points
62 const unsigned int codim = dim - 1;
63 this->unit_face_support_points[0].resize(
64 Utilities::fixed_power<codim>(this->degree + 1));
65
66 if (this->degree == 0)
67 for (unsigned int d = 0; d < codim; ++d)
68 this->unit_face_support_points[0][0][d] = 0.5;
69 else
70 {
71 std::vector<Point<1>> points =
72 internal::FE_FaceQImplementation::get_QGaussLobatto_points(degree);
73
74 unsigned int k = 0;
75 for (unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
76 for (unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
77 for (unsigned int ix = 0; ix <= this->degree; ++ix)
78 {
80
81 p[0] = points[ix][0];
82 if (codim > 1)
83 p[1] = points[iy][0];
84 if (codim > 2)
85 p[2] = points[iz][0];
86
87 this->unit_face_support_points[0][k++] = p;
88 }
89 AssertDimension(k, this->unit_face_support_points[0].size());
90 }
91
92 // initialize unit support points (this makes it possible to assign initial
93 // values to FE_FaceQ)
95 this->unit_face_support_points[0].size());
96 const unsigned int n_face_dofs = this->unit_face_support_points[0].size();
97 for (unsigned int i = 0; i < n_face_dofs; ++i)
98 for (unsigned int d = 0; d < dim; ++d)
99 {
100 for (unsigned int e = 0, c = 0; e < dim; ++e)
101 if (d != e)
102 {
103 // faces in y-direction are oriented differently
104 unsigned int renumber = i;
105 if (dim == 3 && d == 1)
106 renumber = i / (degree + 1) + (degree + 1) * (i % (degree + 1));
107 this->unit_support_points[n_face_dofs * 2 * d + i][e] =
108 this->unit_face_support_points[0][renumber][c];
109 this->unit_support_points[n_face_dofs * (2 * d + 1) + i][e] =
110 this->unit_face_support_points[0][renumber][c];
111 this->unit_support_points[n_face_dofs * (2 * d + 1) + i][d] = 1;
112 ++c;
113 }
114 }
115}
116
117
118
119template <int dim, int spacedim>
120std::unique_ptr<FiniteElement<dim, spacedim>>
122{
123 return std::make_unique<FE_FaceQ<dim, spacedim>>(this->degree);
124}
125
126
127
128template <int dim, int spacedim>
129std::string
131{
132 // note that the FETools::get_fe_by_name function depends on the
133 // particular format of the string this function returns, so they have to be
134 // kept in synch
135 std::ostringstream namebuf;
136 namebuf << "FE_FaceQ<" << Utilities::dim_string(dim, spacedim) << ">("
137 << this->degree << ")";
138
139 return namebuf.str();
140}
141
142
143
144template <int dim, int spacedim>
145void
147 const FiniteElement<dim, spacedim> &source_fe,
148 FullMatrix<double> &interpolation_matrix,
149 const unsigned int face_no) const
150{
151 get_subface_interpolation_matrix(source_fe,
153 interpolation_matrix,
154 face_no);
155}
156
157
158
159template <int dim, int spacedim>
160void
162 const FiniteElement<dim, spacedim> &x_source_fe,
163 const unsigned int subface,
164 FullMatrix<double> &interpolation_matrix,
165 const unsigned int face_no) const
166{
167 // this function is similar to the respective method in FE_Q
168
169 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
170 ExcDimensionMismatch(interpolation_matrix.n(),
171 this->n_dofs_per_face(face_no)));
172 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
173 ExcDimensionMismatch(interpolation_matrix.m(),
174 x_source_fe.n_dofs_per_face(face_no)));
175
176 // see if source is a FaceQ element
177 if (const FE_FaceQ<dim, spacedim> *source_fe =
178 dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&x_source_fe))
179 {
180 // Make sure that the element for which the DoFs should be constrained
181 // is the one with the higher polynomial degree. Actually the procedure
182 // will work also if this assertion is not satisfied. But the matrices
183 // produced in that case might lead to problems in the hp-procedures,
184 // which use this method.
185 Assert(
186 this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
187 (typename FiniteElement<dim,
188 spacedim>::ExcInterpolationNotImplemented()));
189
190 // generate a quadrature with the unit face support points.
191 const Quadrature<dim - 1> face_quadrature(
192 source_fe->get_unit_face_support_points(face_no));
193
194 // Rule of thumb for FP accuracy, that can be expected for a given
195 // polynomial degree. This value is used to cut off values close to
196 // zero.
197 const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
198
199 // compute the interpolation matrix by simply taking the value at the
200 // support points.
201 for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
202 {
203 const Point<dim - 1> p =
205 face_quadrature.point(i) :
207 face_quadrature.point(i), subface);
208
209 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
210 {
211 double matrix_entry = this->poly_space.compute_value(j, p);
212
213 // Correct the interpolated value. I.e. if it is close to 1 or 0,
214 // make it exactly 1 or 0. Unfortunately, this is required to
215 // avoid problems with higher order elements.
216 if (std::fabs(matrix_entry - 1.0) < eps)
217 matrix_entry = 1.0;
218 if (std::fabs(matrix_entry) < eps)
219 matrix_entry = 0.0;
220
221 interpolation_matrix(i, j) = matrix_entry;
222 }
223 }
224
225#ifdef DEBUG
226 // make sure that the row sum of each of the matrices is 1 at this
227 // point. this must be so since the shape functions sum up to 1
228 for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
229 {
230 double sum = 0.;
231
232 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
233 sum += interpolation_matrix(j, i);
234
235 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
236 }
237#endif
238 }
239 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
240 {
241 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
242 }
243 else
245 false,
246 (typename FiniteElement<dim,
247 spacedim>::ExcInterpolationNotImplemented()));
248}
249
250
251
252template <int dim, int spacedim>
253bool
255 const unsigned int shape_index,
256 const unsigned int face_index) const
257{
258 return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
259}
260
261
262
263template <int dim, int spacedim>
264std::vector<unsigned int>
266{
267 std::vector<unsigned int> dpo(dim + 1, 0U);
268 dpo[dim - 1] = deg + 1;
269 for (unsigned int i = 1; i < dim - 1; ++i)
270 dpo[dim - 1] *= deg + 1;
271 return dpo;
272}
273
274
275
276template <int dim, int spacedim>
277bool
282
283
284
285template <int dim, int spacedim>
286std::vector<std::pair<unsigned int, unsigned int>>
288 const FiniteElement<dim, spacedim> & /*fe_other*/) const
289{
290 // this element is always discontinuous at vertices
291 return std::vector<std::pair<unsigned int, unsigned int>>();
292}
293
294
295
296template <int dim, int spacedim>
297std::vector<std::pair<unsigned int, unsigned int>>
299 const FiniteElement<dim, spacedim> &fe_other) const
300{
301 Assert(dim >= 2, ExcInternalError());
302
303 // this element is continuous only for the highest dimensional bounding object
304 if (dim > 2)
305 return std::vector<std::pair<unsigned int, unsigned int>>();
306 else
307 {
308 // this is similar to the FE_Q_Base class
309 if (const FE_FaceQ<dim, spacedim> *fe_q_other =
310 dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
311 {
312 // dofs are located along lines, so two dofs are identical if they are
313 // located at identical positions.
314 // Therefore, read the points in unit_support_points for the
315 // first coordinate direction. We take the lexicographic ordering of
316 // the points in the second direction (i.e., y-direction) since we
317 // know that the first p+1 dofs are located at the left (x=0) face.
318 const unsigned int p = this->degree;
319 const unsigned int q = fe_q_other->degree;
320
321 std::vector<std::pair<unsigned int, unsigned int>> identities;
322
323 const std::vector<unsigned int> &index_map_inverse =
324 this->poly_space.get_numbering_inverse();
325 const std::vector<unsigned int> &index_map_inverse_other =
326 fe_q_other->poly_space.get_numbering_inverse();
327
328 for (unsigned int i = 0; i < p + 1; ++i)
329 for (unsigned int j = 0; j < q + 1; ++j)
330 if (std::fabs(
331 this->unit_support_points[index_map_inverse[i]][dim - 1] -
332 fe_q_other->unit_support_points[index_map_inverse_other[j]]
333 [dim - 1]) < 1e-14)
334 identities.emplace_back(i, j);
335
336 return identities;
337 }
338 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
339 {
340 // the FE_Nothing has no degrees of freedom, so there are no
341 // equivalencies to be recorded
342 return std::vector<std::pair<unsigned int, unsigned int>>();
343 }
344 else if (fe_other.n_unique_faces() == 1 &&
345 fe_other.n_dofs_per_face(0) == 0)
346 {
347 // if the other element has no elements on faces at all,
348 // then it would be impossible to enforce any kind of
349 // continuity even if we knew exactly what kind of element
350 // we have -- simply because the other element declares
351 // that it is discontinuous because it has no DoFs on
352 // its faces. in that case, just state that we have no
353 // constraints to declare
354 return std::vector<std::pair<unsigned int, unsigned int>>();
355 }
356 else
357 {
359 return std::vector<std::pair<unsigned int, unsigned int>>();
360 }
361 }
362}
363
364
365
366template <int dim, int spacedim>
367std::vector<std::pair<unsigned int, unsigned int>>
369 const FiniteElement<dim, spacedim> &fe_other,
370 const unsigned int) const
371{
372 Assert(dim >= 3, ExcInternalError());
373
374 // this element is continuous only for the highest dimensional bounding object
375 if (dim > 3)
376 return std::vector<std::pair<unsigned int, unsigned int>>();
377 else
378 {
379 // this is similar to the FE_Q_Base class
380 if (const FE_FaceQ<dim, spacedim> *fe_q_other =
381 dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
382 {
383 // this works exactly like the line case above, except that now we
384 // have to have two indices i1, i2 and j1, j2 to characterize the dofs
385 // on the face of each of the finite elements. since they are ordered
386 // lexicographically along the first line and we have a tensor
387 // product, the rest is rather straightforward
388 const unsigned int p = this->degree;
389 const unsigned int q = fe_q_other->degree;
390
391 std::vector<std::pair<unsigned int, unsigned int>> identities;
392
393 const std::vector<unsigned int> &index_map_inverse =
394 this->poly_space.get_numbering_inverse();
395 const std::vector<unsigned int> &index_map_inverse_other =
396 fe_q_other->poly_space.get_numbering_inverse();
397
398 std::vector<std::pair<unsigned int, unsigned int>> identities_1d;
399
400 for (unsigned int i = 0; i < p + 1; ++i)
401 for (unsigned int j = 0; j < q + 1; ++j)
402 if (std::fabs(
403 this->unit_support_points[index_map_inverse[i]][dim - 2] -
404 fe_q_other->unit_support_points[index_map_inverse_other[j]]
405 [dim - 2]) < 1e-14)
406 identities_1d.emplace_back(i, j);
407
408 for (unsigned int n1 = 0; n1 < identities_1d.size(); ++n1)
409 for (unsigned int n2 = 0; n2 < identities_1d.size(); ++n2)
410 identities.emplace_back(identities_1d[n1].first * (p + 1) +
411 identities_1d[n2].first,
412 identities_1d[n1].second * (q + 1) +
413 identities_1d[n2].second);
414
415 return identities;
416 }
417 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
418 {
419 // the FE_Nothing has no degrees of freedom, so there are no
420 // equivalencies to be recorded
421 return std::vector<std::pair<unsigned int, unsigned int>>();
422 }
423 else if (fe_other.n_unique_faces() == 1 &&
424 fe_other.n_dofs_per_face(0) == 0)
425 {
426 // if the other element has no elements on faces at all,
427 // then it would be impossible to enforce any kind of
428 // continuity even if we knew exactly what kind of element
429 // we have -- simply because the other element declares
430 // that it is discontinuous because it has no DoFs on
431 // its faces. in that case, just state that we have no
432 // constraints to declare
433 return std::vector<std::pair<unsigned int, unsigned int>>();
434 }
435 else
436 {
438 return std::vector<std::pair<unsigned int, unsigned int>>();
439 }
440 }
441}
442
443
444
445template <int dim, int spacedim>
448 const FiniteElement<dim, spacedim> &fe_other,
449 const unsigned int codim) const
450{
451 Assert(codim <= dim, ExcImpossibleInDim(dim));
452 (void)codim;
453
454 // vertex/line/face/cell domination
455 // --------------------------------
456 if (const FE_FaceQ<dim, spacedim> *fe_faceq_other =
457 dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
458 {
459 if (this->degree < fe_faceq_other->degree)
461 else if (this->degree == fe_faceq_other->degree)
463 else
465 }
466 else if (const FE_Nothing<dim> *fe_nothing =
467 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
468 {
469 if (fe_nothing->is_dominating())
471 else
472 // the FE_Nothing has no degrees of freedom and it is typically used
473 // in a context where we don't require any continuity along the
474 // interface
476 }
477
480}
481
482template <int dim, int spacedim>
483std::pair<Table<2, bool>, std::vector<unsigned int>>
485{
486 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
487 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
488 constant_modes(0, i) = true;
489 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
490 constant_modes, std::vector<unsigned int>(1, 0));
491}
492
493template <int dim, int spacedim>
494void
496 const std::vector<Vector<double>> &support_point_values,
497 std::vector<double> &nodal_values) const
498{
499 AssertDimension(support_point_values.size(),
500 this->get_unit_support_points().size());
501 AssertDimension(support_point_values.size(), nodal_values.size());
502 AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
503
504 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
505 {
506 AssertDimension(support_point_values[i].size(), 1);
507
508 nodal_values[i] = support_point_values[i](0);
509 }
510}
511
512// ----------------------------- FE_FaceQ<1,spacedim> ------------------------
513
514template <int spacedim>
515FE_FaceQ<1, spacedim>::FE_FaceQ(const unsigned int degree)
516 : FiniteElement<1, spacedim>(
517 FiniteElementData<1>(get_dpo_vector(degree),
518 1,
519 degree,
520 FiniteElementData<1>::L2),
521 std::vector<bool>(1, true),
522 std::vector<ComponentMask>(1, ComponentMask(1, true)))
523{
524 this->unit_face_support_points[0].resize(1);
525
526 // initialize unit support points (this makes it possible to assign initial
527 // values to FE_FaceQ)
529 this->unit_support_points[1] = Point<1>(1.);
530}
531
532
533
534template <int spacedim>
535std::unique_ptr<FiniteElement<1, spacedim>>
537{
538 return std::make_unique<FE_FaceQ<1, spacedim>>(this->degree);
539}
540
541
542
543template <int spacedim>
544std::string
546{
547 // note that the FETools::get_fe_by_name function depends on the
548 // particular format of the string this function returns, so they have to be
549 // kept in synch
550 std::ostringstream namebuf;
551 namebuf << "FE_FaceQ<" << Utilities::dim_string(1, spacedim) << ">("
552 << this->degree << ")";
553
554 return namebuf.str();
555}
556
557
558
559template <int spacedim>
560void
562 const FiniteElement<1, spacedim> &source_fe,
563 FullMatrix<double> &interpolation_matrix,
564 const unsigned int face_no) const
565{
566 get_subface_interpolation_matrix(source_fe,
568 interpolation_matrix,
569 face_no);
570}
571
572
573
574template <int spacedim>
575void
577 const FiniteElement<1, spacedim> &x_source_fe,
578 const unsigned int /*subface*/,
579 FullMatrix<double> &interpolation_matrix,
580 const unsigned int face_no) const
581{
582 (void)x_source_fe;
583 (void)face_no;
584
585 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
586 ExcDimensionMismatch(interpolation_matrix.n(),
587 this->n_dofs_per_face(face_no)));
588 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
589 ExcDimensionMismatch(interpolation_matrix.m(),
590 x_source_fe.n_dofs_per_face(face_no)));
591 interpolation_matrix(0, 0) = 1.;
592}
593
594
595
596template <int spacedim>
597bool
598FE_FaceQ<1, spacedim>::has_support_on_face(const unsigned int shape_index,
599 const unsigned int face_index) const
600{
601 AssertIndexRange(shape_index, 2);
602 return (face_index == shape_index);
603}
604
605
606
607template <int spacedim>
608std::vector<unsigned int>
610{
611 std::vector<unsigned int> dpo(2, 0U);
612 dpo[0] = 1;
613 return dpo;
614}
615
616
617
618template <int spacedim>
619bool
624
625template <int spacedim>
626std::vector<std::pair<unsigned int, unsigned int>>
628 const FiniteElement<1, spacedim> & /*fe_other*/) const
629{
630 // this element is always discontinuous at vertices
631 return std::vector<std::pair<unsigned int, unsigned int>>(1,
632 std::make_pair(0U,
633 0U));
634}
635
636
637
638template <int spacedim>
639std::vector<std::pair<unsigned int, unsigned int>>
641 const FiniteElement<1, spacedim> &) const
642{
643 // this element is continuous only for the highest dimensional bounding object
644 return std::vector<std::pair<unsigned int, unsigned int>>();
645}
646
647
648
649template <int spacedim>
650std::vector<std::pair<unsigned int, unsigned int>>
653 const unsigned int) const
654{
655 // this element is continuous only for the highest dimensional bounding object
656 return std::vector<std::pair<unsigned int, unsigned int>>();
657}
658
659
660
661template <int spacedim>
662std::pair<Table<2, bool>, std::vector<unsigned int>>
664{
665 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
666 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
667 constant_modes(0, i) = true;
668 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
669 constant_modes, std::vector<unsigned int>(1, 0));
670}
671
672
673
674template <int spacedim>
688
689
690template <int spacedim>
691void
695 const Quadrature<1> &,
696 const Mapping<1, spacedim> &,
701 spacedim>
702 &) const
703{
704 // Do nothing, since we do not have values in the interior
705}
706
707
708
709template <int spacedim>
710void
713 const unsigned int face,
714 const hp::QCollection<0> &,
715 const Mapping<1, spacedim> &,
718 const typename FiniteElement<1, spacedim>::InternalDataBase &fe_internal,
720 spacedim>
721 &output_data) const
722{
723 const unsigned int foffset = face;
724 if (fe_internal.update_each & update_values)
725 {
726 for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
727 output_data.shape_values(k, 0) = 0.;
728 output_data.shape_values(foffset, 0) = 1;
729 }
730}
731
732
733template <int spacedim>
734void
737 const unsigned int,
738 const unsigned int,
739 const Quadrature<0> &,
740 const Mapping<1, spacedim> &,
745 spacedim>
746 &) const
747{
748 Assert(false, ExcMessage("There are no sub-face values to fill in 1d!"));
749}
750
751
752
753// --------------------------------------- FE_FaceP --------------------------
754
755template <int dim, int spacedim>
756FE_FaceP<dim, spacedim>::FE_FaceP(const unsigned int degree)
757 : FE_PolyFace<PolynomialSpace<dim - 1>, dim, spacedim>(
758 PolynomialSpace<dim - 1>(
759 Polynomials::Legendre::generate_complete_basis(degree)),
760 FiniteElementData<dim>(get_dpo_vector(degree),
761 1,
762 degree,
763 FiniteElementData<dim>::L2),
764 std::vector<bool>(1, true))
765{}
766
767
768
769template <int dim, int spacedim>
770std::unique_ptr<FiniteElement<dim, spacedim>>
772{
773 return std::make_unique<FE_FaceP<dim, spacedim>>(this->degree);
774}
775
776
777
778template <int dim, int spacedim>
779std::string
781{
782 // note that the FETools::get_fe_by_name function depends on the
783 // particular format of the string this function returns, so they have to be
784 // kept in synch
785 std::ostringstream namebuf;
786 namebuf << "FE_FaceP<" << Utilities::dim_string(dim, spacedim) << ">("
787 << this->degree << ")";
788
789 return namebuf.str();
790}
791
792
793
794template <int dim, int spacedim>
795bool
797 const unsigned int shape_index,
798 const unsigned int face_index) const
799{
800 return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
801}
802
803
804
805template <int dim, int spacedim>
806std::vector<unsigned int>
808{
809 std::vector<unsigned int> dpo(dim + 1, 0U);
810 dpo[dim - 1] = deg + 1;
811 for (unsigned int i = 1; i < dim - 1; ++i)
812 {
813 dpo[dim - 1] *= deg + 1 + i;
814 dpo[dim - 1] /= i + 1;
815 }
816 return dpo;
817}
818
819
820
821template <int dim, int spacedim>
822bool
827
828
829
830template <int dim, int spacedim>
833 const FiniteElement<dim, spacedim> &fe_other,
834 const unsigned int codim) const
835{
836 Assert(codim <= dim, ExcImpossibleInDim(dim));
837 (void)codim;
838
839 // vertex/line/face/cell domination
840 // --------------------------------
841 if (const FE_FaceP<dim, spacedim> *fe_facep_other =
842 dynamic_cast<const FE_FaceP<dim, spacedim> *>(&fe_other))
843 {
844 if (this->degree < fe_facep_other->degree)
846 else if (this->degree == fe_facep_other->degree)
848 else
850 }
851 else if (const FE_Nothing<dim> *fe_nothing =
852 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
853 {
854 if (fe_nothing->is_dominating())
856 else
857 // the FE_Nothing has no degrees of freedom and it is typically used
858 // in a context where we don't require any continuity along the
859 // interface
861 }
862
865}
866
867
868
869template <int dim, int spacedim>
870void
872 const FiniteElement<dim, spacedim> &source_fe,
873 FullMatrix<double> &interpolation_matrix,
874 const unsigned int face_no) const
875{
876 get_subface_interpolation_matrix(source_fe,
878 interpolation_matrix,
879 face_no);
880}
881
882
883
884template <int dim, int spacedim>
885void
887 const FiniteElement<dim, spacedim> &x_source_fe,
888 const unsigned int subface,
889 FullMatrix<double> &interpolation_matrix,
890 const unsigned int face_no) const
891{
892 // this function is similar to the respective method in FE_Q
893
894 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
895 ExcDimensionMismatch(interpolation_matrix.n(),
896 this->n_dofs_per_face(face_no)));
897 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
898 ExcDimensionMismatch(interpolation_matrix.m(),
899 x_source_fe.n_dofs_per_face(face_no)));
900
901 // see if source is a FaceP element
902 if (const FE_FaceP<dim, spacedim> *source_fe =
903 dynamic_cast<const FE_FaceP<dim, spacedim> *>(&x_source_fe))
904 {
905 // Make sure that the element for which the DoFs should be constrained
906 // is the one with the higher polynomial degree. Actually the procedure
907 // will work also if this assertion is not satisfied. But the matrices
908 // produced in that case might lead to problems in the hp-procedures,
909 // which use this method.
910 Assert(
911 this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
912 (typename FiniteElement<dim,
913 spacedim>::ExcInterpolationNotImplemented()));
914
915 // do this as in FETools by solving a least squares problem where we
916 // force the source FE polynomial to be equal the given FE on all
917 // quadrature points
918 const QGauss<dim - 1> face_quadrature(source_fe->degree + 1);
919
920 // Rule of thumb for FP accuracy, that can be expected for a given
921 // polynomial degree. This value is used to cut off values close to
922 // zero.
923 const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
924
925 FullMatrix<double> mass(face_quadrature.size(),
926 source_fe->n_dofs_per_face(face_no));
927
928 for (unsigned int k = 0; k < face_quadrature.size(); ++k)
929 {
930 const Point<dim - 1> p =
932 face_quadrature.point(k) :
934 face_quadrature.point(k), subface);
935
936 for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
937 mass(k, j) = source_fe->poly_space.compute_value(j, p);
938 }
939
940 Householder<double> H(mass);
941 Vector<double> v_in(face_quadrature.size());
942 Vector<double> v_out(source_fe->n_dofs_per_face(face_no));
943
944
945 // compute the interpolation matrix by evaluating on the fine side and
946 // then solving the least squares problem
947 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
948 {
949 for (unsigned int k = 0; k < face_quadrature.size(); ++k)
950 {
951 const Point<dim - 1> p =
953 face_quadrature.point(k) :
955 face_quadrature.point(k), subface);
956 v_in(k) = this->poly_space.compute_value(i, p);
957 }
958 const double result = H.least_squares(v_out, v_in);
959 (void)result;
960 Assert(result < 1e-12, FETools::ExcLeastSquaresError(result));
961
962 for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
963 {
964 double matrix_entry = v_out(j);
965
966 // Correct the interpolated value. I.e. if it is close to 1 or 0,
967 // make it exactly 1 or 0. Unfortunately, this is required to
968 // avoid problems with higher order elements.
969 if (std::fabs(matrix_entry - 1.0) < eps)
970 matrix_entry = 1.0;
971 if (std::fabs(matrix_entry) < eps)
972 matrix_entry = 0.0;
973
974 interpolation_matrix(j, i) = matrix_entry;
975 }
976 }
977 }
978 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
979 {
980 // nothing to do here, the FE_Nothing has no degrees of freedom anyway
981 }
982 else
984 false,
985 (typename FiniteElement<dim,
986 spacedim>::ExcInterpolationNotImplemented()));
987}
988
989
990
991template <int dim, int spacedim>
992std::pair<Table<2, bool>, std::vector<unsigned int>>
994{
995 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
996 for (const unsigned int face : GeometryInfo<dim>::face_indices())
997 constant_modes(0, face * this->n_dofs_per_face(face)) = true;
998 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
999 constant_modes, std::vector<unsigned int>(1, 0));
1000}
1001
1002
1003
1004template <int spacedim>
1005FE_FaceP<1, spacedim>::FE_FaceP(const unsigned int degree)
1006 : FE_FaceQ<1, spacedim>(degree)
1007{}
1008
1009
1010
1011template <int spacedim>
1012std::string
1014{
1015 // note that the FETools::get_fe_by_name function depends on the
1016 // particular format of the string this function returns, so they have to be
1017 // kept in synch
1018 std::ostringstream namebuf;
1019 namebuf << "FE_FaceP<" << Utilities::dim_string(1, spacedim) << ">("
1020 << this->degree << ")";
1021
1022 return namebuf.str();
1023}
1024
1025
1026
1027// explicit instantiations
1028#include "fe_face.inst"
1029
1030
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_face.cc:993
virtual bool hp_constraints_are_implemented() const override
Definition fe_face.cc:823
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition fe_face.cc:796
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition fe_face.cc:832
virtual std::string get_name() const override
Definition fe_face.cc:780
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition fe_face.cc:807
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_face.cc:886
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_face.cc:771
FE_FaceP(unsigned int p)
Definition fe_face.cc:756
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_face.cc:871
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition fe_face.cc:447
virtual std::string get_name() const override
Definition fe_face.cc:130
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition fe_face.cc:265
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
Definition fe_face.cc:495
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_face.cc:484
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_face.cc:146
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition fe_face.cc:298
FE_FaceQ(const unsigned int p)
Definition fe_face.cc:50
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_face.cc:121
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_face.cc:161
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition fe_face.cc:254
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition fe_face.cc:287
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
Definition fe_face.cc:368
virtual bool hp_constraints_are_implemented() const override
Definition fe_face.cc:278
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
virtual void fill_fe_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const CellSimilarity::Similarity cell_similarity, const Quadrature< dim > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const=0
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition fe.h:2458
virtual void fill_fe_face_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const hp::QCollection< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const
virtual void fill_fe_subface_values(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const unsigned int face_no, const unsigned int sub_no, const Quadrature< dim - 1 > &quadrature, const Mapping< dim, spacedim > &mapping, const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, const internal::FEValuesImplementation::MappingRelatedData< dim, spacedim > &mapping_data, const InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const=0
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const=0
std::vector< Point< dim > > unit_support_points
Definition fe.h:2451
size_type n() const
size_type m() const
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
Point< 2 > second
Definition grid_out.cc:4624
#define Assert(cond, exc)
static ::ExceptionBase & ExcLeastSquaresError(double arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_covariant_transformation
Covariant transformation.
@ update_gradients
Shape function gradients.
#define DEAL_II_NOT_IMPLEMENTED()
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:555
constexpr T fixed_power(const T t)
Definition utilities.h:942
static const unsigned int invalid_unsigned_int
Definition types.h:220
STL namespace.
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static Point< dim > child_to_cell_coordinates(const Point< dim > &p, const unsigned int child_index, const RefinementCase< dim > refine_case=RefinementCase< dim >::isotropic_refinement)