Reference documentation for deal.II version GIT 574c7c8486 2023-09-22 10:20: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_bernstein.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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 
21 
23 #include <deal.II/fe/fe_dgq.h>
24 #include <deal.II/fe/fe_nothing.h>
25 #include <deal.II/fe/fe_q.h>
26 #include <deal.II/fe/fe_tools.h>
27 
28 #include <memory>
29 #include <sstream>
30 #include <vector>
31 
32 
34 
35 
36 
37 template <int dim, int spacedim>
39  : FE_Q_Base<dim, spacedim>(this->renumber_bases(degree),
41  degree),
42  1,
43  degree,
44  FiniteElementData<dim>::H1),
45  std::vector<bool>(1, false))
46 {}
47 
48 
49 
50 template <int dim, int spacedim>
51 void
54  FullMatrix<double> &) const
55 {
56  // no interpolation possible. throw exception, as documentation says
58  false,
60 }
61 
62 
63 
64 template <int dim, int spacedim>
65 const FullMatrix<double> &
67  const unsigned int,
68  const RefinementCase<dim> &) const
69 {
70  AssertThrow(false,
72  // return dummy, nothing will happen because the base class FE_Q_Base
73  // implements lazy evaluation of those matrices
74  return this->restriction[0][0];
75 }
76 
77 
78 
79 template <int dim, int spacedim>
80 const FullMatrix<double> &
82  const unsigned int,
83  const RefinementCase<dim> &) const
84 {
85  AssertThrow(false,
87  // return dummy, nothing will happen because the base class FE_Q_Base
88  // implements lazy evaluation of those matrices
89  return this->prolongation[0][0];
90 }
91 
92 
93 
94 template <int dim, int spacedim>
95 void
97  const FiniteElement<dim, spacedim> &source_fe,
98  FullMatrix<double> &interpolation_matrix,
99  const unsigned int face_no) const
100 {
101  get_subface_interpolation_matrix(source_fe,
103  interpolation_matrix,
104  face_no);
105 }
106 
107 
108 template <int dim, int spacedim>
109 void
111  const FiniteElement<dim, spacedim> &x_source_fe,
112  const unsigned int subface,
113  FullMatrix<double> &interpolation_matrix,
114  const unsigned int face_no) const
115 {
116  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
117  ExcDimensionMismatch(interpolation_matrix.m(),
118  x_source_fe.n_dofs_per_face(face_no)));
119 
120  // see if source is a Bernstein element
121  if (const FE_Bernstein<dim, spacedim> *source_fe =
122  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&x_source_fe))
123  {
124  // have this test in here since a table of size 2x0 reports its size as
125  // 0x0
126  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
127  ExcDimensionMismatch(interpolation_matrix.n(),
128  this->n_dofs_per_face(face_no)));
129 
130  // Make sure that the element for which the DoFs should be constrained
131  // is the one with the higher polynomial degree. Actually the procedure
132  // will work also if this assertion is not satisfied. But the matrices
133  // produced in that case might lead to problems in the hp-procedures,
134  // which use this method.
135  Assert(
136  this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
137  (typename FiniteElement<dim,
138  spacedim>::ExcInterpolationNotImplemented()));
139 
140  const Quadrature<dim - 1> quad_face_support(
141  FE_Q<dim, spacedim>(QIterated<1>(QTrapezoid<1>(), source_fe->degree))
142  .get_unit_face_support_points(face_no));
143 
144  // Rule of thumb for FP accuracy, that can be expected for a given
145  // polynomial degree. This value is used to cut off values close to
146  // zero.
147  const double eps = 2e-13 * std::max(this->degree, source_fe->degree) *
148  std::max(dim - 1, 1);
149 
150  // compute the interpolation matrix by simply taking the value at the
151  // support points.
152  // TODO: Verify that all faces are the same with respect to
153  // these support points. Furthermore, check if something has to
154  // be done for the face orientation flag in 3d.
155  const Quadrature<dim> subface_quadrature =
156  subface == numbers::invalid_unsigned_int ?
158  quad_face_support,
159  0) :
161  quad_face_support,
162  0,
163  subface);
164 
165  for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
166  {
167  const Point<dim> &p = subface_quadrature.point(i);
168  for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
169  {
170  double matrix_entry =
171  this->shape_value(this->face_to_cell_index(j, 0), p);
172 
173  // Correct the interpolated value. I.e. if it is close to 1 or
174  // 0, make it exactly 1 or 0. Unfortunately, this is required to
175  // avoid problems with higher order elements.
176  if (std::fabs(matrix_entry - 1.0) < eps)
177  matrix_entry = 1.0;
178  if (std::fabs(matrix_entry) < eps)
179  matrix_entry = 0.0;
180 
181  interpolation_matrix(i, j) = matrix_entry;
182  }
183  }
184 
185 #ifdef DEBUG
186  // make sure that the row sum of each of the matrices is 1 at this
187  // point. this must be so since the shape functions sum up to 1
188  for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
189  {
190  double sum = 0.;
191 
192  for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
193  sum += interpolation_matrix(j, i);
194 
196  }
197 #endif
198  }
199  else
200  {
201  // When the incoming element is not FE_Bernstein we can just delegate to
202  // the base class to create the interpolation matrix.
204  x_source_fe, subface, interpolation_matrix, face_no);
205  }
206 }
207 
208 
209 
210 template <int dim, int spacedim>
211 bool
213 {
214  return true;
215 }
216 
217 
218 template <int dim, int spacedim>
219 std::vector<std::pair<unsigned int, unsigned int>>
221  const FiniteElement<dim, spacedim> &fe_other) const
222 {
223  // we can presently only compute these identities if both FEs are
224  // FE_Bernsteins or if the other one is an FE_Nothing. in the first case,
225  // there should be exactly one single DoF of each FE at a vertex, and they
226  // should have identical value
227  if (dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other) != nullptr)
228  {
229  return std::vector<std::pair<unsigned int, unsigned int>>(
230  1, std::make_pair(0U, 0U));
231  }
232  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
233  {
234  // the FE_Nothing has no degrees of freedom, so there are no
235  // equivalencies to be recorded
236  return std::vector<std::pair<unsigned int, unsigned int>>();
237  }
238  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
239  {
240  // if the other element has no elements on faces at all,
241  // then it would be impossible to enforce any kind of
242  // continuity even if we knew exactly what kind of element
243  // we have -- simply because the other element declares
244  // that it is discontinuous because it has no DoFs on
245  // its faces. in that case, just state that we have no
246  // constraints to declare
247  return std::vector<std::pair<unsigned int, unsigned int>>();
248  }
249  else
250  {
251  Assert(false, ExcNotImplemented());
252  return std::vector<std::pair<unsigned int, unsigned int>>();
253  }
254 }
255 
256 
257 template <int dim, int spacedim>
258 std::vector<std::pair<unsigned int, unsigned int>>
260  const FiniteElement<dim, spacedim> &) const
261 {
262  // Since this FE is not interpolatory but on the vertices, we can
263  // not identify dofs on lines and on quads even if there are dofs
264  // on lines and on quads.
265  //
266  // we also have nothing to say about interpolation to other finite
267  // elements. consequently, we never have anything to say at all
268  return std::vector<std::pair<unsigned int, unsigned int>>();
269 }
270 
271 
272 template <int dim, int spacedim>
273 std::vector<std::pair<unsigned int, unsigned int>>
276  const unsigned int) const
277 {
278  // Since this FE is not interpolatory but on the vertices, we can
279  // not identify dofs on lines and on quads even if there are dofs
280  // on lines and on quads.
281  //
282  // we also have nothing to say about interpolation to other finite
283  // elements. consequently, we never have anything to say at all
284  return std::vector<std::pair<unsigned int, unsigned int>>();
285 }
286 
287 
288 template <int dim, int spacedim>
291  const FiniteElement<dim, spacedim> &fe_other,
292  const unsigned int codim) const
293 {
294  Assert(codim <= dim, ExcImpossibleInDim(dim));
295 
296  // vertex/line/face domination
297  // (if fe_other is derived from FE_DGQ)
298  // ------------------------------------
299  if (codim > 0)
300  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
301  // there are no requirements between continuous and discontinuous elements
303 
304  // vertex/line/face domination
305  // (if fe_other is not derived from FE_DGQ)
306  // & cell domination
307  // ----------------------------------------
308  if (const FE_Bernstein<dim, spacedim> *fe_b_other =
309  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
310  {
311  if (this->degree < fe_b_other->degree)
313  else if (this->degree == fe_b_other->degree)
315  else
317  }
318  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
319  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
320  {
321  if (fe_nothing->is_dominating())
323  else
324  // the FE_Nothing has no degrees of freedom and it is typically used
325  // in a context where we don't require any continuity along the
326  // interface
328  }
329 
330  Assert(false, ExcNotImplemented());
332 }
333 
334 
335 template <int dim, int spacedim>
336 std::string
338 {
339  // note that the FETools::get_fe_by_name function depends on the
340  // particular format of the string this function returns, so they have to be
341  // kept in synch
342 
343  std::ostringstream namebuf;
344  namebuf << "FE_Bernstein<" << dim << ">(" << this->degree << ")";
345  return namebuf.str();
346 }
347 
348 
349 template <int dim, int spacedim>
350 std::unique_ptr<FiniteElement<dim, spacedim>>
352 {
353  return std::make_unique<FE_Bernstein<dim, spacedim>>(*this);
354 }
355 
356 
360 template <int dim, int spacedim>
361 std::vector<unsigned int>
363 {
364  AssertThrow(deg > 0, ExcMessage("FE_Bernstein needs to be of degree > 0."));
365  std::vector<unsigned int> dpo(dim + 1, 1U);
366  for (unsigned int i = 1; i < dpo.size(); ++i)
367  dpo[i] = dpo[i - 1] * (deg - 1);
368  return dpo;
369 }
370 
371 
372 template <int dim, int spacedim>
375 {
377  ::generate_complete_bernstein_basis<double>(deg));
378  tpp.set_numbering(FETools::hierarchic_to_lexicographic_numbering<dim>(deg));
379  return tpp;
380 }
381 
382 
383 // explicit instantiations
384 #include "fe_bernstein.inst"
385 
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_Bernstein(const unsigned int p)
Definition: fe_bernstein.cc:38
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_bernstein.cc:81
TensorProductPolynomials< dim > renumber_bases(const unsigned int degree)
virtual std::string get_name() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_bernstein.cc:52
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_bernstein.cc:66
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition: fe_bernstein.cc:96
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_dgq.h:113
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition: fe_q_base.cc:628
Definition: fe_q.h:551
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
const std::vector< Point< dim - 1 > > & get_unit_face_support_points(const unsigned int face_no=0) const
size_type n() const
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)
const Point< dim > & point(const unsigned int i) const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
void set_numbering(const std::vector< unsigned int > &renumber)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
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)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
Expression fabs(const Expression &x)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
static const char U
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
static const unsigned int invalid_unsigned_int
Definition: types.h:213