Reference documentation for deal.II version GIT f9a1fd2718 2022-08-19 02:45: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 - 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_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  // Piecewise constants are a special case: use a support point at the
95  // centroid and only the centroid
96  if (degree == 0)
97  {
98  Point<dim> centroid;
99  std::fill(centroid.begin_raw(),
100  centroid.end_raw(),
101  1.0 / double(dim + 1));
102  unit_points.emplace_back(centroid);
103  return unit_points;
104  }
105 
106  if (dim == 1)
107  {
108  // We don't really have dim = 1 support for simplex elements yet, but
109  // its convenient for populating the face array
110  Assert(degree <= 2, ExcNotImplemented());
111  if (degree >= 1)
112  {
113  unit_points.emplace_back(0.0);
114  unit_points.emplace_back(1.0);
115 
116  if (degree == 2)
117  unit_points.emplace_back(0.5);
118  }
119  }
120  else if (dim == 2)
121  {
122  Assert(degree <= 2, ExcNotImplemented());
123  if (degree >= 1)
124  {
125  unit_points.emplace_back(0.0, 0.0);
126  unit_points.emplace_back(1.0, 0.0);
127  unit_points.emplace_back(0.0, 1.0);
128 
129  if (degree == 2)
130  {
131  unit_points.emplace_back(0.5, 0.0);
132  unit_points.emplace_back(0.5, 0.5);
133  unit_points.emplace_back(0.0, 0.5);
134  }
135  }
136  }
137  else if (dim == 3)
138  {
139  Assert(degree <= 2, ExcNotImplemented());
140  if (degree >= 1)
141  {
142  unit_points.emplace_back(0.0, 0.0, 0.0);
143  unit_points.emplace_back(1.0, 0.0, 0.0);
144  unit_points.emplace_back(0.0, 1.0, 0.0);
145  unit_points.emplace_back(0.0, 0.0, 1.0);
146 
147  if (degree == 2)
148  {
149  unit_points.emplace_back(0.5, 0.0, 0.0);
150  unit_points.emplace_back(0.5, 0.5, 0.0);
151  unit_points.emplace_back(0.0, 0.5, 0.0);
152  unit_points.emplace_back(0.0, 0.0, 0.5);
153  unit_points.emplace_back(0.5, 0.0, 0.5);
154  unit_points.emplace_back(0.0, 0.5, 0.5);
155  }
156  }
157  }
158  else
159  {
160  Assert(false, ExcNotImplemented());
161  }
162 
163  return unit_points;
164  }
165 
170  template <int dim>
171  std::vector<std::vector<Point<dim - 1>>>
172  unit_face_support_points_fe_p(
173  const unsigned int degree,
174  typename FiniteElementData<dim>::Conformity conformity)
175  {
176  // Discontinuous elements don't have face support points
177  if (conformity == FiniteElementData<dim>::Conformity::L2)
178  return {};
179 
180  // this concept doesn't exist in 1D so just return an empty vector
181  if (dim == 1)
182  return {};
183 
184  std::vector<std::vector<Point<dim - 1>>> unit_face_points;
185 
186  // all faces have the same support points
187  for (auto face_n : ReferenceCells::get_simplex<dim>().face_indices())
188  {
189  (void)face_n;
190  unit_face_points.emplace_back(
191  unit_support_points_fe_p<dim - 1>(degree));
192  }
193 
194  return unit_face_points;
195  }
196 
202  template <int dim>
204  constraints_fe_p(const unsigned int /*degree*/)
205  {
206  // no constraints in 1d
207  // constraints in 3d not implemented yet
208  return FullMatrix<double>();
209  }
210 
211  template <>
213  constraints_fe_p<2>(const unsigned int degree)
214  {
215  const unsigned int dim = 2;
216 
217  Assert(degree <= 2, ExcNotImplemented());
218 
219  // the following implements the 2d case
220  // (the 3d case is not implemented yet)
221  //
222  // consult FE_Q_Base::Implementation::initialize_constraints()
223  // for more information
224 
225  std::vector<Point<dim - 1>> constraint_points;
226  // midpoint
227  constraint_points.emplace_back(0.5);
228  if (degree == 2)
229  {
230  // midpoint on subface 0
231  constraint_points.emplace_back(0.25);
232  // midpoint on subface 1
233  constraint_points.emplace_back(0.75);
234  }
235 
236  // Now construct relation between destination (child) and source (mother)
237  // dofs.
238 
239  const unsigned int n_dofs_constrained = constraint_points.size();
240  unsigned int n_dofs_per_face = degree + 1;
241  FullMatrix<double> interface_constraints(n_dofs_constrained,
242  n_dofs_per_face);
243 
244  const auto poly = BarycentricPolynomials<dim - 1>::get_fe_p_basis(degree);
245 
246  for (unsigned int i = 0; i < n_dofs_constrained; ++i)
247  for (unsigned int j = 0; j < n_dofs_per_face; ++j)
248  {
249  interface_constraints(i, j) =
250  poly.compute_value(j, constraint_points[i]);
251 
252  // if the value is small up to round-off, then simply set it to zero
253  // to avoid unwanted fill-in of the constraint matrices (which would
254  // then increase the number of other DoFs a constrained DoF would
255  // couple to)
256  if (std::fabs(interface_constraints(i, j)) < 1e-13)
257  interface_constraints(i, j) = 0;
258  }
259  return interface_constraints;
260  }
261 
262 
263 
268  std::vector<unsigned int>
269  get_dpo_vector_fe_dgp(const unsigned int dim, const unsigned int degree)
270  {
271  // First treat the case of piecewise constant elements:
272  if (degree == 0)
273  {
274  std::vector<unsigned int> dpo(dim + 1, 0U);
275  dpo[dim] = 1;
276  return dpo;
277  }
278  else
279  {
280  // This element has the same degrees of freedom as the continuous one,
281  // but they are all counted for the interior of the cell because
282  // it is continuous. Rather than hard-code how many DoFs the element
283  // has, we just get the numbers from the continuous case and add them
284  // up
285  const auto continuous_dpo = get_dpo_vector_fe_p(dim, degree);
286 
287  switch (dim)
288  {
289  case 1:
290  return {0U,
291  ReferenceCells::Line.n_vertices() * continuous_dpo[0] +
292  continuous_dpo[dim]};
293 
294  case 2:
295  return {0U,
296  0U,
298  continuous_dpo[0] +
299  ReferenceCells::Triangle.n_lines() * continuous_dpo[1] +
300  continuous_dpo[dim]};
301 
302  case 3:
303  return {
304  0U,
305  0U,
306  0U,
307  ReferenceCells::Tetrahedron.n_vertices() * continuous_dpo[0] +
308  ReferenceCells::Tetrahedron.n_lines() * continuous_dpo[1] +
309  ReferenceCells::Tetrahedron.n_faces() * continuous_dpo[2] +
310  continuous_dpo[dim]};
311  }
312 
313  Assert(false, ExcNotImplemented());
314  return {};
315  }
316  }
317 } // namespace
318 
319 
320 
321 template <int dim, int spacedim>
323  const BarycentricPolynomials<dim> polynomials,
324  const FiniteElementData<dim> & fe_data,
325  const std::vector<Point<dim>> & unit_support_points,
326  const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
327  const FullMatrix<double> & interface_constraints)
328  : ::FE_Poly<dim, spacedim>(
329  polynomials,
330  fe_data,
331  std::vector<bool>(fe_data.dofs_per_cell),
332  std::vector<ComponentMask>(fe_data.dofs_per_cell,
333  std::vector<bool>(1, true)))
334 {
337  this->interface_constraints = interface_constraints;
338 }
339 
340 
341 
342 template <int dim, int spacedim>
343 std::pair<Table<2, bool>, std::vector<unsigned int>>
345 {
346  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
347  constant_modes.fill(true);
348  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
349  constant_modes, std::vector<unsigned int>(1, 0));
350 }
351 
352 
353 
354 template <int dim, int spacedim>
355 const FullMatrix<double> &
357  const unsigned int child,
358  const RefinementCase<dim> &refinement_case) const
359 {
362  AssertDimension(dim, spacedim);
363 
364  // initialization upon first request
365  if (this->prolongation[refinement_case - 1][child].n() == 0)
366  {
367  std::lock_guard<std::mutex> lock(this->mutex);
368 
369  // if matrix got updated while waiting for the lock
370  if (this->prolongation[refinement_case - 1][child].n() ==
371  this->n_dofs_per_cell())
372  return this->prolongation[refinement_case - 1][child];
373 
374  // now do the work. need to get a non-const version of data in order to
375  // be able to modify them inside a const function
376  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
377 
378  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
380  isotropic_matrices.back().resize(
382  FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
383 
384  FETools::compute_embedding_matrices(*this, isotropic_matrices, true);
385 
386  this_nonconst.prolongation[refinement_case - 1].swap(
387  isotropic_matrices.back());
388  }
389 
390  // finally return the matrix
391  return this->prolongation[refinement_case - 1][child];
392 }
393 
394 
395 
396 template <int dim, int spacedim>
397 const FullMatrix<double> &
399  const unsigned int child,
400  const RefinementCase<dim> &refinement_case) const
401 {
404  AssertDimension(dim, spacedim);
405 
406  // initialization upon first request
407  if (this->restriction[refinement_case - 1][child].n() == 0)
408  {
409  std::lock_guard<std::mutex> lock(this->mutex);
410 
411  // if matrix got updated while waiting for the lock
412  if (this->restriction[refinement_case - 1][child].n() ==
413  this->n_dofs_per_cell())
414  return this->restriction[refinement_case - 1][child];
415 
416  // now do the work. need to get a non-const version of data in order to
417  // be able to modify them inside a const function
418  auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
419 
420  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
422  isotropic_matrices.back().resize(
424  FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
425 
426  FETools::compute_projection_matrices(*this, isotropic_matrices, true);
427 
428  this_nonconst.restriction[refinement_case - 1].swap(
429  isotropic_matrices.back());
430  }
431 
432  // finally return the matrix
433  return this->restriction[refinement_case - 1][child];
434 }
435 
436 
437 
438 template <int dim, int spacedim>
439 void
441  const FiniteElement<dim, spacedim> &source_fe,
442  FullMatrix<double> & interpolation_matrix,
443  const unsigned int face_no) const
444 {
445  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
446  ExcDimensionMismatch(interpolation_matrix.m(),
447  source_fe.n_dofs_per_face(face_no)));
448 
449  // see if source is a P or Q element
450  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
451  nullptr) ||
452  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
453  {
454  const Quadrature<dim - 1> quad_face_support(
455  source_fe.get_unit_face_support_points(face_no));
456 
457  const double eps = 2e-13 * this->degree * (dim - 1);
458 
459  std::vector<Point<dim>> face_quadrature_points(quad_face_support.size());
461  quad_face_support,
462  face_no,
463  face_quadrature_points);
464 
465  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
466  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
467  {
468  double matrix_entry =
469  this->shape_value(this->face_to_cell_index(j, 0),
470  face_quadrature_points[i]);
471 
472  // Correct the interpolated value. I.e. if it is close to 1 or
473  // 0, make it exactly 1 or 0. Unfortunately, this is required to
474  // avoid problems with higher order elements.
475  if (std::fabs(matrix_entry - 1.0) < eps)
476  matrix_entry = 1.0;
477  if (std::fabs(matrix_entry) < eps)
478  matrix_entry = 0.0;
479 
480  interpolation_matrix(i, j) = matrix_entry;
481  }
482 
483 #ifdef DEBUG
484  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
485  {
486  double sum = 0.;
487 
488  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
489  sum += interpolation_matrix(j, i);
490 
492  }
493 #endif
494  }
495  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
496  {
497  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
498  }
499  else
500  AssertThrow(
501  false,
502  (typename FiniteElement<dim,
503  spacedim>::ExcInterpolationNotImplemented()));
504 }
505 
506 
507 
508 template <int dim, int spacedim>
509 void
511  const FiniteElement<dim, spacedim> &source_fe,
512  const unsigned int subface,
513  FullMatrix<double> & interpolation_matrix,
514  const unsigned int face_no) const
515 {
516  Assert(interpolation_matrix.m() == source_fe.n_dofs_per_face(face_no),
517  ExcDimensionMismatch(interpolation_matrix.m(),
518  source_fe.n_dofs_per_face(face_no)));
519 
520  // see if source is a P or Q element
521  if ((dynamic_cast<const FE_SimplexPoly<dim, spacedim> *>(&source_fe) !=
522  nullptr) ||
523  (dynamic_cast<const FE_Q_Base<dim, spacedim> *>(&source_fe) != nullptr))
524  {
525  const Quadrature<dim - 1> quad_face_support(
526  source_fe.get_unit_face_support_points(face_no));
527 
528  const double eps = 2e-13 * this->degree * (dim - 1);
529 
530  std::vector<Point<dim>> subface_quadrature_points(
531  quad_face_support.size());
533  quad_face_support,
534  face_no,
535  subface,
536  subface_quadrature_points);
537 
538  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
539  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
540  {
541  double matrix_entry =
542  this->shape_value(this->face_to_cell_index(j, 0),
543  subface_quadrature_points[i]);
544 
545  // Correct the interpolated value. I.e. if it is close to 1 or
546  // 0, make it exactly 1 or 0. Unfortunately, this is required to
547  // avoid problems with higher order elements.
548  if (std::fabs(matrix_entry - 1.0) < eps)
549  matrix_entry = 1.0;
550  if (std::fabs(matrix_entry) < eps)
551  matrix_entry = 0.0;
552 
553  interpolation_matrix(i, j) = matrix_entry;
554  }
555 
556 #ifdef DEBUG
557  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
558  {
559  double sum = 0.;
560 
561  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
562  sum += interpolation_matrix(j, i);
563 
565  }
566 #endif
567  }
568  else if (dynamic_cast<const FE_Nothing<dim> *>(&source_fe) != nullptr)
569  {
570  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
571  }
572  else
573  AssertThrow(
574  false,
575  (typename FiniteElement<dim,
576  spacedim>::ExcInterpolationNotImplemented()));
577 }
578 
579 
580 
581 template <int dim, int spacedim>
582 bool
584 {
585  return true;
586 }
587 
588 
589 
590 template <int dim, int spacedim>
591 void
594  const std::vector<Vector<double>> &support_point_values,
595  std::vector<double> & nodal_values) const
596 {
597  AssertDimension(support_point_values.size(),
598  this->get_unit_support_points().size());
599  AssertDimension(support_point_values.size(), nodal_values.size());
600  AssertDimension(this->dofs_per_cell, nodal_values.size());
601 
602  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
603  {
604  AssertDimension(support_point_values[i].size(), 1);
605 
606  nodal_values[i] = support_point_values[i](0);
607  }
608 }
609 
610 
611 
612 template <int dim, int spacedim>
613 FE_SimplexP<dim, spacedim>::FE_SimplexP(const unsigned int degree)
614  : FE_SimplexPoly<dim, spacedim>(
615  BarycentricPolynomials<dim>::get_fe_p_basis(degree),
616  FiniteElementData<dim>(get_dpo_vector_fe_p(dim, degree),
617  ReferenceCells::get_simplex<dim>(),
618  1,
619  degree,
620  FiniteElementData<dim>::H1),
621  unit_support_points_fe_p<dim>(degree),
622  unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
623  constraints_fe_p<dim>(degree))
624 {}
625 
626 
627 
628 template <int dim, int spacedim>
629 std::unique_ptr<FiniteElement<dim, spacedim>>
631 {
632  return std::make_unique<FE_SimplexP<dim, spacedim>>(*this);
633 }
634 
635 
636 
637 template <int dim, int spacedim>
638 std::string
640 {
641  std::ostringstream namebuf;
642  namebuf << "FE_SimplexP<" << dim << ">(" << this->degree << ")";
643 
644  return namebuf.str();
645 }
646 
647 
648 
649 template <int dim, int spacedim>
652  const FiniteElement<dim, spacedim> &fe_other,
653  const unsigned int codim) const
654 {
655  Assert(codim <= dim, ExcImpossibleInDim(dim));
656 
657  // vertex/line/face domination
658  // (if fe_other is derived from FE_SimplexDGP)
659  // ------------------------------------
660  if (codim > 0)
661  if (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other) !=
662  nullptr)
663  // there are no requirements between continuous and discontinuous
664  // elements
666 
667  // vertex/line/face domination
668  // (if fe_other is not derived from FE_SimplexDGP)
669  // & cell domination
670  // ----------------------------------------
671  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
672  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
673  {
674  if (this->degree < fe_p_other->degree)
676  else if (this->degree == fe_p_other->degree)
678  else
680  }
681  else if (const FE_Q<dim, spacedim> *fe_q_other =
682  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
683  {
684  if (this->degree < fe_q_other->degree)
686  else if (this->degree == fe_q_other->degree)
688  else
690  }
691  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
692  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
693  {
694  if (fe_nothing->is_dominating())
696  else
697  // the FE_Nothing has no degrees of freedom and it is typically used
698  // in a context where we don't require any continuity along the
699  // interface
701  }
702 
703  Assert(false, ExcNotImplemented());
705 }
706 
707 
708 
709 template <int dim, int spacedim>
710 std::vector<std::pair<unsigned int, unsigned int>>
712  const FiniteElement<dim, spacedim> &fe_other) const
713 {
714  AssertDimension(dim, 2);
715 
716  if (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other) != nullptr)
717  {
718  // there should be exactly one single DoF of each FE at a vertex, and
719  // they should have identical value
720  return {{0U, 0U}};
721  }
722  else if (dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other) != nullptr)
723  {
724  // there should be exactly one single DoF of each FE at a vertex, and
725  // they should have identical value
726  return {{0U, 0U}};
727  }
728  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
729  {
730  // the FE_Nothing has no degrees of freedom, so there are no
731  // equivalencies to be recorded
732  return {};
733  }
734  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
735  {
736  // if the other element has no elements on faces at all,
737  // then it would be impossible to enforce any kind of
738  // continuity even if we knew exactly what kind of element
739  // we have -- simply because the other element declares
740  // that it is discontinuous because it has no DoFs on
741  // its faces. in that case, just state that we have no
742  // constraints to declare
743  return {};
744  }
745  else
746  {
747  Assert(false, ExcNotImplemented());
748  return {};
749  }
750 }
751 
752 
753 
754 template <int dim, int spacedim>
755 std::vector<std::pair<unsigned int, unsigned int>>
757  const FiniteElement<dim, spacedim> &fe_other) const
758 {
759  AssertDimension(dim, 2);
760  Assert(this->degree <= 2, ExcNotImplemented());
761 
762  if (const FE_SimplexP<dim, spacedim> *fe_p_other =
763  dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe_other))
764  {
765  // dofs are located along lines, so two dofs are identical if they are
766  // located at identical positions.
767  // Therefore, read the points in unit_support_points for the
768  // first coordinate direction. For FE_SimplexP, they are currently
769  // hard-coded and we iterate over points on the first line which begin
770  // after the 3 vertex points in the complete list of unit support points
771 
772  Assert(fe_p_other->degree <= 2, ExcNotImplemented());
773 
774  std::vector<std::pair<unsigned int, unsigned int>> identities;
775 
776  for (unsigned int i = 0; i < this->degree - 1; ++i)
777  for (unsigned int j = 0; j < fe_p_other->degree - 1; ++j)
778  if (std::fabs(this->unit_support_points[i + 3][0] -
779  fe_p_other->unit_support_points[i + 3][0]) < 1e-14)
780  identities.emplace_back(i, j);
781 
782  return identities;
783  }
784  else if (const FE_Q<dim, spacedim> *fe_q_other =
785  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
786  {
787  // dofs are located along lines, so two dofs are identical if they are
788  // located at identical positions. if we had only equidistant points, we
789  // could simply check for similarity like (i+1)*q == (j+1)*p, but we
790  // might have other support points (e.g. Gauss-Lobatto
791  // points). Therefore, read the points in unit_support_points for the
792  // first coordinate direction. For FE_Q, we take the lexicographic
793  // ordering of the line support points in the first direction (i.e.,
794  // x-direction), which we access between index 1 and p-1 (index 0 and p
795  // are vertex dofs). For FE_SimplexP, they are currently hard-coded and we
796  // iterate over points on the first line which begin after the 3 vertex
797  // points in the complete list of unit support points
798 
799  const std::vector<unsigned int> &index_map_inverse_q_other =
800  fe_q_other->get_poly_space_numbering_inverse();
801 
802  std::vector<std::pair<unsigned int, unsigned int>> identities;
803 
804  for (unsigned int i = 0; i < this->degree - 1; ++i)
805  for (unsigned int j = 0; j < fe_q_other->degree - 1; ++j)
806  if (std::fabs(this->unit_support_points[i + 3][0] -
807  fe_q_other->get_unit_support_points()
808  [index_map_inverse_q_other[j + 1]][0]) < 1e-14)
809  identities.emplace_back(i, j);
810 
811  return identities;
812  }
813  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
814  {
815  // the FE_Nothing has no degrees of freedom, so there are no
816  // equivalencies to be recorded
817  return {};
818  }
819  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
820  {
821  // if the other element has no elements on faces at all,
822  // then it would be impossible to enforce any kind of
823  // continuity even if we knew exactly what kind of element
824  // we have -- simply because the other element declares
825  // that it is discontinuous because it has no DoFs on
826  // its faces. in that case, just state that we have no
827  // constraints to declare
828  return {};
829  }
830  else
831  {
832  Assert(false, ExcNotImplemented());
833  return {};
834  }
835 }
836 
837 
838 
839 template <int dim, int spacedim>
841  : FE_SimplexPoly<dim, spacedim>(
842  BarycentricPolynomials<dim>::get_fe_p_basis(degree),
843  FiniteElementData<dim>(get_dpo_vector_fe_dgp(dim, degree),
844  ReferenceCells::get_simplex<dim>(),
845  1,
846  degree,
847  FiniteElementData<dim>::L2),
848  unit_support_points_fe_p<dim>(degree),
849  unit_face_support_points_fe_p<dim>(degree, FiniteElementData<dim>::H1),
850  constraints_fe_p<dim>(degree))
851 {}
852 
853 
854 
855 template <int dim, int spacedim>
856 std::unique_ptr<FiniteElement<dim, spacedim>>
858 {
859  return std::make_unique<FE_SimplexDGP<dim, spacedim>>(*this);
860 }
861 
862 
863 
864 template <int dim, int spacedim>
865 std::string
867 {
868  std::ostringstream namebuf;
869  namebuf << "FE_SimplexDGP<" << dim << ">(" << this->degree << ")";
870 
871  return namebuf.str();
872 }
873 
874 
875 template <int dim, int spacedim>
878  const FiniteElement<dim, spacedim> &fe_other,
879  const unsigned int codim) const
880 {
881  Assert(codim <= dim, ExcImpossibleInDim(dim));
882 
883  // vertex/line/face domination
884  // ---------------------------
885  if (codim > 0)
886  // this is a discontinuous element, so by definition there will
887  // be no constraints wherever this element comes together with
888  // any other kind of element
890 
891  // cell domination
892  // ---------------
893  if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
894  dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
895  {
896  if (this->degree < fe_dgp_other->degree)
898  else if (this->degree == fe_dgp_other->degree)
900  else
902  }
903  else if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
904  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
905  {
906  if (this->degree < fe_dgq_other->degree)
908  else if (this->degree == fe_dgq_other->degree)
910  else
912  }
913  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
914  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
915  {
916  if (fe_nothing->is_dominating())
918  else
919  // the FE_Nothing has no degrees of freedom and it is typically used
920  // in a context where we don't require any continuity along the
921  // interface
923  }
924 
925  Assert(false, ExcNotImplemented());
927 }
928 
929 
930 
931 template <int dim, int spacedim>
932 std::vector<std::pair<unsigned int, unsigned int>>
934  const FiniteElement<dim, spacedim> &fe_other) const
935 {
936  (void)fe_other;
937 
938  return {};
939 }
940 
941 
942 
943 template <int dim, int spacedim>
944 std::vector<std::pair<unsigned int, unsigned int>>
946  const FiniteElement<dim, spacedim> &fe_other) const
947 {
948  (void)fe_other;
949 
950  return {};
951 }
952 
953 // explicit instantiations
954 #include "fe_simplex_p.inst"
955 
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:2431
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2424
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:2412
size_type m() const
Definition: point.h:111
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)
Number * end_raw()
Number * begin_raw()
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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)