Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_face.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2018 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 
17 #include <deal.II/base/quadrature_lib.h>
18 #include <deal.II/base/std_cxx14/memory.h>
19 
20 #include <deal.II/fe/fe_face.h>
21 #include <deal.II/fe/fe_nothing.h>
22 #include <deal.II/fe/fe_poly_face.templates.h>
23 #include <deal.II/fe/fe_tools.h>
24 
25 #include <deal.II/lac/householder.h>
26 
27 #include <sstream>
28 
29 DEAL_II_NAMESPACE_OPEN
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>(
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.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][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[k++] = p;
89  }
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.size());
97  const unsigned int n_face_dofs = this->unit_face_support_points.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[renumber][c];
110  this->unit_support_points[n_face_dofs * (2 * d + 1) + i][e] =
111  this->unit_face_support_points[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_cxx14::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) const
150 {
151  get_subface_interpolation_matrix(source_fe,
153  interpolation_matrix);
154 }
155 
156 
157 
158 template <int dim, int spacedim>
159 void
161  const FiniteElement<dim, spacedim> &x_source_fe,
162  const unsigned int subface,
163  FullMatrix<double> & interpolation_matrix) const
164 {
165  // this function is similar to the respective method in FE_Q
166 
167  Assert(interpolation_matrix.n() == this->dofs_per_face,
168  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
169  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
170  ExcDimensionMismatch(interpolation_matrix.m(),
171  x_source_fe.dofs_per_face));
172 
173  // see if source is a FaceQ element
174  if (const FE_FaceQ<dim, spacedim> *source_fe =
175  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&x_source_fe))
176  {
177  // Make sure that the element for which the DoFs should be constrained
178  // is the one with the higher polynomial degree. Actually the procedure
179  // will work also if this assertion is not satisfied. But the matrices
180  // produced in that case might lead to problems in the hp procedures,
181  // which use this method.
182  Assert(
183  this->dofs_per_face <= source_fe->dofs_per_face,
184  (typename FiniteElement<dim,
185  spacedim>::ExcInterpolationNotImplemented()));
186 
187  // generate a quadrature with the unit face support points.
188  const Quadrature<dim - 1> face_quadrature(
189  source_fe->get_unit_face_support_points());
190 
191  // Rule of thumb for FP accuracy, that can be expected for a given
192  // polynomial degree. This value is used to cut off values close to
193  // zero.
194  const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
195 
196  // compute the interpolation matrix by simply taking the value at the
197  // support points.
198  for (unsigned int i = 0; i < source_fe->dofs_per_face; ++i)
199  {
200  const Point<dim - 1> p =
201  subface == numbers::invalid_unsigned_int ?
202  face_quadrature.point(i) :
204  face_quadrature.point(i), subface);
205 
206  for (unsigned int j = 0; j < this->dofs_per_face; ++j)
207  {
208  double matrix_entry = this->poly_space.compute_value(j, p);
209 
210  // Correct the interpolated value. I.e. if it is close to 1 or 0,
211  // make it exactly 1 or 0. Unfortunately, this is required to
212  // avoid problems with higher order elements.
213  if (std::fabs(matrix_entry - 1.0) < eps)
214  matrix_entry = 1.0;
215  if (std::fabs(matrix_entry) < eps)
216  matrix_entry = 0.0;
217 
218  interpolation_matrix(i, j) = matrix_entry;
219  }
220  }
221 
222  // make sure that the row sum of each of the matrices is 1 at this
223  // point. this must be so since the shape functions sum up to 1
224  for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
225  {
226  double sum = 0.;
227 
228  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
229  sum += interpolation_matrix(j, i);
230 
231  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
232  }
233  }
234  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
235  {
236  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
237  }
238  else
239  AssertThrow(
240  false,
241  (typename FiniteElement<dim,
242  spacedim>::ExcInterpolationNotImplemented()));
243 }
244 
245 
246 
247 template <int dim, int spacedim>
248 bool
250  const unsigned int shape_index,
251  const unsigned int face_index) const
252 {
253  return (face_index == (shape_index / this->dofs_per_face));
254 }
255 
256 
257 
258 template <int dim, int spacedim>
259 std::vector<unsigned int>
261 {
262  std::vector<unsigned int> dpo(dim + 1, 0U);
263  dpo[dim - 1] = deg + 1;
264  for (unsigned int i = 1; i < dim - 1; ++i)
265  dpo[dim - 1] *= deg + 1;
266  return dpo;
267 }
268 
269 
270 
271 template <int dim, int spacedim>
272 bool
274 {
275  return true;
276 }
277 
278 
279 
280 template <int dim, int spacedim>
281 std::vector<std::pair<unsigned int, unsigned int>>
283  const FiniteElement<dim, spacedim> & /*fe_other*/) const
284 {
285  // this element is always discontinuous at vertices
286  return std::vector<std::pair<unsigned int, unsigned int>>();
287 }
288 
289 
290 
291 template <int dim, int spacedim>
292 std::vector<std::pair<unsigned int, unsigned int>>
294  const FiniteElement<dim, spacedim> &fe_other) const
295 {
296  Assert(dim >= 2, ExcInternalError());
297 
298  // this element is continuous only for the highest dimensional bounding object
299  if (dim > 2)
300  return std::vector<std::pair<unsigned int, unsigned int>>();
301  else
302  {
303  // this is similar to the FE_Q_Base class
304  if (const FE_FaceQ<dim, spacedim> *fe_q_other =
305  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
306  {
307  // dofs are located along lines, so two dofs are identical if they are
308  // located at identical positions.
309  // Therefore, read the points in unit_support_points for the
310  // first coordinate direction. We take the lexicographic ordering of
311  // the points in the second direction (i.e., y-direction) since we
312  // know that the first p+1 dofs are located at the left (x=0) face.
313  const unsigned int p = this->degree;
314  const unsigned int q = fe_q_other->degree;
315 
316  std::vector<std::pair<unsigned int, unsigned int>> identities;
317 
318  const std::vector<unsigned int> &index_map_inverse =
319  this->poly_space.get_numbering_inverse();
320  const std::vector<unsigned int> &index_map_inverse_other =
321  fe_q_other->poly_space.get_numbering_inverse();
322 
323  for (unsigned int i = 0; i < p + 1; ++i)
324  for (unsigned int j = 0; j < q + 1; ++j)
325  if (std::fabs(
326  this->unit_support_points[index_map_inverse[i]][dim - 1] -
327  fe_q_other->unit_support_points[index_map_inverse_other[j]]
328  [dim - 1]) < 1e-14)
329  identities.emplace_back(i, j);
330 
331  return identities;
332  }
333  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
334  {
335  // the FE_Nothing has no degrees of freedom, so there are no
336  // equivalencies to be recorded
337  return std::vector<std::pair<unsigned int, unsigned int>>();
338  }
339  else if (fe_other.dofs_per_face == 0)
340  {
341  // if the other element has no elements on faces at all,
342  // then it would be impossible to enforce any kind of
343  // continuity even if we knew exactly what kind of element
344  // we have -- simply because the other element declares
345  // that it is discontinuous because it has no DoFs on
346  // its faces. in that case, just state that we have no
347  // constraints to declare
348  return std::vector<std::pair<unsigned int, unsigned int>>();
349  }
350  else
351  {
352  Assert(false, ExcNotImplemented());
353  return std::vector<std::pair<unsigned int, unsigned int>>();
354  }
355  }
356 }
357 
358 
359 
360 template <int dim, int spacedim>
361 std::vector<std::pair<unsigned int, unsigned int>>
363  const FiniteElement<dim, spacedim> &fe_other) const
364 {
365  Assert(dim >= 3, ExcInternalError());
366 
367  // this element is continuous only for the highest dimensional bounding object
368  if (dim > 3)
369  return std::vector<std::pair<unsigned int, unsigned int>>();
370  else
371  {
372  // this is similar to the FE_Q_Base class
373  if (const FE_FaceQ<dim, spacedim> *fe_q_other =
374  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
375  {
376  // this works exactly like the line case above, except that now we
377  // have to have two indices i1, i2 and j1, j2 to characterize the dofs
378  // on the face of each of the finite elements. since they are ordered
379  // lexicographically along the first line and we have a tensor
380  // product, the rest is rather straightforward
381  const unsigned int p = this->degree;
382  const unsigned int q = fe_q_other->degree;
383 
384  std::vector<std::pair<unsigned int, unsigned int>> identities;
385 
386  const std::vector<unsigned int> &index_map_inverse =
387  this->poly_space.get_numbering_inverse();
388  const std::vector<unsigned int> &index_map_inverse_other =
389  fe_q_other->poly_space.get_numbering_inverse();
390 
391  std::vector<std::pair<unsigned int, unsigned int>> identities_1d;
392 
393  for (unsigned int i = 0; i < p + 1; ++i)
394  for (unsigned int j = 0; j < q + 1; ++j)
395  if (std::fabs(
396  this->unit_support_points[index_map_inverse[i]][dim - 2] -
397  fe_q_other->unit_support_points[index_map_inverse_other[j]]
398  [dim - 2]) < 1e-14)
399  identities_1d.emplace_back(i, j);
400 
401  for (unsigned int n1 = 0; n1 < identities_1d.size(); ++n1)
402  for (unsigned int n2 = 0; n2 < identities_1d.size(); ++n2)
403  identities.emplace_back(identities_1d[n1].first * (p + 1) +
404  identities_1d[n2].first,
405  identities_1d[n1].second * (q + 1) +
406  identities_1d[n2].second);
407 
408  return identities;
409  }
410  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
411  {
412  // the FE_Nothing has no degrees of freedom, so there are no
413  // equivalencies to be recorded
414  return std::vector<std::pair<unsigned int, unsigned int>>();
415  }
416  else if (fe_other.dofs_per_face == 0)
417  {
418  // if the other element has no elements on faces at all,
419  // then it would be impossible to enforce any kind of
420  // continuity even if we knew exactly what kind of element
421  // we have -- simply because the other element declares
422  // that it is discontinuous because it has no DoFs on
423  // its faces. in that case, just state that we have no
424  // constraints to declare
425  return std::vector<std::pair<unsigned int, unsigned int>>();
426  }
427  else
428  {
429  Assert(false, ExcNotImplemented());
430  return std::vector<std::pair<unsigned int, unsigned int>>();
431  }
432  }
433 }
434 
435 
436 
437 template <int dim, int spacedim>
440  const FiniteElement<dim, spacedim> &fe_other,
441  const unsigned int codim) const
442 {
443  Assert(codim <= dim, ExcImpossibleInDim(dim));
444  (void)codim;
445 
446  // vertex/line/face/cell domination
447  // --------------------------------
448  if (const FE_FaceQ<dim, spacedim> *fe_faceq_other =
449  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
450  {
451  if (this->degree < fe_faceq_other->degree)
453  else if (this->degree == fe_faceq_other->degree)
455  else
457  }
458  else if (const FE_Nothing<dim> *fe_nothing =
459  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
460  {
461  if (fe_nothing->is_dominating())
463  else
464  // the FE_Nothing has no degrees of freedom and it is typically used
465  // in a context where we don't require any continuity along the
466  // interface
468  }
469 
470  Assert(false, ExcNotImplemented());
472 }
473 
474 template <int dim, int spacedim>
475 std::pair<Table<2, bool>, std::vector<unsigned int>>
477 {
478  Table<2, bool> constant_modes(1, this->dofs_per_cell);
479  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
480  constant_modes(0, i) = true;
481  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
482  constant_modes, std::vector<unsigned int>(1, 0));
483 }
484 
485 template <int dim, int spacedim>
486 void
488  const std::vector<Vector<double>> &support_point_values,
489  std::vector<double> & nodal_values) const
490 {
491  AssertDimension(support_point_values.size(),
492  this->get_unit_support_points().size());
493  AssertDimension(support_point_values.size(), nodal_values.size());
494  AssertDimension(this->dofs_per_cell, nodal_values.size());
495 
496  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
497  {
498  AssertDimension(support_point_values[i].size(), 1);
499 
500  nodal_values[i] = support_point_values[i](0);
501  }
502 }
503 
504 // ----------------------------- FE_FaceQ<1,spacedim> ------------------------
505 
506 template <int spacedim>
507 FE_FaceQ<1, spacedim>::FE_FaceQ(const unsigned int degree)
508  : FiniteElement<1, spacedim>(
509  FiniteElementData<1>(get_dpo_vector(degree),
510  1,
511  degree,
512  FiniteElementData<1>::L2),
513  std::vector<bool>(1, true),
514  std::vector<ComponentMask>(1, ComponentMask(1, true)))
515 {
516  this->unit_face_support_points.resize(1);
517 
518  // initialize unit support points (this makes it possible to assign initial
519  // values to FE_FaceQ)
521  this->unit_support_points[1] = Point<1>(1.);
522 }
523 
524 
525 
526 template <int spacedim>
527 std::unique_ptr<FiniteElement<1, spacedim>>
529 {
530  return std_cxx14::make_unique<FE_FaceQ<1, spacedim>>(this->degree);
531 }
532 
533 
534 
535 template <int spacedim>
536 std::string
538 {
539  // note that the FETools::get_fe_by_name function depends on the
540  // particular format of the string this function returns, so they have to be
541  // kept in synch
542  std::ostringstream namebuf;
543  namebuf << "FE_FaceQ<" << Utilities::dim_string(1, spacedim) << ">("
544  << this->degree << ")";
545 
546  return namebuf.str();
547 }
548 
549 
550 
551 template <int spacedim>
552 void
554  const FiniteElement<1, spacedim> &source_fe,
555  FullMatrix<double> & interpolation_matrix) const
556 {
559  interpolation_matrix);
560 }
561 
562 
563 
564 template <int spacedim>
565 void
567  const FiniteElement<1, spacedim> &x_source_fe,
568  const unsigned int /*subface*/,
569  FullMatrix<double> &interpolation_matrix) const
570 {
571  (void)x_source_fe;
572  Assert(interpolation_matrix.n() == this->dofs_per_face,
573  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
574  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
575  ExcDimensionMismatch(interpolation_matrix.m(),
576  x_source_fe.dofs_per_face));
577  interpolation_matrix(0, 0) = 1.;
578 }
579 
580 
581 
582 template <int spacedim>
583 bool
584 FE_FaceQ<1, spacedim>::has_support_on_face(const unsigned int shape_index,
585  const unsigned int face_index) const
586 {
587  AssertIndexRange(shape_index, 2);
588  return (face_index == shape_index);
589 }
590 
591 
592 
593 template <int spacedim>
594 std::vector<unsigned int>
596 {
597  std::vector<unsigned int> dpo(2, 0U);
598  dpo[0] = 1;
599  return dpo;
600 }
601 
602 
603 
604 template <int spacedim>
605 bool
607 {
608  return true;
609 }
610 
611 template <int spacedim>
612 std::vector<std::pair<unsigned int, unsigned int>>
614  const FiniteElement<1, spacedim> & /*fe_other*/) const
615 {
616  // this element is always discontinuous at vertices
617  return std::vector<std::pair<unsigned int, unsigned int>>(1,
618  std::make_pair(0U,
619  0U));
620 }
621 
622 
623 
624 template <int spacedim>
625 std::vector<std::pair<unsigned int, unsigned int>>
627  const FiniteElement<1, spacedim> &) const
628 {
629  // this element is continuous only for the highest dimensional bounding object
630  return std::vector<std::pair<unsigned int, unsigned int>>();
631 }
632 
633 
634 
635 template <int spacedim>
636 std::vector<std::pair<unsigned int, unsigned int>>
638  const FiniteElement<1, spacedim> &) const
639 {
640  // this element is continuous only for the highest dimensional bounding object
641  return std::vector<std::pair<unsigned int, unsigned int>>();
642 }
643 
644 
645 
646 template <int spacedim>
647 std::pair<Table<2, bool>, std::vector<unsigned int>>
649 {
650  Table<2, bool> constant_modes(1, this->dofs_per_cell);
651  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
652  constant_modes(0, i) = true;
653  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
654  constant_modes, std::vector<unsigned int>(1, 0));
655 }
656 
657 
658 
659 template <int spacedim>
662 {
663  UpdateFlags out = flags & update_values;
664  if (flags & update_gradients)
666  if (flags & update_hessians)
668  if (flags & update_normal_vectors)
670 
671  return out;
672 }
673 
674 
675 template <int spacedim>
676 void
680  const Quadrature<1> &,
681  const Mapping<1, spacedim> &,
683  const ::internal::FEValuesImplementation::MappingRelatedData<1,
684  spacedim>
685  &,
688  spacedim>
689  &) const
690 {
691  // Do nothing, since we do not have values in the interior
692 }
693 
694 
695 
696 template <int spacedim>
697 void
700  const unsigned int face,
701  const Quadrature<0> &,
702  const Mapping<1, spacedim> &,
704  const ::internal::FEValuesImplementation::MappingRelatedData<1,
705  spacedim>
706  &,
707  const typename FiniteElement<1, spacedim>::InternalDataBase &fe_internal,
709  spacedim>
710  &output_data) const
711 {
712  const unsigned int foffset = face;
713  if (fe_internal.update_each & update_values)
714  {
715  for (unsigned int k = 0; k < this->dofs_per_cell; ++k)
716  output_data.shape_values(k, 0) = 0.;
717  output_data.shape_values(foffset, 0) = 1;
718  }
719 }
720 
721 
722 template <int spacedim>
723 void
726  const unsigned int,
727  const unsigned int,
728  const Quadrature<0> &,
729  const Mapping<1, spacedim> &,
731  const ::internal::FEValuesImplementation::MappingRelatedData<1,
732  spacedim>
733  &,
736  spacedim>
737  &) const
738 {
739  Assert(false, ExcMessage("There are no sub-face values to fill in 1D!"));
740 }
741 
742 
743 
744 // --------------------------------------- FE_FaceP --------------------------
745 
746 template <int dim, int spacedim>
747 FE_FaceP<dim, spacedim>::FE_FaceP(const unsigned int degree)
748  : FE_PolyFace<PolynomialSpace<dim - 1>, dim, spacedim>(
749  PolynomialSpace<dim - 1>(
750  Polynomials::Legendre::generate_complete_basis(degree)),
751  FiniteElementData<dim>(get_dpo_vector(degree),
752  1,
753  degree,
754  FiniteElementData<dim>::L2),
755  std::vector<bool>(1, true))
756 {}
757 
758 
759 
760 template <int dim, int spacedim>
761 std::unique_ptr<FiniteElement<dim, spacedim>>
763 {
764  return std_cxx14::make_unique<FE_FaceP<dim, spacedim>>(this->degree);
765 }
766 
767 
768 
769 template <int dim, int spacedim>
770 std::string
772 {
773  // note that the FETools::get_fe_by_name function depends on the
774  // particular format of the string this function returns, so they have to be
775  // kept in synch
776  std::ostringstream namebuf;
777  namebuf << "FE_FaceP<" << Utilities::dim_string(dim, spacedim) << ">("
778  << this->degree << ")";
779 
780  return namebuf.str();
781 }
782 
783 
784 
785 template <int dim, int spacedim>
786 bool
788  const unsigned int shape_index,
789  const unsigned int face_index) const
790 {
791  return (face_index == (shape_index / this->dofs_per_face));
792 }
793 
794 
795 
796 template <int dim, int spacedim>
797 std::vector<unsigned int>
799 {
800  std::vector<unsigned int> dpo(dim + 1, 0U);
801  dpo[dim - 1] = deg + 1;
802  for (unsigned int i = 1; i < dim - 1; ++i)
803  {
804  dpo[dim - 1] *= deg + 1 + i;
805  dpo[dim - 1] /= i + 1;
806  }
807  return dpo;
808 }
809 
810 
811 
812 template <int dim, int spacedim>
813 bool
815 {
816  return true;
817 }
818 
819 
820 
821 template <int dim, int spacedim>
824  const FiniteElement<dim, spacedim> &fe_other,
825  const unsigned int codim) const
826 {
827  Assert(codim <= dim, ExcImpossibleInDim(dim));
828  (void)codim;
829 
830  // vertex/line/face/cell domination
831  // --------------------------------
832  if (const FE_FaceP<dim, spacedim> *fe_facep_other =
833  dynamic_cast<const FE_FaceP<dim, spacedim> *>(&fe_other))
834  {
835  if (this->degree < fe_facep_other->degree)
837  else if (this->degree == fe_facep_other->degree)
839  else
841  }
842  else if (const FE_Nothing<dim> *fe_nothing =
843  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
844  {
845  if (fe_nothing->is_dominating())
847  else
848  // the FE_Nothing has no degrees of freedom and it is typically used
849  // in a context where we don't require any continuity along the
850  // interface
852  }
853 
854  Assert(false, ExcNotImplemented());
856 }
857 
858 
859 
860 template <int dim, int spacedim>
861 void
863  const FiniteElement<dim, spacedim> &source_fe,
864  FullMatrix<double> & interpolation_matrix) const
865 {
866  get_subface_interpolation_matrix(source_fe,
868  interpolation_matrix);
869 }
870 
871 
872 
873 template <int dim, int spacedim>
874 void
876  const FiniteElement<dim, spacedim> &x_source_fe,
877  const unsigned int subface,
878  FullMatrix<double> & interpolation_matrix) const
879 {
880  // this function is similar to the respective method in FE_Q
881 
882  Assert(interpolation_matrix.n() == this->dofs_per_face,
883  ExcDimensionMismatch(interpolation_matrix.n(), this->dofs_per_face));
884  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
885  ExcDimensionMismatch(interpolation_matrix.m(),
886  x_source_fe.dofs_per_face));
887 
888  // see if source is a FaceP element
889  if (const FE_FaceP<dim, spacedim> *source_fe =
890  dynamic_cast<const FE_FaceP<dim, spacedim> *>(&x_source_fe))
891  {
892  // Make sure that the element for which the DoFs should be constrained
893  // is the one with the higher polynomial degree. Actually the procedure
894  // will work also if this assertion is not satisfied. But the matrices
895  // produced in that case might lead to problems in the hp procedures,
896  // which use this method.
897  Assert(
898  this->dofs_per_face <= source_fe->dofs_per_face,
899  (typename FiniteElement<dim,
900  spacedim>::ExcInterpolationNotImplemented()));
901 
902  // do this as in FETools by solving a least squares problem where we
903  // force the source FE polynomial to be equal the given FE on all
904  // quadrature points
905  const QGauss<dim - 1> face_quadrature(source_fe->degree + 1);
906 
907  // Rule of thumb for FP accuracy, that can be expected for a given
908  // polynomial degree. This value is used to cut off values close to
909  // zero.
910  const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
911 
912  FullMatrix<double> mass(face_quadrature.size(), source_fe->dofs_per_face);
913 
914  for (unsigned int k = 0; k < face_quadrature.size(); ++k)
915  {
916  const Point<dim - 1> p =
917  subface == numbers::invalid_unsigned_int ?
918  face_quadrature.point(k) :
920  face_quadrature.point(k), subface);
921 
922  for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
923  mass(k, j) = source_fe->poly_space.compute_value(j, p);
924  }
925 
926  Householder<double> H(mass);
927  Vector<double> v_in(face_quadrature.size());
928  Vector<double> v_out(source_fe->dofs_per_face);
929 
930 
931  // compute the interpolation matrix by evaluating on the fine side and
932  // then solving the least squares problem
933  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
934  {
935  for (unsigned int k = 0; k < face_quadrature.size(); ++k)
936  {
937  const Point<dim - 1> p =
939  face_quadrature.point(k) :
941  face_quadrature.point(k), subface);
942  v_in(k) = this->poly_space.compute_value(i, p);
943  }
944  const double result = H.least_squares(v_out, v_in);
945  (void)result;
946  Assert(result < 1e-12, FETools::ExcLeastSquaresError(result));
947 
948  for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
949  {
950  double matrix_entry = v_out(j);
951 
952  // Correct the interpolated value. I.e. if it is close to 1 or 0,
953  // make it exactly 1 or 0. Unfortunately, this is required to
954  // avoid problems with higher order elements.
955  if (std::fabs(matrix_entry - 1.0) < eps)
956  matrix_entry = 1.0;
957  if (std::fabs(matrix_entry) < eps)
958  matrix_entry = 0.0;
959 
960  interpolation_matrix(j, i) = matrix_entry;
961  }
962  }
963  }
964  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
965  {
966  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
967  }
968  else
969  AssertThrow(
970  false,
971  (typename FiniteElement<dim,
972  spacedim>::ExcInterpolationNotImplemented()));
973 }
974 
975 
976 
977 template <int dim, int spacedim>
978 std::pair<Table<2, bool>, std::vector<unsigned int>>
980 {
981  Table<2, bool> constant_modes(1, this->dofs_per_cell);
982  for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
983  constant_modes(0, face * this->dofs_per_face) = true;
984  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
985  constant_modes, std::vector<unsigned int>(1, 0));
986 }
987 
988 
989 
990 template <int spacedim>
991 FE_FaceP<1, spacedim>::FE_FaceP(const unsigned int degree)
992  : FE_FaceQ<1, spacedim>(degree)
993 {}
994 
995 
996 
997 template <int spacedim>
998 std::string
1000 {
1001  // note that the FETools::get_fe_by_name function depends on the
1002  // particular format of the string this function returns, so they have to be
1003  // kept in synch
1004  std::ostringstream namebuf;
1005  namebuf << "FE_FaceP<" << Utilities::dim_string(1, spacedim) << ">("
1006  << this->degree << ")";
1007 
1008  return namebuf.str();
1009 }
1010 
1011 
1012 
1013 // explicit instantiations
1014 #include "fe_face.inst"
1015 
1016 
1017 DEAL_II_NAMESPACE_CLOSE
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)
Transformed quadrature weights.
Shape function values.
size_type m() const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_face.cc:862
static const unsigned int invalid_unsigned_int
Definition: types.h:173
static ::ExceptionBase & ExcLeastSquaresError(double arg1)
virtual UpdateFlags requires_update_flags(const UpdateFlags update_flags) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_face.cc:476
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
const unsigned int degree
Definition: fe_base.h:298
std::vector< Point< dim - 1 > > unit_face_support_points
Definition: fe.h:2487
STL namespace.
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_face.cc:249
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_face.cc:762
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_face.cc:122
Definition: point.h:110
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2480
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_face.cc:875
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_face.cc:439
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:293
#define Assert(cond, exc)
Definition: exceptions.h:1407
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: dof_tools.h:57
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_face.cc:160
FE_FaceQ(const unsigned int p)
Definition: fe_face.cc:51
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:282
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_face.cc:787
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:798
virtual bool hp_constraints_are_implemented() const override
Definition: fe_face.cc:814
Second derivatives of shape functions.
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:457
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_face.cc:362
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
const unsigned int dofs_per_cell
Definition: fe_base.h:282
Shape function gradients.
Normal vectors.
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition: fe_face.cc:260
const unsigned int dofs_per_face
Definition: fe_base.h:275
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_face.cc:979
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:37
FE_FaceP(unsigned int p)
Definition: fe_face.cc:747
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_face.cc:147
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_face.cc:823
virtual std::string get_name() const override
Definition: fe_face.cc:771
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:487
Covariant transformation.
static ::ExceptionBase & ExcInternalError()
virtual bool hp_constraints_are_implemented() const override
Definition: fe_face.cc:273