Reference documentation for deal.II version 9.3.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\}}\)
fe_face.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2009 - 2020 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 {
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,
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  // make sure that the row sum of each of the matrices is 1 at this
227  // point. this must be so since the shape functions sum up to 1
228  for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
229  {
230  double sum = 0.;
231 
232  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
233  sum += interpolation_matrix(j, i);
234 
235  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
236  }
237  }
238  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
239  {
240  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
241  }
242  else
243  AssertThrow(
244  false,
245  (typename FiniteElement<dim,
246  spacedim>::ExcInterpolationNotImplemented()));
247 }
248 
249 
250 
251 template <int dim, int spacedim>
252 bool
254  const unsigned int shape_index,
255  const unsigned int face_index) const
256 {
257  return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
258 }
259 
260 
261 
262 template <int dim, int spacedim>
263 std::vector<unsigned int>
265 {
266  std::vector<unsigned int> dpo(dim + 1, 0U);
267  dpo[dim - 1] = deg + 1;
268  for (unsigned int i = 1; i < dim - 1; ++i)
269  dpo[dim - 1] *= deg + 1;
270  return dpo;
271 }
272 
273 
274 
275 template <int dim, int spacedim>
276 bool
278 {
279  return true;
280 }
281 
282 
283 
284 template <int dim, int spacedim>
285 std::vector<std::pair<unsigned int, unsigned int>>
287  const FiniteElement<dim, spacedim> & /*fe_other*/) const
288 {
289  // this element is always discontinuous at vertices
290  return std::vector<std::pair<unsigned int, unsigned int>>();
291 }
292 
293 
294 
295 template <int dim, int spacedim>
296 std::vector<std::pair<unsigned int, unsigned int>>
298  const FiniteElement<dim, spacedim> &fe_other) const
299 {
300  Assert(dim >= 2, ExcInternalError());
301 
302  // this element is continuous only for the highest dimensional bounding object
303  if (dim > 2)
304  return std::vector<std::pair<unsigned int, unsigned int>>();
305  else
306  {
307  // this is similar to the FE_Q_Base class
308  if (const FE_FaceQ<dim, spacedim> *fe_q_other =
309  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
310  {
311  // dofs are located along lines, so two dofs are identical if they are
312  // located at identical positions.
313  // Therefore, read the points in unit_support_points for the
314  // first coordinate direction. We take the lexicographic ordering of
315  // the points in the second direction (i.e., y-direction) since we
316  // know that the first p+1 dofs are located at the left (x=0) face.
317  const unsigned int p = this->degree;
318  const unsigned int q = fe_q_other->degree;
319 
320  std::vector<std::pair<unsigned int, unsigned int>> identities;
321 
322  const std::vector<unsigned int> &index_map_inverse =
324  const std::vector<unsigned int> &index_map_inverse_other =
325  fe_q_other->poly_space.get_numbering_inverse();
326 
327  for (unsigned int i = 0; i < p + 1; ++i)
328  for (unsigned int j = 0; j < q + 1; ++j)
329  if (std::fabs(
330  this->unit_support_points[index_map_inverse[i]][dim - 1] -
331  fe_q_other->unit_support_points[index_map_inverse_other[j]]
332  [dim - 1]) < 1e-14)
333  identities.emplace_back(i, j);
334 
335  return identities;
336  }
337  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
338  {
339  // the FE_Nothing has no degrees of freedom, so there are no
340  // equivalencies to be recorded
341  return std::vector<std::pair<unsigned int, unsigned int>>();
342  }
343  else if (fe_other.n_unique_faces() == 1 &&
344  fe_other.n_dofs_per_face(0) == 0)
345  {
346  // if the other element has no elements on faces at all,
347  // then it would be impossible to enforce any kind of
348  // continuity even if we knew exactly what kind of element
349  // we have -- simply because the other element declares
350  // that it is discontinuous because it has no DoFs on
351  // its faces. in that case, just state that we have no
352  // constraints to declare
353  return std::vector<std::pair<unsigned int, unsigned int>>();
354  }
355  else
356  {
357  Assert(false, ExcNotImplemented());
358  return std::vector<std::pair<unsigned int, unsigned int>>();
359  }
360  }
361 }
362 
363 
364 
365 template <int dim, int spacedim>
366 std::vector<std::pair<unsigned int, unsigned int>>
368  const FiniteElement<dim, spacedim> &fe_other,
369  const unsigned int) const
370 {
371  Assert(dim >= 3, ExcInternalError());
372 
373  // this element is continuous only for the highest dimensional bounding object
374  if (dim > 3)
375  return std::vector<std::pair<unsigned int, unsigned int>>();
376  else
377  {
378  // this is similar to the FE_Q_Base class
379  if (const FE_FaceQ<dim, spacedim> *fe_q_other =
380  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
381  {
382  // this works exactly like the line case above, except that now we
383  // have to have two indices i1, i2 and j1, j2 to characterize the dofs
384  // on the face of each of the finite elements. since they are ordered
385  // lexicographically along the first line and we have a tensor
386  // product, the rest is rather straightforward
387  const unsigned int p = this->degree;
388  const unsigned int q = fe_q_other->degree;
389 
390  std::vector<std::pair<unsigned int, unsigned int>> identities;
391 
392  const std::vector<unsigned int> &index_map_inverse =
394  const std::vector<unsigned int> &index_map_inverse_other =
395  fe_q_other->poly_space.get_numbering_inverse();
396 
397  std::vector<std::pair<unsigned int, unsigned int>> identities_1d;
398 
399  for (unsigned int i = 0; i < p + 1; ++i)
400  for (unsigned int j = 0; j < q + 1; ++j)
401  if (std::fabs(
402  this->unit_support_points[index_map_inverse[i]][dim - 2] -
403  fe_q_other->unit_support_points[index_map_inverse_other[j]]
404  [dim - 2]) < 1e-14)
405  identities_1d.emplace_back(i, j);
406 
407  for (unsigned int n1 = 0; n1 < identities_1d.size(); ++n1)
408  for (unsigned int n2 = 0; n2 < identities_1d.size(); ++n2)
409  identities.emplace_back(identities_1d[n1].first * (p + 1) +
410  identities_1d[n2].first,
411  identities_1d[n1].second * (q + 1) +
412  identities_1d[n2].second);
413 
414  return identities;
415  }
416  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
417  {
418  // the FE_Nothing has no degrees of freedom, so there are no
419  // equivalencies to be recorded
420  return std::vector<std::pair<unsigned int, unsigned int>>();
421  }
422  else if (fe_other.n_unique_faces() == 1 &&
423  fe_other.n_dofs_per_face(0) == 0)
424  {
425  // if the other element has no elements on faces at all,
426  // then it would be impossible to enforce any kind of
427  // continuity even if we knew exactly what kind of element
428  // we have -- simply because the other element declares
429  // that it is discontinuous because it has no DoFs on
430  // its faces. in that case, just state that we have no
431  // constraints to declare
432  return std::vector<std::pair<unsigned int, unsigned int>>();
433  }
434  else
435  {
436  Assert(false, ExcNotImplemented());
437  return std::vector<std::pair<unsigned int, unsigned int>>();
438  }
439  }
440 }
441 
442 
443 
444 template <int dim, int spacedim>
447  const FiniteElement<dim, spacedim> &fe_other,
448  const unsigned int codim) const
449 {
450  Assert(codim <= dim, ExcImpossibleInDim(dim));
451  (void)codim;
452 
453  // vertex/line/face/cell domination
454  // --------------------------------
455  if (const FE_FaceQ<dim, spacedim> *fe_faceq_other =
456  dynamic_cast<const FE_FaceQ<dim, spacedim> *>(&fe_other))
457  {
458  if (this->degree < fe_faceq_other->degree)
460  else if (this->degree == fe_faceq_other->degree)
462  else
464  }
465  else if (const FE_Nothing<dim> *fe_nothing =
466  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
467  {
468  if (fe_nothing->is_dominating())
470  else
471  // the FE_Nothing has no degrees of freedom and it is typically used
472  // in a context where we don't require any continuity along the
473  // interface
475  }
476 
477  Assert(false, ExcNotImplemented());
479 }
480 
481 template <int dim, int spacedim>
482 std::pair<Table<2, bool>, std::vector<unsigned int>>
484 {
485  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
486  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
487  constant_modes(0, i) = true;
488  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
489  constant_modes, std::vector<unsigned int>(1, 0));
490 }
491 
492 template <int dim, int spacedim>
493 void
495  const std::vector<Vector<double>> &support_point_values,
496  std::vector<double> & nodal_values) const
497 {
498  AssertDimension(support_point_values.size(),
499  this->get_unit_support_points().size());
500  AssertDimension(support_point_values.size(), nodal_values.size());
501  AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
502 
503  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
504  {
505  AssertDimension(support_point_values[i].size(), 1);
506 
507  nodal_values[i] = support_point_values[i](0);
508  }
509 }
510 
511 // ----------------------------- FE_FaceQ<1,spacedim> ------------------------
512 
513 template <int spacedim>
515  : FiniteElement<1, spacedim>(
517  1,
518  degree,
519  FiniteElementData<1>::L2),
520  std::vector<bool>(1, true),
521  std::vector<ComponentMask>(1, ComponentMask(1, true)))
522 {
523  this->unit_face_support_points[0].resize(1);
524 
525  // initialize unit support points (this makes it possible to assign initial
526  // values to FE_FaceQ)
528  this->unit_support_points[1] = Point<1>(1.);
529 }
530 
531 
532 
533 template <int spacedim>
534 std::unique_ptr<FiniteElement<1, spacedim>>
536 {
537  return std::make_unique<FE_FaceQ<1, spacedim>>(this->degree);
538 }
539 
540 
541 
542 template <int spacedim>
543 std::string
545 {
546  // note that the FETools::get_fe_by_name function depends on the
547  // particular format of the string this function returns, so they have to be
548  // kept in synch
549  std::ostringstream namebuf;
550  namebuf << "FE_FaceQ<" << Utilities::dim_string(1, spacedim) << ">("
551  << this->degree << ")";
552 
553  return namebuf.str();
554 }
555 
556 
557 
558 template <int spacedim>
559 void
561  const FiniteElement<1, spacedim> &source_fe,
562  FullMatrix<double> & interpolation_matrix,
563  const unsigned int face_no) const
564 {
567  interpolation_matrix,
568  face_no);
569 }
570 
571 
572 
573 template <int spacedim>
574 void
576  const FiniteElement<1, spacedim> &x_source_fe,
577  const unsigned int /*subface*/,
578  FullMatrix<double> &interpolation_matrix,
579  const unsigned int face_no) const
580 {
581  (void)x_source_fe;
582  (void)face_no;
583 
584  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
585  ExcDimensionMismatch(interpolation_matrix.n(),
586  this->n_dofs_per_face(face_no)));
587  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
588  ExcDimensionMismatch(interpolation_matrix.m(),
589  x_source_fe.n_dofs_per_face(face_no)));
590  interpolation_matrix(0, 0) = 1.;
591 }
592 
593 
594 
595 template <int spacedim>
596 bool
597 FE_FaceQ<1, spacedim>::has_support_on_face(const unsigned int shape_index,
598  const unsigned int face_index) const
599 {
600  AssertIndexRange(shape_index, 2);
601  return (face_index == shape_index);
602 }
603 
604 
605 
606 template <int spacedim>
607 std::vector<unsigned int>
609 {
610  std::vector<unsigned int> dpo(2, 0U);
611  dpo[0] = 1;
612  return dpo;
613 }
614 
615 
616 
617 template <int spacedim>
618 bool
620 {
621  return true;
622 }
623 
624 template <int spacedim>
625 std::vector<std::pair<unsigned int, unsigned int>>
627  const FiniteElement<1, spacedim> & /*fe_other*/) const
628 {
629  // this element is always discontinuous at vertices
630  return std::vector<std::pair<unsigned int, unsigned int>>(1,
631  std::make_pair(0U,
632  0U));
633 }
634 
635 
636 
637 template <int spacedim>
638 std::vector<std::pair<unsigned int, unsigned int>>
640  const FiniteElement<1, spacedim> &) const
641 {
642  // this element is continuous only for the highest dimensional bounding object
643  return std::vector<std::pair<unsigned int, unsigned int>>();
644 }
645 
646 
647 
648 template <int spacedim>
649 std::vector<std::pair<unsigned int, unsigned int>>
652  const unsigned int) const
653 {
654  // this element is continuous only for the highest dimensional bounding object
655  return std::vector<std::pair<unsigned int, unsigned int>>();
656 }
657 
658 
659 
660 template <int spacedim>
661 std::pair<Table<2, bool>, std::vector<unsigned int>>
663 {
664  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
665  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
666  constant_modes(0, i) = true;
667  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
668  constant_modes, std::vector<unsigned int>(1, 0));
669 }
670 
671 
672 
673 template <int spacedim>
676 {
677  UpdateFlags out = flags & update_values;
678  if (flags & update_gradients)
679  out |= update_gradients | update_covariant_transformation;
680  if (flags & update_hessians)
681  out |= update_hessians | update_covariant_transformation;
682  if (flags & update_normal_vectors)
683  out |= update_normal_vectors | update_JxW_values;
684 
685  return out;
686 }
687 
688 
689 template <int spacedim>
690 void
694  const Quadrature<1> &,
695  const Mapping<1, spacedim> &,
697  const ::internal::FEValuesImplementation::MappingRelatedData<1,
698  spacedim>
699  &,
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> &,
718  const ::internal::FEValuesImplementation::MappingRelatedData<1,
719  spacedim>
720  &,
721  const typename FiniteElement<1, spacedim>::InternalDataBase &fe_internal,
723  spacedim>
724  &output_data) const
725 {
726  const unsigned int foffset = face;
727  if (fe_internal.update_each & update_values)
728  {
729  for (unsigned int k = 0; k < this->n_dofs_per_cell(); ++k)
730  output_data.shape_values(k, 0) = 0.;
731  output_data.shape_values(foffset, 0) = 1;
732  }
733 }
734 
735 
736 template <int spacedim>
737 void
740  const unsigned int,
741  const unsigned int,
742  const Quadrature<0> &,
743  const Mapping<1, spacedim> &,
745  const ::internal::FEValuesImplementation::MappingRelatedData<1,
746  spacedim>
747  &,
750  spacedim>
751  &) const
752 {
753  Assert(false, ExcMessage("There are no sub-face values to fill in 1D!"));
754 }
755 
756 
757 
758 // --------------------------------------- FE_FaceP --------------------------
759 
760 template <int dim, int spacedim>
762  : FE_PolyFace<PolynomialSpace<dim - 1>, dim, spacedim>(
763  PolynomialSpace<dim - 1>(
764  Polynomials::Legendre::generate_complete_basis(degree)),
765  FiniteElementData<dim>(get_dpo_vector(degree),
766  1,
767  degree,
768  FiniteElementData<dim>::L2),
769  std::vector<bool>(1, true))
770 {}
771 
772 
773 
774 template <int dim, int spacedim>
775 std::unique_ptr<FiniteElement<dim, spacedim>>
777 {
778  return std::make_unique<FE_FaceP<dim, spacedim>>(this->degree);
779 }
780 
781 
782 
783 template <int dim, int spacedim>
784 std::string
786 {
787  // note that the FETools::get_fe_by_name function depends on the
788  // particular format of the string this function returns, so they have to be
789  // kept in synch
790  std::ostringstream namebuf;
791  namebuf << "FE_FaceP<" << Utilities::dim_string(dim, spacedim) << ">("
792  << this->degree << ")";
793 
794  return namebuf.str();
795 }
796 
797 
798 
799 template <int dim, int spacedim>
800 bool
802  const unsigned int shape_index,
803  const unsigned int face_index) const
804 {
805  return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
806 }
807 
808 
809 
810 template <int dim, int spacedim>
811 std::vector<unsigned int>
813 {
814  std::vector<unsigned int> dpo(dim + 1, 0U);
815  dpo[dim - 1] = deg + 1;
816  for (unsigned int i = 1; i < dim - 1; ++i)
817  {
818  dpo[dim - 1] *= deg + 1 + i;
819  dpo[dim - 1] /= i + 1;
820  }
821  return dpo;
822 }
823 
824 
825 
826 template <int dim, int spacedim>
827 bool
829 {
830  return true;
831 }
832 
833 
834 
835 template <int dim, int spacedim>
838  const FiniteElement<dim, spacedim> &fe_other,
839  const unsigned int codim) const
840 {
841  Assert(codim <= dim, ExcImpossibleInDim(dim));
842  (void)codim;
843 
844  // vertex/line/face/cell domination
845  // --------------------------------
846  if (const FE_FaceP<dim, spacedim> *fe_facep_other =
847  dynamic_cast<const FE_FaceP<dim, spacedim> *>(&fe_other))
848  {
849  if (this->degree < fe_facep_other->degree)
851  else if (this->degree == fe_facep_other->degree)
853  else
855  }
856  else if (const FE_Nothing<dim> *fe_nothing =
857  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
858  {
859  if (fe_nothing->is_dominating())
861  else
862  // the FE_Nothing has no degrees of freedom and it is typically used
863  // in a context where we don't require any continuity along the
864  // interface
866  }
867 
868  Assert(false, ExcNotImplemented());
870 }
871 
872 
873 
874 template <int dim, int spacedim>
875 void
877  const FiniteElement<dim, spacedim> &source_fe,
878  FullMatrix<double> & interpolation_matrix,
879  const unsigned int face_no) const
880 {
883  interpolation_matrix,
884  face_no);
885 }
886 
887 
888 
889 template <int dim, int spacedim>
890 void
892  const FiniteElement<dim, spacedim> &x_source_fe,
893  const unsigned int subface,
894  FullMatrix<double> & interpolation_matrix,
895  const unsigned int face_no) const
896 {
897  // this function is similar to the respective method in FE_Q
898 
899  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
900  ExcDimensionMismatch(interpolation_matrix.n(),
901  this->n_dofs_per_face(face_no)));
902  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
903  ExcDimensionMismatch(interpolation_matrix.m(),
904  x_source_fe.n_dofs_per_face(face_no)));
905 
906  // see if source is a FaceP element
907  if (const FE_FaceP<dim, spacedim> *source_fe =
908  dynamic_cast<const FE_FaceP<dim, spacedim> *>(&x_source_fe))
909  {
910  // Make sure that the element for which the DoFs should be constrained
911  // is the one with the higher polynomial degree. Actually the procedure
912  // will work also if this assertion is not satisfied. But the matrices
913  // produced in that case might lead to problems in the hp-procedures,
914  // which use this method.
915  Assert(
916  this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
917  (typename FiniteElement<dim,
919 
920  // do this as in FETools by solving a least squares problem where we
921  // force the source FE polynomial to be equal the given FE on all
922  // quadrature points
923  const QGauss<dim - 1> face_quadrature(source_fe->degree + 1);
924 
925  // Rule of thumb for FP accuracy, that can be expected for a given
926  // polynomial degree. This value is used to cut off values close to
927  // zero.
928  const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
929 
930  FullMatrix<double> mass(face_quadrature.size(),
931  source_fe->n_dofs_per_face(face_no));
932 
933  for (unsigned int k = 0; k < face_quadrature.size(); ++k)
934  {
935  const Point<dim - 1> p =
936  subface == numbers::invalid_unsigned_int ?
937  face_quadrature.point(k) :
939  face_quadrature.point(k), subface);
940 
941  for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
942  mass(k, j) = source_fe->poly_space.compute_value(j, p);
943  }
944 
945  Householder<double> H(mass);
946  Vector<double> v_in(face_quadrature.size());
947  Vector<double> v_out(source_fe->n_dofs_per_face(face_no));
948 
949 
950  // compute the interpolation matrix by evaluating on the fine side and
951  // then solving the least squares problem
952  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
953  {
954  for (unsigned int k = 0; k < face_quadrature.size(); ++k)
955  {
956  const Point<dim - 1> p =
958  face_quadrature.point(k) :
960  face_quadrature.point(k), subface);
961  v_in(k) = this->poly_space.compute_value(i, p);
962  }
963  const double result = H.least_squares(v_out, v_in);
964  (void)result;
965  Assert(result < 1e-12, FETools::ExcLeastSquaresError(result));
966 
967  for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
968  {
969  double matrix_entry = v_out(j);
970 
971  // Correct the interpolated value. I.e. if it is close to 1 or 0,
972  // make it exactly 1 or 0. Unfortunately, this is required to
973  // avoid problems with higher order elements.
974  if (std::fabs(matrix_entry - 1.0) < eps)
975  matrix_entry = 1.0;
976  if (std::fabs(matrix_entry) < eps)
977  matrix_entry = 0.0;
978 
979  interpolation_matrix(j, i) = matrix_entry;
980  }
981  }
982  }
983  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
984  {
985  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
986  }
987  else
988  AssertThrow(
989  false,
990  (typename FiniteElement<dim,
991  spacedim>::ExcInterpolationNotImplemented()));
992 }
993 
994 
995 
996 template <int dim, int spacedim>
997 std::pair<Table<2, bool>, std::vector<unsigned int>>
999 {
1000  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1001  for (unsigned int face : GeometryInfo<dim>::face_indices())
1002  constant_modes(0, face * this->n_dofs_per_face(face)) = true;
1003  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1004  constant_modes, std::vector<unsigned int>(1, 0));
1005 }
1006 
1007 
1008 
1009 template <int spacedim>
1011  : FE_FaceQ<1, spacedim>(degree)
1012 {}
1013 
1014 
1015 
1016 template <int spacedim>
1017 std::string
1019 {
1020  // note that the FETools::get_fe_by_name function depends on the
1021  // particular format of the string this function returns, so they have to be
1022  // kept in synch
1023  std::ostringstream namebuf;
1024  namebuf << "FE_FaceP<" << Utilities::dim_string(1, spacedim) << ">("
1025  << this->degree << ")";
1026 
1027  return namebuf.str();
1028 }
1029 
1030 
1031 
1032 // explicit instantiations
1033 #include "fe_face.inst"
1034 
1035 
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
static const unsigned int invalid_unsigned_int
Definition: types.h:196
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:483
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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:876
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
const unsigned int degree
Definition: fe_base.h:435
STL namespace.
static const char U
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_face.cc:253
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_face.cc:776
static ::ExceptionBase & ExcInterpolationNotImplemented()
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_face.cc:122
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2442
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:367
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
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
T sum(const T &t, const MPI_Comm &mpi_communicator)
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_face.cc:446
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:297
#define Assert(cond, exc)
Definition: exceptions.h:1465
UpdateFlags
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Abstract base class for mapping classes.
Definition: mapping.h:303
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
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2449
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:286
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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 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 bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_face.cc:801
virtual std::string get_name() const override
Definition: fe_face.cc:131
Expression fabs(const Expression &x)
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition: fe_face.cc:812
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
virtual bool hp_constraints_are_implemented() const override
Definition: fe_face.cc:828
Second derivatives of shape functions.
unsigned int n_unique_faces() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:558
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:1049
double compute_value(const unsigned int i, const Point< dim > &p) const override
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:891
double compute_value(const unsigned int i, const Point< dim > &p) const override
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
Shape function gradients.
Normal vectors.
static std::vector< unsigned int > get_dpo_vector(const unsigned int deg)
Definition: fe_face.cc:264
const std::vector< unsigned int > & get_numbering_inverse() const
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:998
static ::ExceptionBase & ExcNotImplemented()
FE_FaceP(unsigned int p)
Definition: fe_face.cc:761
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_face.cc:837
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
virtual std::string get_name() const override
Definition: fe_face.cc:785
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:494
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:701
Covariant transformation.
static ::ExceptionBase & ExcInternalError()
virtual bool hp_constraints_are_implemented() const override
Definition: fe_face.cc:277