Reference documentation for deal.II version GIT 7b2de2f2f9 2023-09-24 11:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_simplex_p.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #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_p(const unsigned int degree)
85  {
86  std::vector<Point<dim>> unit_points;
87  // If we do dim - 1 we can get here in dim = 0:
88  if (dim == 0)
89  {
90  unit_points.emplace_back();
91  return unit_points;
92  }
93 
94  const auto reference_cell = ReferenceCells::get_simplex<dim>();
95  // Piecewise constants are a special case: use a support point at the
96  // centroid and only the centroid
97  if (degree == 0)
98  {
99  unit_points.emplace_back(reference_cell.template barycenter<dim>());
100  return unit_points;
101  }
102 
103  Assert(degree <= 2, ExcNotImplemented());
104  for (const auto &vertex_no : reference_cell.vertex_indices())
105  unit_points.emplace_back(reference_cell.template vertex<dim>(vertex_no));
106 
107  if (degree == 2)
108  for (const auto &line_no : reference_cell.line_indices())
109  {
110  const auto v0 = reference_cell.template vertex<dim>(
111  reference_cell.line_to_cell_vertices(line_no, 0));
112  const auto v1 = reference_cell.template vertex<dim>(
113  reference_cell.line_to_cell_vertices(line_no, 1));
114  unit_points.emplace_back((v0 + v1) / 2.0);
115  }
116 
117  return unit_points;
118  }
119 
124  template <int dim>
125  std::vector<std::vector<Point<dim - 1>>>
126  unit_face_support_points_fe_p(
127  const unsigned int degree,
128  typename FiniteElementData<dim>::Conformity conformity)
129  {
130  // Discontinuous elements don't have face support points
131  if (conformity == FiniteElementData<dim>::Conformity::L2)
132  return {};
133 
134  // this concept doesn't exist in 1d so just return an empty vector
135  if (dim == 1)
136  return {};
137 
138  std::vector<std::vector<Point<dim - 1>>> unit_face_points;
139 
140  // all faces have the same support points
141  for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
142  {
143  (void)face_n;
144  unit_face_points.emplace_back(
145  unit_support_points_fe_p<dim - 1>(degree));
146  }
147 
148  return unit_face_points;
149  }
150 
156  template <int dim>
158  constraints_fe_p(const unsigned int /*degree*/)
159  {
160  // no constraints in 1d
161  // constraints in 3d not implemented yet
162  return FullMatrix<double>();
163  }
164 
165  template <>
167  constraints_fe_p<2>(const unsigned int degree)
168  {
169  const unsigned int dim = 2;
170 
171  Assert(degree <= 2, ExcNotImplemented());
172 
173  // the following implements the 2d case
174  // (the 3d case is not implemented yet)
175  //
176  // consult FE_Q_Base::Implementation::initialize_constraints()
177  // for more information
178 
179  std::vector<Point<dim - 1>> constraint_points;
180  // midpoint
181  constraint_points.emplace_back(0.5);
182  if (degree == 2)
183  {
184  // midpoint on subface 0
185  constraint_points.emplace_back(0.25);
186  // midpoint on subface 1
187  constraint_points.emplace_back(0.75);
188  }
189 
190  // Now construct relation between destination (child) and source (mother)
191  // dofs.
192 
193  const unsigned int n_dofs_constrained = constraint_points.size();
194  unsigned int n_dofs_per_face = degree + 1;
195  FullMatrix<double> interface_constraints(n_dofs_constrained,
196  n_dofs_per_face);
197 
198  const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
199 
200  for (unsigned int i = 0; i < n_dofs_constrained; ++i)
201  for (unsigned int j = 0; j < n_dofs_per_face; ++j)
202  {
203  interface_constraints(i, j) =
204  poly.compute_value(j, constraint_points[i]);
205 
206  // if the value is small up to round-off, then simply set it to zero
207  // to avoid unwanted fill-in of the constraint matrices (which would
208  // then increase the number of other DoFs a constrained DoF would
209  // couple to)
210  if (std::fabs(interface_constraints(i, j)) < 1e-13)
211  interface_constraints(i, j) = 0;
212  }
213  return interface_constraints;
214  }
215 
216 
217 
222  std::vector<unsigned int>
223  get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
224  {
225  // First treat the case of piecewise constant elements:
226  if (degree == 0)
227  {
228  std::vector<unsigned int> dpo(dim + 1, 0U);
229  dpo[dim] = 1;
230  return dpo;
231  }
232  else
233  {
234  // This element has the same degrees of freedom as the continuous one,
235  // but they are all counted for the interior of the cell because
236  // it is continuous. Rather than hard-code how many DoFs the element
237  // has, we just get the numbers from the continuous case and add them
238  // up
239  const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
240 
241  switch (dim)
242  {
243  case 1:
244  return {0U,
245  ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
246  continuous_dpo[dim]};
247 
248  case 2:
249  return {0U,
250  0U,
252  continuous_dpo[0] +
253  ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
254  continuous_dpo[dim]};
255 
256  case 3:
257  return {
258  0U,
259  0U,
260  0U,
261  ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
262  ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
263  ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
264  continuous_dpo[dim]};
265  }
266 
267  Assert(false, ExcNotImplemented());
268  return {};
269  }
270  }
271 } // namespace
272 
273 
274 
275 template <int dim, int spacedim>
277  const BarycentricPolynomials<dim> polynomials,
278  const FiniteElementData<dim> &fe_data,
279  const std::vector<Point<dim>> &unit_support_points,
280  const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
281  const FullMatrix<double> &interface_constraints)
282  : ::FE_Poly<dim, spacedim>(
283  polynomials,
284  fe_data,
285  std::vector<bool>(fe_data.dofs_per_cell),
286  std::vector<ComponentMask>(fe_data.dofs_per_cell,
287  ComponentMask(std::vector<bool>(1, true))))
288 {
291  this->interface_constraints = interface_constraints;
292 }
293 
294 
295 
296 template <int dim, int spacedim>
297 std::pair<Table<2, bool>, std::vector<unsigned int>>
299 {
300  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
301  constant_modes.fill(true);
302  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
303  constant_modes, std::vector<unsigned int>(1, 0));
304 }
305 
306 
307 
308 template <int dim, int spacedim>
309 const FullMatrix<double> &
311  const unsigned int child,
312  const RefinementCase<dim> &refinement_case) const
313 {
316  AssertDimension(dim, spacedim);
317 
318  // initialization upon first request
319  if (this->prolongation[refinement_case - 1][child].n() == 0)
320  {
321  std::lock_guard<std::mutex> lock(this->mutex);
322 
323  // if matrix got updated while waiting for the lock
324  if (this->prolongation[refinement_case - 1][child].n() ==
325  this->n_dofs_per_cell())
326  return this->prolongation[refinement_case - 1][child];
327 
328  // now do the work. need to get a non-const version of data in order to
329  // be able to modify them inside a const function
330  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
331 
332  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
334  isotropic_matrices.back().resize(
336  FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
337 
338  FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
339 
340  this_nonconst.prolongation[refinement_case - 1].swap(
341  isotropic_matrices.back());
342  }
343 
344  // finally return the matrix
345  return this->prolongation[refinement_case - 1][child];
346 }
347 
348 
349 
350 template <int dim, int spacedim>
351 const FullMatrix<double> &
353  const unsigned int child,
354  const RefinementCase<dim> &refinement_case) const
355 {
358  AssertDimension(dim, spacedim);
359 
360  // initialization upon first request
361  if (this->restriction[refinement_case - 1][child].n() == 0)
362  {
363  std::lock_guard<std::mutex> lock(this->mutex);
364 
365  // if matrix got updated while waiting for the lock
366  if (this->restriction[refinement_case - 1][child].n() ==
367  this->n_dofs_per_cell())
368  return this->restriction[refinement_case - 1][child];
369 
370  // now do the work. need to get a non-const version of data in order to
371  // be able to modify them inside a const function
372  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
373 
374  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
376  isotropic_matrices.back().resize(
378  FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
379 
380  FETools::compute_projection_matrices(*this, isotropic_matrices, true);
381 
382  this_nonconst.restriction[refinement_case - 1].swap(
383  isotropic_matrices.back());
384  }
385 
386  // finally return the matrix
387  return this->restriction[refinement_case - 1][child];
388 }
389 
390 
391 
392 template <int dim, int spacedim>
393 void
395  const FiniteElement<dim, spacedim> &source_fe,
396  FullMatrix<double> &interpolation_matrix,
397  const unsigned int face_no) const
398 {
399  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
400  ExcDimensionMismatch(interpolation_matrix.m(),
401  source_fe.n_dofs_per_face(face_no)));
402 
403  // see if source is a P or Q element
404  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
405  nullptr) ||
406  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
407  {
408  const Quadrature<dim - 1> quad_face_support(
409  source_fe.get_unit_face_support_points(face_no));
410 
411  const double eps = 2e-13 * this->degree * (dim - 1);
412 
413  std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
415  quad_face_support,
416  face_no,
417  face_quadrature_points);
418 
419  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
420  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
421  {
422  double matrix_entry =
423  this->shape_value(this->face_to_cell_index(j, 0),
424  face_quadrature_points[i]);
425 
426  // Correct the interpolated value. I.e. if it is close to 1 or
427  // 0, make it exactly 1 or 0. Unfortunately, this is required to
428  // avoid problems with higher order elements.
429  if (std::fabs(matrix_entry - 1.0) < eps)
430  matrix_entry = 1.0;
431  if (std::fabs(matrix_entry) < eps)
432  matrix_entry = 0.0;
433 
434  interpolation_matrix(i, j) = matrix_entry;
435  }
436 
437 #ifdef DEBUG
438  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
439  {
440  double sum = 0.;
441 
442  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
443  sum += interpolation_matrix(j, i);
444 
446  }
447 #endif
448  }
449  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
450  {
451  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
452  }
453  else
454  AssertThrow(
455  false,
456  (typename FiniteElement<dim,
457  spacedim>::ExcInterpolationNotImplemented()));
458 }
459 
460 
461 
462 template <int dim, int spacedim>
463 void
465  const FiniteElement<dim, spacedim> &source_fe,
466  const unsigned int subface,
467  FullMatrix<double> &interpolation_matrix,
468  const unsigned int face_no) const
469 {
470  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
471  ExcDimensionMismatch(interpolation_matrix.m(),
472  source_fe.n_dofs_per_face(face_no)));
473 
474  // see if source is a P or Q element
475  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
476  nullptr) ||
477  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
478  {
479  const Quadrature<dim - 1> quad_face_support(
480  source_fe.get_unit_face_support_points(face_no));
481 
482  const double eps = 2e-13 * this->degree * (dim - 1);
483 
484  std::vector<Point<dim>> subface_quadrature_points(
485  quad_face_support.size());
487  quad_face_support,
488  face_no,
489  subface,
490  subface_quadrature_points);
491 
492  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
493  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
494  {
495  double matrix_entry =
496  this->shape_value(this->face_to_cell_index(j, 0),
497  subface_quadrature_points[i]);
498 
499  // Correct the interpolated value. I.e. if it is close to 1 or
500  // 0, make it exactly 1 or 0. Unfortunately, this is required to
501  // avoid problems with higher order elements.
502  if (std::fabs(matrix_entry - 1.0) < eps)
503  matrix_entry = 1.0;
504  if (std::fabs(matrix_entry) < eps)
505  matrix_entry = 0.0;
506 
507  interpolation_matrix(i, j) = matrix_entry;
508  }
509 
510 #ifdef DEBUG
511  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
512  {
513  double sum = 0.;
514 
515  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
516  sum += interpolation_matrix(j, i);
517 
519  }
520 #endif
521  }
522  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
523  {
524  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
525  }
526  else
527  AssertThrow(
528  false,
529  (typename FiniteElement<dim,
530  spacedim>::ExcInterpolationNotImplemented()));
531 }
532 
533 
534 
535 template <int dim, int spacedim>
536 bool
538 {
539  return true;
540 }
541 
542 
543 
544 template <int dim, int spacedim>
545 void
548  const std::vector<Vector<double>> &support_point_values,
549  std::vector<double> &nodal_values) const
550 {
551  AssertDimension(support_point_values.size(),
552  this->get_unit_support_points().size());
553  AssertDimension(support_point_values.size(), nodal_values.size());
554  AssertDimension(this->dofs_per_cell, nodal_values.size());
555 
556  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
557  {
558  AssertDimension(support_point_values[i].size(), 1);
559 
560  nodal_values[i] = support_point_values[i](0);
561  }
562 }
563 
564 
565 
566 template <int dim, int spacedim>
567 FE_SimplexP<dim, spacedim>::FE_SimplexP(const unsigned int degree)
568  : FE_SimplexPoly<dim, spacedim>(
569  BarycentricPolynomials<dim>::get_fe_p_basis(degree),
570  FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
571  ReferenceCells::get_simplex<dim>(),
572  1,
573  degree,
574  FiniteElementData<dim>::H1),
575  unit_support_points_fe_p<dim>(degree),
576  unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
577  constraints_fe_p<dim>(degree))
578 {}
579 
580 
581 
582 template <int dim, int spacedim>
583 std::unique_ptr<FiniteElement<dim, spacedim>>
585 {
586  return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
587 }
588 
589 
590 
591 template <int dim, int spacedim>
592 std::string
594 {
595  std::ostringstream namebuf;
596  namebuf << "FE_SimplexP<" << dim << ">(" << this->degree << ")";
597 
598  return namebuf.str();
599 }
600 
601 
602 
603 template <int dim, int spacedim>
606  const FiniteElement<dim, spacedim> &fe_other,
607  const unsigned int codim) const
608 {
609  Assert(codim <= dim, ExcImpossibleInDim(dim));
610 
611  // vertex/line/face domination
612  // (if fe_other is derived from FE_SimplexDGP)
613  // ------------------------------------
614  if (codim > 0)
615  if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
616  nullptr)
617  // there are no requirements between continuous and discontinuous
618  // elements
620 
621  // vertex/line/face domination
622  // (if fe_other is not derived from FE_SimplexDGP)
623  // & cell domination
624  // ----------------------------------------
625  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
626  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
627  {
628  if (this->degree < fe_p_other->degree)
630  else if (this->degree == fe_p_other->degree)
632  else
634  }
635  else if (const FE_Q<dim, spacedim> *fe_q_other =
636  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
637  {
638  if (this->degree < fe_q_other->degree)
640  else if (this->degree == fe_q_other->degree)
642  else
644  }
645  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
646  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
647  {
648  if (fe_nothing->is_dominating())
650  else
651  // the FE_Nothing has no degrees of freedom and it is typically used
652  // in a context where we don't require any continuity along the
653  // interface
655  }
656 
657  Assert(false, ExcNotImplemented());
659 }
660 
661 
662 
663 template <int dim, int spacedim>
664 std::vector<std::pair<unsigned int, unsigned int>>
666  const FiniteElement<dim, spacedim> &fe_other) const
667 {
668  AssertDimension(dim, 2);
669 
670  if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
671  {
672  // there should be exactly one single DoF of each FE at a vertex, and
673  // they should have identical value
674  return {{0U, 0U}};
675  }
676  else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
677  {
678  // there should be exactly one single DoF of each FE at a vertex, and
679  // they should have identical value
680  return {{0U, 0U}};
681  }
682  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
683  {
684  // the FE_Nothing has no degrees of freedom, so there are no
685  // equivalencies to be recorded
686  return {};
687  }
688  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
689  {
690  // if the other element has no DoFs on faces at all,
691  // then it would be impossible to enforce any kind of
692  // continuity even if we knew exactly what kind of element
693  // we have -- simply because the other element declares
694  // that it is discontinuous because it has no DoFs on
695  // its faces. in that case, just state that we have no
696  // constraints to declare
697  return {};
698  }
699  else
700  {
701  Assert(false, ExcNotImplemented());
702  return {};
703  }
704 }
705 
706 
707 
708 template <int dim, int spacedim>
709 std::vector<std::pair<unsigned int, unsigned int>>
711  const FiniteElement<dim, spacedim> &fe_other) const
712 {
713  AssertDimension(dim, 2);
714 
715  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
716  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
717  {
718  // dofs are located along lines, so two dofs are identical if they are
719  // located at identical positions.
720  // Therefore, read the points in unit_support_points for the
721  // first coordinate direction. For FE_SimplexP, they are currently
722  // hard-coded and we iterate over points on the first line which begin
723  // after the 3 vertex points in the complete list of unit support points
724 
725  std::vector<std::pair<unsigned int, unsigned int>> identities;
726 
727  for (unsigned int i = 0; i < this->degree - 1; ++i)
728  for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
729  if (std::fabs(this->unit_support_points[i + 3][0] -
730  fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
731  identities.emplace_back(i, j);
732  else
733  {
734  // If nodes are not located in the same place, we have to
735  // interpolate. This is then not handled through the
736  // current function, but via interpolation matrices that
737  // result in constraints, rather than identities. Since
738  // that happens in a different function, there is nothing
739  // for us to do here.
740  }
741 
742  return identities;
743  }
744  else if (const FE_Q<dim, spacedim> *fe_q_other =
745  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
746  {
747  // dofs are located along lines, so two dofs are identical if they are
748  // located at identical positions. if we had only equidistant points, we
749  // could simply check for similarity like (i+1)*q == (j+1)*p, but we
750  // might have other support points (e.g. Gauss-Lobatto
751  // points). Therefore, read the points in unit_support_points for the
752  // first coordinate direction. For FE_Q, we take the lexicographic
753  // ordering of the line support points in the first direction (i.e.,
754  // x-direction), which we access between index 1 and p-1 (index 0 and p
755  // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
756  // iterate over points on the first line which begin after the 3 vertex
757  // points in the complete list of unit support points
758 
759  const std::vector<unsigned int> &index_map_inverse_q_other =
760  fe_q_other->get_poly_space_numbering_inverse();
761 
762  std::vector<std::pair<unsigned int, unsigned int>> identities;
763 
764  for (unsigned int i = 0; i < this->degree - 1; ++i)
765  for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
766  if (std::fabs(this->unit_support_points[i + 3][0] -
767  fe_q_other->get_unit_support_points()
768  [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
769  identities.emplace_back(i, j);
770  else
771  {
772  // If nodes are not located in the same place, we have to
773  // interpolate. This will then also
774  // capture the case where the FE_Q has a different polynomial
775  // degree than the current element. In either case, the resulting
776  // constraints are computed elsewhere, rather than via the
777  // identities this function returns: Since
778  // that happens in a different function, there is nothing
779  // for us to do here.
780  }
781 
782  return identities;
783  }
784  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
785  {
786  // The FE_Nothing has no degrees of freedom, so there are no
787  // equivalencies to be recorded. (If the FE_Nothing is dominating,
788  // then this will also leads to constraints, but we are not concerned
789  // with this here.)
790  return {};
791  }
792  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
793  {
794  // if the other element has no elements on faces at all,
795  // then it would be impossible to enforce any kind of
796  // continuity even if we knew exactly what kind of element
797  // we have -- simply because the other element declares
798  // that it is discontinuous because it has no DoFs on
799  // its faces. in that case, just state that we have no
800  // constraints to declare
801  return {};
802  }
803  else
804  {
805  Assert(false, ExcNotImplemented());
806  return {};
807  }
808 }
809 
810 
811 
812 template <int dim, int spacedim>
814  : FE_SimplexPoly<dim, spacedim>(
815  BarycentricPolynomials<dim>::get_fe_p_basis(degree),
816  FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
817  ReferenceCells::get_simplex<dim>(),
818  1,
819  degree,
820  FiniteElementData<dim>::L2),
821  unit_support_points_fe_p<dim>(degree),
822  unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::L2),
823  constraints_fe_p<dim>(degree))
824 {}
825 
826 
827 
828 template <int dim, int spacedim>
829 std::unique_ptr<FiniteElement<dim, spacedim>>
831 {
832  return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
833 }
834 
835 
836 
837 template <int dim, int spacedim>
838 std::string
840 {
841  std::ostringstream namebuf;
842  namebuf << "FE_SimplexDGP<" << dim << ">(" << this->degree << ")";
843 
844  return namebuf.str();
845 }
846 
847 
848 template <int dim, int spacedim>
851  const FiniteElement<dim, spacedim> &fe_other,
852  const unsigned int codim) const
853 {
854  Assert(codim <= dim, ExcImpossibleInDim(dim));
855 
856  // vertex/line/face domination
857  // ---------------------------
858  if (codim > 0)
859  // this is a discontinuous element, so by definition there will
860  // be no constraints wherever this element comes together with
861  // any other kind of element
863 
864  // cell domination
865  // ---------------
866  if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
867  dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
868  {
869  if (this->degree < fe_dgp_other->degree)
871  else if (this->degree == fe_dgp_other->degree)
873  else
875  }
876  else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
877  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
878  {
879  if (this->degree < fe_dgq_other->degree)
881  else if (this->degree == fe_dgq_other->degree)
883  else
885  }
886  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
887  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
888  {
889  if (fe_nothing->is_dominating())
891  else
892  // the FE_Nothing has no degrees of freedom and it is typically used
893  // in a context where we don't require any continuity along the
894  // interface
896  }
897 
898  Assert(false, ExcNotImplemented());
900 }
901 
902 
903 
904 template <int dim, int spacedim>
905 std::vector<std::pair<unsigned int, unsigned int>>
907  const FiniteElement<dim, spacedim> &fe_other) const
908 {
909  (void)fe_other;
910 
911  return {};
912 }
913 
914 
915 
916 template <int dim, int spacedim>
917 std::vector<std::pair<unsigned int, unsigned int>>
919  const FiniteElement<dim, spacedim> &fe_other) const
920 {
921  (void)fe_other;
922 
923  return {};
924 }
925 
926 // explicit instantiations
927 #include "fe_simplex_p.inst"
928 
static BarycentricPolynomials< dim > get_fe_p_basis(const unsigned int degree)
Definition: fe_dgq.h:113
Definition: fe_q.h:551
FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim) const override
std::string get_name() const override
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_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
FE_SimplexDGP(const unsigned int degree)
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_SimplexP(const unsigned int degree)
std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
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
std::string get_name() const override
std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
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
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
FE_SimplexPoly(const BarycentricPolynomials< dim > polynomials, const FiniteElementData< dim > &fe_data, const std::vector< Point< dim >> &unit_support_points, const std::vector< std::vector< Point< dim - 1 >>> unit_face_support_points, const FullMatrix< double > &interface_constraints)
bool hp_constraints_are_implemented() const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) 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
void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source_fe, FullMatrix< double > &interpolation_matrix, const unsigned int face_no) const override
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
std::vector< std::vector< Point< dim - 1 > > > unit_face_support_points
Definition: fe.h:2447
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2440
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
FullMatrix< double > interface_constraints
Definition: fe.h:2428
size_type m() const
Definition: point.h:112
static void project_to_subface(const ReferenceCell &reference_cell, 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)
static void project_to_face(const ReferenceCell &reference_cell, const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
unsigned int n_vertices() const
unsigned int n_faces() const
unsigned int n_lines() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
const unsigned int v0
Definition: grid_tools.cc:1063
const unsigned int v1
Definition: grid_tools.cc:1063
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression fabs(const Expression &x)
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)
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char U
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
constexpr const ReferenceCell Tetrahedron
constexpr const ReferenceCell Triangle
constexpr const ReferenceCell Line
constexpr const ReferenceCell & get_simplex()
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)