Reference documentation for deal.II version GIT 9c182271f7 2023-03-28 14:30:01+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\}}\)
fe_face.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2022 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>
20 #include <deal.II/fe/fe_nothing.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 
32 namespace 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 
50 template <int dim, int spacedim>
51 FE_FaceQ<dim, spacedim>::FE_FaceQ(const unsigned int degree)
52  : FE_PolyFace<TensorProductPolynomials<dim - 1>, dim, spacedim>(
53  TensorProductPolynomials<dim - 1>(
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  {
80  Point<codim> p;
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 
120 template <int dim, int spacedim>
121 std::unique_ptr<FiniteElement<dim, spacedim>>
123 {
124  return std::make_unique<FE_FaceQ<dim, spacedim>>(this->degree);
125 }
126 
127 
128 
129 template <int dim, int spacedim>
130 std::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 
145 template <int dim, int spacedim>
146 void
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 
160 template <int dim, int spacedim>
161 void
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 =
205  subface == numbers::invalid_unsigned_int ?
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 
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
245  AssertThrow(
246  false,
247  (typename FiniteElement<dim,
248  spacedim>::ExcInterpolationNotImplemented()));
249 }
250 
251 
252 
253 template <int dim, int spacedim>
254 bool
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 
264 template <int dim, int spacedim>
265 std::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 
277 template <int dim, int spacedim>
278 bool
280 {
281  return true;
282 }
283 
284 
285 
286 template <int dim, int spacedim>
287 std::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 
297 template <int dim, int spacedim>
298 std::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 
367 template <int dim, int spacedim>
368 std::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 
446 template <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 
483 template <int dim, int spacedim>
484 std::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 
494 template <int dim, int spacedim>
495 void
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 
515 template <int spacedim>
516 FE_FaceQ<1, spacedim>::FE_FaceQ(const unsigned int degree)
517  : FiniteElement<1, spacedim>(
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 
535 template <int spacedim>
536 std::unique_ptr<FiniteElement<1, spacedim>>
538 {
539  return std::make_unique<FE_FaceQ<1, spacedim>>(this->degree);
540 }
541 
542 
543 
544 template <int spacedim>
545 std::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 
560 template <int spacedim>
561 void
563  const FiniteElement<1, spacedim> &source_fe,
564  FullMatrix<double> & interpolation_matrix,
565  const unsigned int face_no) const
566 {
569  interpolation_matrix,
570  face_no);
571 }
572 
573 
574 
575 template <int spacedim>
576 void
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 
597 template <int spacedim>
598 bool
599 FE_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 
608 template <int spacedim>
609 std::vector<unsigned int>
611 {
612  std::vector<unsigned int> dpo(2, 0U);
613  dpo[0] = 1;
614  return dpo;
615 }
616 
617 
618 
619 template <int spacedim>
620 bool
622 {
623  return true;
624 }
625 
626 template <int spacedim>
627 std::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 
639 template <int spacedim>
640 std::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 
650 template <int spacedim>
651 std::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 
662 template <int spacedim>
663 std::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 
675 template <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 
691 template <int spacedim>
692 void
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 
710 template <int spacedim>
711 void
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 
734 template <int spacedim>
735 void
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 
756 template <int dim, int spacedim>
757 FE_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 
770 template <int dim, int spacedim>
771 std::unique_ptr<FiniteElement<dim, spacedim>>
773 {
774  return std::make_unique<FE_FaceP<dim, spacedim>>(this->degree);
775 }
776 
777 
778 
779 template <int dim, int spacedim>
780 std::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 
795 template <int dim, int spacedim>
796 bool
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 
806 template <int dim, int spacedim>
807 std::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 
822 template <int dim, int spacedim>
823 bool
825 {
826  return true;
827 }
828 
829 
830 
831 template <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 
870 template <int dim, int spacedim>
871 void
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 
885 template <int dim, int spacedim>
886 void
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 =
932  subface == numbers::invalid_unsigned_int ?
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 =
953  subface == numbers::invalid_unsigned_int ?
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
984  AssertThrow(
985  false,
986  (typename FiniteElement<dim,
987  spacedim>::ExcInterpolationNotImplemented()));
988 }
989 
990 
991 
992 template <int dim, int spacedim>
993 std::pair<Table<2, bool>, std::vector<unsigned int>>
995 {
996  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
997  for (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 
1005 template <int spacedim>
1006 FE_FaceP<1, spacedim>::FE_FaceP(const unsigned int degree)
1007  : FE_FaceQ<1, spacedim>(degree)
1008 {}
1009 
1010 
1011 
1012 template <int spacedim>
1013 std::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 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 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 bool hp_constraints_are_implemented() const override
Definition: fe_face.cc:279
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 typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
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 typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
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 typename FiniteElement< dim, spacedim >::InternalDataBase &fe_internal, ::internal::FEValuesImplementation::FiniteElementRelatedData< dim, spacedim > &output_data) const override
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_cell() 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:2447
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:314
const std::vector< Point< dim > > & get_points() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
Point< 2 > second
Definition: grid_out.cc:4616
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcLeastSquaresError(double arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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.
Expression fabs(const Expression &x)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
static const char U
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:702
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:556
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
static const unsigned int invalid_unsigned_int
Definition: types.h:213
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)