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