Reference documentation for deal.II version Git 70360f6955 2021-06-22 11:20:06 -0600
\(\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_simplex_p.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2021 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 #include <deal.II/base/config.h>
17 
20 
21 #include <deal.II/fe/fe_dgq.h>
22 #include <deal.II/fe/fe_nothing.h>
23 #include <deal.II/fe/fe_q.h>
25 #include <deal.II/fe/fe_tools.h>
26 
28 
29 namespace
30 {
35  std::vector<unsigned int>
36  get_dpo_vector_fe_p(const unsigned int dim, const unsigned int degree)
37  {
38  switch (dim)
39  {
40  case 1:
41  switch (degree)
42  {
43  case 1:
44  return {1, 0};
45  case 2:
46  return {1, 1};
47  default:
48  Assert(false, ExcNotImplemented());
49  }
50  case 2:
51  switch (degree)
52  {
53  case 1:
54  return {1, 0, 0};
55  case 2:
56  return {1, 1, 0};
57  default:
58  Assert(false, ExcNotImplemented());
59  }
60  case 3:
61  switch (degree)
62  {
63  case 1:
64  return {1, 0, 0, 0};
65  case 2:
66  return {1, 1, 0, 0};
67  default:
68  Assert(false, ExcNotImplemented());
69  }
70  }
71 
72  Assert(false, ExcNotImplemented());
73  return {};
74  }
75 
76 
77 
82  template <int dim>
83  std::vector<Point<dim>>
84  unit_support_points_fe_poly(const unsigned int degree)
85  {
86  std::vector<Point<dim>> unit_points;
87 
88  // Piecewise constants are a special case: use a support point at the
89  // centroid and only the centroid
90  if (degree == 0)
91  {
92  Point<dim> centroid;
93  std::fill(centroid.begin_raw(),
94  centroid.end_raw(),
95  1.0 / double(dim + 1));
96  unit_points.emplace_back(centroid);
97  return unit_points;
98  }
99 
100  if (dim == 1)
101  {
102  // We don't really have dim = 1 support for simplex elements yet, but
103  // its convenient for populating the face array
104  Assert(degree <= 2, ExcNotImplemented());
105  if (degree >= 1)
106  {
107  unit_points.emplace_back(0.0);
108  unit_points.emplace_back(1.0);
109 
110  if (degree == 2)
111  unit_points.emplace_back(0.5);
112  }
113  }
114  else if (dim == 2)
115  {
116  Assert(degree <= 2, ExcNotImplemented());
117  if (degree >= 1)
118  {
119  unit_points.emplace_back(0.0, 0.0);
120  unit_points.emplace_back(1.0, 0.0);
121  unit_points.emplace_back(0.0, 1.0);
122 
123  if (degree == 2)
124  {
125  unit_points.emplace_back(0.5, 0.0);
126  unit_points.emplace_back(0.5, 0.5);
127  unit_points.emplace_back(0.0, 0.5);
128  }
129  }
130  }
131  else if (dim == 3)
132  {
133  Assert(degree <= 2, ExcNotImplemented());
134  if (degree >= 1)
135  {
136  unit_points.emplace_back(0.0, 0.0, 0.0);
137  unit_points.emplace_back(1.0, 0.0, 0.0);
138  unit_points.emplace_back(0.0, 1.0, 0.0);
139  unit_points.emplace_back(0.0, 0.0, 1.0);
140 
141  if (degree == 2)
142  {
143  unit_points.emplace_back(0.5, 0.0, 0.0);
144  unit_points.emplace_back(0.5, 0.5, 0.0);
145  unit_points.emplace_back(0.0, 0.5, 0.0);
146  unit_points.emplace_back(0.0, 0.0, 0.5);
147  unit_points.emplace_back(0.5, 0.0, 0.5);
148  unit_points.emplace_back(0.0, 0.5, 0.5);
149  }
150  }
151  }
152  else
153  {
154  Assert(false, ExcNotImplemented());
155  }
156 
157  return unit_points;
158  }
159 
164  template <int dim>
165  std::vector<std::vector<Point<dim - 1>>>
166  unit_face_support_points_fe_poly(const unsigned int degree)
167  {
168  // this concept doesn't exist in 1D so just return an empty vector
169  if (dim == 1)
170  return {};
171 
172  std::vector<std::vector<Point<dim - 1>>> unit_face_points;
173 
174  // all faces have the same support points
175  for (auto face_n :
177  .face_indices())
178  {
179  (void)face_n;
180  unit_face_points.emplace_back(
181  unit_support_points_fe_poly<dim - 1>(degree));
182  }
183 
184  return unit_face_points;
185  }
186 
192  template <int dim>
194  constraints_fe_poly(const unsigned int /*degree*/)
195  {
196  // no constraints in 1d
197  // constraints in 3d not implemented yet
198  return FullMatrix<double>();
199  }
200 
201  template <>
203  constraints_fe_poly<2>(const unsigned int degree)
204  {
205  const unsigned int dim = 2;
206 
207  Assert(degree <= 2, ExcNotImplemented());
208 
209  // the following implements the 2d case
210  // (the 3d case is not implemented yet)
211  //
212  // consult FE_Q_Base::Implementation::initialize_constraints()
213  // for more information
214 
215  std::vector<Point<dim - 1>> constraint_points;
216  // midpoint
217  constraint_points.emplace_back(0.5);
218  if (degree == 2)
219  {
220  // midpoint on subface 0
221  constraint_points.emplace_back(0.25);
222  // midpoint on subface 1
223  constraint_points.emplace_back(0.75);
224  }
225 
226  // Now construct relation between destination (child) and source (mother)
227  // dofs.
228 
229  const unsigned int n_dofs_constrained = constraint_points.size();
230  unsigned int n_dofs_per_face = degree + 1;
231  FullMatrix<double> interface_constraints(n_dofs_constrained,
232  n_dofs_per_face);
233 
234  const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
235 
236  for (unsigned int i = 0; i < n_dofs_constrained; ++i)
237  for (unsigned int j = 0; j < n_dofs_per_face; ++j)
238  {
239  interface_constraints(i, j) =
240  poly.compute_value(j, constraint_points[i]);
241 
242  // if the value is small up to round-off, then simply set it to zero
243  // to avoid unwanted fill-in of the constraint matrices (which would
244  // then increase the number of other DoFs a constrained DoF would
245  // couple to)
246  if (std::fabs(interface_constraints(i, j)) < 1e-13)
247  interface_constraints(i, j) = 0;
248  }
249  return interface_constraints;
250  }
251 
252 
253 
258  std::vector<unsigned int>
259  get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
260  {
261  // First treat the case of piecewise constant elements:
262  if (degree == 0)
263  {
264  std::vector<unsigned int> dpo(dim + 1, 0U);
265  dpo[dim] = 1;
266  return dpo;
267  }
268  else
269  {
270  // This element has the same degrees of freedom as the continuous one,
271  // but they are all counted for the interior of the cell because
272  // it is continuous. Rather than hard-code how many DoFs the element
273  // has, we just get the numbers from the continuous case and add them
274  // up
275  const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
276 
277  switch (dim)
278  {
279  case 1:
280  return {0U,
281  ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
282  continuous_dpo[dim]};
283 
284  case 2:
285  return {0U,
286  0U,
288  continuous_dpo[0] +
289  ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
290  continuous_dpo[dim]};
291 
292  case 3:
293  return {
294  0U,
295  0U,
296  0U,
297  ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
298  ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
299  ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
300  continuous_dpo[dim]};
301  }
302 
303  Assert(false, ExcNotImplemented());
304  return {};
305  }
306  }
307 } // namespace
308 
309 
310 
311 template <int dim, int spacedim>
313  const unsigned int degree,
314  const std::vector<unsigned int> & dpo_vector,
315  const typename FiniteElementData<dim>::Conformity conformity)
316  : ::FE_Poly<dim, spacedim>(
317  BarycentricPolynomials<dim>::get_fe_p_basis(degree),
318  FiniteElementData<dim>(dpo_vector,
319  dim == 2 ? ReferenceCells::Triangle :
321  1,
322  degree,
323  conformity),
324  std::vector<bool>(FiniteElementData<dim>(dpo_vector,
325  dim == 2 ?
328  1,
329  degree)
330  .dofs_per_cell,
331  true),
332  std::vector<ComponentMask>(
333  FiniteElementData<dim>(dpo_vector,
334  dim == 2 ? ReferenceCells::Triangle :
336  1,
337  degree)
338  .dofs_per_cell,
339  std::vector<bool>(1, true)))
340 {
341  this->unit_support_points = unit_support_points_fe_poly<dim>(degree);
342  // Discontinuous elements don't have face support points
343  if (conformity == FiniteElementData<dim>::Conformity::H1)
345  unit_face_support_points_fe_poly<dim>(degree);
346  this->interface_constraints = constraints_fe_poly<dim>(degree);
347 }
348 
349 
350 
351 template <int dim, int spacedim>
352 std::pair<Table<2, bool>, std::vector<unsigned int>>
354 {
355  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
356  constant_modes.fill(true);
357  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
358  constant_modes, std::vector<unsigned int>(1, 0));
359 }
360 
361 
362 
363 template <int dim, int spacedim>
364 const FullMatrix<double> &
366  const unsigned int child,
367  const RefinementCase<dim> &refinement_case) const
368 {
371  AssertDimension(dim, spacedim);
372 
373  // initialization upon first request
374  if (this->prolongation[refinement_case - 1][child].n() == 0)
375  {
376  std::lock_guard<std::mutex> lock(this->mutex);
377 
378  // if matrix got updated while waiting for the lock
379  if (this->prolongation[refinement_case - 1][child].n() ==
380  this->n_dofs_per_cell())
381  return this->prolongation[refinement_case - 1][child];
382 
383  // now do the work. need to get a non-const version of data in order to
384  // be able to modify them inside a const function
385  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
386 
387  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
389  isotropic_matrices.back().resize(
392 
393  FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
394 
395  this_nonconst.prolongation[refinement_case - 1].swap(
396  isotropic_matrices.back());
397  }
398 
399  // finally return the matrix
400  return this->prolongation[refinement_case - 1][child];
401 }
402 
403 
404 
405 template <int dim, int spacedim>
406 const FullMatrix<double> &
408  const unsigned int child,
409  const RefinementCase<dim> &refinement_case) const
410 {
413  AssertDimension(dim, spacedim);
414 
415  // initialization upon first request
416  if (this->restriction[refinement_case - 1][child].n() == 0)
417  {
418  std::lock_guard<std::mutex> lock(this->mutex);
419 
420  // if matrix got updated while waiting for the lock
421  if (this->restriction[refinement_case - 1][child].n() ==
422  this->n_dofs_per_cell())
423  return this->restriction[refinement_case - 1][child];
424 
425  // now do the work. need to get a non-const version of data in order to
426  // be able to modify them inside a const function
427  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
428 
429  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
431  isotropic_matrices.back().resize(
434 
435  FETools::compute_projection_matrices(*this, isotropic_matrices, true);
436 
437  this_nonconst.restriction[refinement_case - 1].swap(
438  isotropic_matrices.back());
439  }
440 
441  // finally return the matrix
442  return this->restriction[refinement_case - 1][child];
443 }
444 
445 
446 
447 template <int dim, int spacedim>
448 void
450  const FiniteElement<dim, spacedim> &source_fe,
451  FullMatrix<double> & interpolation_matrix,
452  const unsigned int face_no) const
453 {
454  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
455  ExcDimensionMismatch(interpolation_matrix.m(),
456  source_fe.n_dofs_per_face(face_no)));
457 
458  // see if source is a P or Q element
459  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
460  nullptr) ||
461  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
462  {
463  const Quadrature<dim - 1> quad_face_support(
464  source_fe.get_unit_face_support_points(face_no));
465 
466  const double eps = 2e-13 * this->degree * (dim - 1);
467 
468  std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
470  quad_face_support,
471  face_no,
472  face_quadrature_points);
473 
474  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
475  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
476  {
477  double matrix_entry =
478  this->shape_value(this->face_to_cell_index(j, 0),
479  face_quadrature_points[i]);
480 
481  // Correct the interpolated value. I.e. if it is close to 1 or
482  // 0, make it exactly 1 or 0. Unfortunately, this is required to
483  // avoid problems with higher order elements.
484  if (std::fabs(matrix_entry - 1.0) < eps)
485  matrix_entry = 1.0;
486  if (std::fabs(matrix_entry) < eps)
487  matrix_entry = 0.0;
488 
489  interpolation_matrix(i, j) = matrix_entry;
490  }
491 
492 #ifdef DEBUG
493  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
494  {
495  double sum = 0.;
496 
497  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
498  sum += interpolation_matrix(j, i);
499 
500  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
501  }
502 #endif
503  }
504  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
505  {
506  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
507  }
508  else
509  AssertThrow(
510  false,
511  (typename FiniteElement<dim,
512  spacedim>::ExcInterpolationNotImplemented()));
513 }
514 
515 
516 
517 template <int dim, int spacedim>
518 void
520  const FiniteElement<dim, spacedim> &source_fe,
521  const unsigned int subface,
522  FullMatrix<double> & interpolation_matrix,
523  const unsigned int face_no) const
524 {
525  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
526  ExcDimensionMismatch(interpolation_matrix.m(),
527  source_fe.n_dofs_per_face(face_no)));
528 
529  // see if source is a P or Q element
530  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
531  nullptr) ||
532  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
533  {
534  const Quadrature<dim - 1> quad_face_support(
535  source_fe.get_unit_face_support_points(face_no));
536 
537  const double eps = 2e-13 * this->degree * (dim - 1);
538 
539  std::vector<Point<dim>> subface_quadrature_points(
540  quad_face_support.size());
542  quad_face_support,
543  face_no,
544  subface,
545  subface_quadrature_points);
546 
547  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
548  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
549  {
550  double matrix_entry =
551  this->shape_value(this->face_to_cell_index(j, 0),
552  subface_quadrature_points[i]);
553 
554  // Correct the interpolated value. I.e. if it is close to 1 or
555  // 0, make it exactly 1 or 0. Unfortunately, this is required to
556  // avoid problems with higher order elements.
557  if (std::fabs(matrix_entry - 1.0) < eps)
558  matrix_entry = 1.0;
559  if (std::fabs(matrix_entry) < eps)
560  matrix_entry = 0.0;
561 
562  interpolation_matrix(i, j) = matrix_entry;
563  }
564 
565 #ifdef DEBUG
566  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
567  {
568  double sum = 0.;
569 
570  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
571  sum += interpolation_matrix(j, i);
572 
573  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
574  }
575 #endif
576  }
577  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
578  {
579  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
580  }
581  else
582  AssertThrow(
583  false,
584  (typename FiniteElement<dim,
585  spacedim>::ExcInterpolationNotImplemented()));
586 }
587 
588 
589 
590 template <int dim, int spacedim>
591 bool
593 {
594  return true;
595 }
596 
597 
598 
599 template <int dim, int spacedim>
600 void
603  const std::vector<Vector<double>> &support_point_values,
604  std::vector<double> & nodal_values) const
605 {
606  AssertDimension(support_point_values.size(),
607  this->get_unit_support_points().size());
608  AssertDimension(support_point_values.size(), nodal_values.size());
609  AssertDimension(this->dofs_per_cell, nodal_values.size());
610 
611  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
612  {
613  AssertDimension(support_point_values[i].size(), 1);
614 
615  nodal_values[i] = support_point_values[i](0);
616  }
617 }
618 
619 
620 
621 template <int dim, int spacedim>
622 FE_SimplexP<dim, spacedim>::FE_SimplexP(const unsigned int degree)
623  : FE_SimplexPoly<dim, spacedim>(degree,
624  get_dpo_vector_fe_p(dim, degree),
625  FiniteElementData<dim>::H1)
626 {}
627 
628 
629 
630 template <int dim, int spacedim>
631 std::unique_ptr<FiniteElement<dim, spacedim>>
633 {
634  return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
635 }
636 
637 
638 
639 template <int dim, int spacedim>
640 std::string
642 {
643  std::ostringstream namebuf;
644  namebuf << "FE_SimplexP<" << dim << ">(" << this->degree << ")";
645 
646  return namebuf.str();
647 }
648 
649 
650 
651 template <int dim, int spacedim>
654  const FiniteElement<dim, spacedim> &fe_other,
655  const unsigned int codim) const
656 {
657  Assert(codim <= dim, ExcImpossibleInDim(dim));
658 
659  // vertex/line/face domination
660  // (if fe_other is derived from FE_SimplexDGP)
661  // ------------------------------------
662  if (codim > 0)
663  if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
664  nullptr)
665  // there are no requirements between continuous and discontinuous
666  // elements
668 
669  // vertex/line/face domination
670  // (if fe_other is not derived from FE_SimplexDGP)
671  // & cell domination
672  // ----------------------------------------
673  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
674  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
675  {
676  if (this->degree < fe_p_other->degree)
678  else if (this->degree == fe_p_other->degree)
680  else
682  }
683  else if (const FE_Q<dim, spacedim> *fe_q_other =
684  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
685  {
686  if (this->degree < fe_q_other->degree)
688  else if (this->degree == fe_q_other->degree)
690  else
692  }
693  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
694  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
695  {
696  if (fe_nothing->is_dominating())
698  else
699  // the FE_Nothing has no degrees of freedom and it is typically used
700  // in a context where we don't require any continuity along the
701  // interface
703  }
704 
705  Assert(false, ExcNotImplemented());
707 }
708 
709 
710 
711 template <int dim, int spacedim>
712 std::vector<std::pair<unsigned int, unsigned int>>
714  const FiniteElement<dim, spacedim> &fe_other) const
715 {
716  AssertDimension(dim, 2);
717 
718  if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
719  {
720  // there should be exactly one single DoF of each FE at a vertex, and
721  // they should have identical value
722  return {{0U, 0U}};
723  }
724  else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
725  {
726  // there should be exactly one single DoF of each FE at a vertex, and
727  // they should have identical value
728  return {{0U, 0U}};
729  }
730  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
731  {
732  // the FE_Nothing has no degrees of freedom, so there are no
733  // equivalencies to be recorded
734  return {};
735  }
736  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
737  {
738  // if the other element has no elements on faces at all,
739  // then it would be impossible to enforce any kind of
740  // continuity even if we knew exactly what kind of element
741  // we have -- simply because the other element declares
742  // that it is discontinuous because it has no DoFs on
743  // its faces. in that case, just state that we have no
744  // constraints to declare
745  return {};
746  }
747  else
748  {
749  Assert(false, ExcNotImplemented());
750  return {};
751  }
752 }
753 
754 
755 
756 template <int dim, int spacedim>
757 std::vector<std::pair<unsigned int, unsigned int>>
759  const FiniteElement<dim, spacedim> &fe_other) const
760 {
761  AssertDimension(dim, 2);
762  Assert(this->degree <= 2, ExcNotImplemented());
763 
764  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
765  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
766  {
767  // dofs are located along lines, so two dofs are identical if they are
768  // located at identical positions.
769  // Therefore, read the points in unit_support_points for the
770  // first coordinate direction. For FE_SimplexP, they are currently
771  // hard-coded and we iterate over points on the first line which begin
772  // after the 3 vertex points in the complete list of unit support points
773 
774  Assert(fe_p_other->degree <= 2, ExcNotImplemented());
775 
776  std::vector<std::pair<unsigned int, unsigned int>> identities;
777 
778  for (unsigned int i = 0; i < this->degree - 1; ++i)
779  for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
780  if (std::fabs(this->unit_support_points[i + 3][0] -
781  fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
782  identities.emplace_back(i, j);
783 
784  return identities;
785  }
786  else if (const FE_Q<dim, spacedim> *fe_q_other =
787  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
788  {
789  // dofs are located along lines, so two dofs are identical if they are
790  // located at identical positions. if we had only equidistant points, we
791  // could simply check for similarity like (i+1)*q == (j+1)*p, but we
792  // might have other support points (e.g. Gauss-Lobatto
793  // points). Therefore, read the points in unit_support_points for the
794  // first coordinate direction. For FE_Q, we take the lexicographic
795  // ordering of the line support points in the first direction (i.e.,
796  // x-direction), which we access between index 1 and p-1 (index 0 and p
797  // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
798  // iterate over points on the first line which begin after the 3 vertex
799  // points in the complete list of unit support points
800 
801  const std::vector<unsigned int> &index_map_inverse_q_other =
802  fe_q_other->get_poly_space_numbering_inverse();
803 
804  std::vector<std::pair<unsigned int, unsigned int>> identities;
805 
806  for (unsigned int i = 0; i < this->degree - 1; ++i)
807  for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
808  if (std::fabs(this->unit_support_points[i + 3][0] -
809  fe_q_other->get_unit_support_points()
810  [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
811  identities.emplace_back(i, j);
812 
813  return identities;
814  }
815  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
816  {
817  // the FE_Nothing has no degrees of freedom, so there are no
818  // equivalencies to be recorded
819  return {};
820  }
821  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
822  {
823  // if the other element has no elements on faces at all,
824  // then it would be impossible to enforce any kind of
825  // continuity even if we knew exactly what kind of element
826  // we have -- simply because the other element declares
827  // that it is discontinuous because it has no DoFs on
828  // its faces. in that case, just state that we have no
829  // constraints to declare
830  return {};
831  }
832  else
833  {
834  Assert(false, ExcNotImplemented());
835  return {};
836  }
837 }
838 
839 
840 
841 template <int dim, int spacedim>
843  : FE_SimplexPoly<dim, spacedim>(degree,
844  get_dpo_vector_fe_dgp(dim, degree),
845  FiniteElementData<dim>::L2)
846 {}
847 
848 
849 
850 template <int dim, int spacedim>
851 std::unique_ptr<FiniteElement<dim, spacedim>>
853 {
854  return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
855 }
856 
857 
858 
859 template <int dim, int spacedim>
860 std::string
862 {
863  std::ostringstream namebuf;
864  namebuf << "FE_SimplexDGP<" << dim << ">(" << this->degree << ")";
865 
866  return namebuf.str();
867 }
868 
869 
870 template <int dim, int spacedim>
873  const FiniteElement<dim, spacedim> &fe_other,
874  const unsigned int codim) const
875 {
876  Assert(codim <= dim, ExcImpossibleInDim(dim));
877 
878  // vertex/line/face domination
879  // ---------------------------
880  if (codim > 0)
881  // this is a discontinuous element, so by definition there will
882  // be no constraints wherever this element comes together with
883  // any other kind of element
885 
886  // cell domination
887  // ---------------
888  if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
889  dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
890  {
891  if (this->degree < fe_dgp_other->degree)
893  else if (this->degree == fe_dgp_other->degree)
895  else
897  }
898  else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
899  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
900  {
901  if (this->degree < fe_dgq_other->degree)
903  else if (this->degree == fe_dgq_other->degree)
905  else
907  }
908  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
909  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
910  {
911  if (fe_nothing->is_dominating())
913  else
914  // the FE_Nothing has no degrees of freedom and it is typically used
915  // in a context where we don't require any continuity along the
916  // interface
918  }
919 
920  Assert(false, ExcNotImplemented());
922 }
923 
924 
925 
926 template <int dim, int spacedim>
927 std::vector<std::pair<unsigned int, unsigned int>>
929  const FiniteElement<dim, spacedim> &fe_other) const
930 {
931  (void)fe_other;
932 
933  return {};
934 }
935 
936 
937 
938 template <int dim, int spacedim>
939 std::vector<std::pair<unsigned int, unsigned int>>
941  const FiniteElement<dim, spacedim> &fe_other) const
942 {
943  (void)fe_other;
944 
945  return {};
946 }
947 
948 // explicit instantiations
949 #include "fe_simplex_p.inst"
950 
size_type m() const
Number * begin_raw()
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2404
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1657
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
FullMatrix< double > interface_constraints
Definition: fe.h:2430
static void project_to_subface(const SubQuadrature &quadrature, const unsigned int face_no, const unsigned int subface_no, std::vector< Point< dim >> &q_points, const RefinementCase< dim - 1 > &ref_case=RefinementCase< dim - 1 >::isotropic_refinement)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Definition: fe_dgq.h:110
constexpr const ReferenceCell Triangle
const unsigned int degree
Definition: fe_base.h:435
STL namespace.
static const char U
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
static ::ExceptionBase & ExcInterpolationNotImplemented()
FE_SimplexPoly(const unsigned int degree, const std::vector< unsigned int > &dpo_vector, const typename FiniteElementData< dim >::Conformity conformity)
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2442
constexpr const ReferenceCell Line
Number * end_raw()
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
Definition: fe_q.h:548
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_SimplexP(const unsigned int degree)
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
constexpr const ReferenceCell Tetrahedron
std::string get_name() const override
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2418
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
T sum(const T &t, const MPI_Comm &mpi_communicator)
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
Definition: fe.cc:1108
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2449
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:399
unsigned int n_vertices() const
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Expression fabs(const Expression &x)
virtual unsigned int face_to_cell_index(const unsigned int face_dof_index, const unsigned int face, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false) const
Definition: fe.cc:568
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source_fe, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
const unsigned int dofs_per_cell
Definition: fe_base.h:419
std::string get_name() const override
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
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:1049
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
unsigned int n_lines() const
Threads::Mutex mutex
Definition: fe_simplex_p.h:106
void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &x_source_fe, const unsigned int subface, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
unsigned int n_dofs_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:398
bool hp_constraints_are_implemented() const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_SimplexDGP(const unsigned int degree)
static ::ExceptionBase & ExcNotImplemented()
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
ReferenceCell reference_cell() const
static ::ExceptionBase & ExcInternalError()