Reference documentation for deal.II version Git a285ab78cc 2020-08-04 19:09:24 +0200
\(\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 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
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<TensorProductPolynomials<dim>, dim, spacedim>(
40  this->renumber_bases(degree),
41  FiniteElementData<dim>(this->get_dpo_vector(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  Assert(dim > 1, ExcImpossibleInDim(1));
104  interpolation_matrix,
105  face_no);
106 }
107 
108 
109 template <int dim, int spacedim>
110 void
112  const FiniteElement<dim, spacedim> &x_source_fe,
113  const unsigned int subface,
114  FullMatrix<double> & interpolation_matrix,
115  const unsigned int) const
116 {
117  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
118  ExcDimensionMismatch(interpolation_matrix.m(),
119  x_source_fe.n_dofs_per_face()));
120 
121  // see if source is a Bernstein element
122  if (const FE_Bernstein<dim, spacedim> *source_fe =
123  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&x_source_fe))
124  {
125  // have this test in here since a table of size 2x0 reports its size as
126  // 0x0
127  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
128  ExcDimensionMismatch(interpolation_matrix.n(),
129  this->n_dofs_per_face()));
130 
131  // Make sure that the element for which the DoFs should be constrained
132  // is the one with the higher polynomial degree. Actually the procedure
133  // will work also if this assertion is not satisfied. But the matrices
134  // produced in that case might lead to problems in the hp procedures,
135  // which use this method.
136  Assert(
137  this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
138  (typename FiniteElement<dim,
140 
141  const Quadrature<dim - 1> quad_face_support(
142  FE_Q<dim, spacedim>(QIterated<1>(QTrapez<1>(), source_fe->degree))
144 
145  // Rule of thumb for FP accuracy, that can be expected for a given
146  // polynomial degree. This value is used to cut off values close to
147  // zero.
148  double eps =
149  2e-13 * std::max(this->degree, source_fe->degree) * (dim - 1);
150 
151  // compute the interpolation matrix by simply taking the value at the
152  // support points.
153  // TODO: Verify that all faces are the same with respect to
154  // these support points. Furthermore, check if something has to
155  // be done for the face orientation flag in 3D.
156  const Quadrature<dim> subface_quadrature =
157  subface == numbers::invalid_unsigned_int ?
159  quad_face_support,
160  0) :
162  quad_face_support,
163  0,
164  subface);
165 
166  for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
167  {
168  const Point<dim> &p = subface_quadrature.point(i);
169  for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
170  {
171  double matrix_entry =
172  this->shape_value(this->face_to_cell_index(j, 0), p);
173 
174  // Correct the interpolated value. I.e. if it is close to 1 or
175  // 0, make it exactly 1 or 0. Unfortunately, this is required to
176  // avoid problems with higher order elements.
177  if (std::fabs(matrix_entry - 1.0) < eps)
178  matrix_entry = 1.0;
179  if (std::fabs(matrix_entry) < eps)
180  matrix_entry = 0.0;
181 
182  interpolation_matrix(i, j) = matrix_entry;
183  }
184  }
185 
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(); ++j)
189  {
190  double sum = 0.;
191 
192  for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
193  sum += interpolation_matrix(j, i);
194 
195  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
196  }
197  }
198  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
199  {
200  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
201  }
202  else
203  AssertThrow(
204  false,
205  (typename FiniteElement<dim,
206  spacedim>::ExcInterpolationNotImplemented()));
207 }
208 
209 
210 
211 template <int dim, int spacedim>
212 bool
214 {
215  return true;
216 }
217 
218 
219 template <int dim, int spacedim>
220 std::vector<std::pair<unsigned int, unsigned int>>
222  const FiniteElement<dim, spacedim> &fe_other) const
223 {
224  // we can presently only compute these identities if both FEs are
225  // FE_Bernsteins or if the other one is an FE_Nothing. in the first case,
226  // there should be exactly one single DoF of each FE at a vertex, and they
227  // should have identical value
228  if (dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other) != nullptr)
229  {
230  return std::vector<std::pair<unsigned int, unsigned int>>(
231  1, std::make_pair(0U, 0U));
232  }
233  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
234  {
235  // the FE_Nothing has no degrees of freedom, so there are no
236  // equivalencies to be recorded
237  return std::vector<std::pair<unsigned int, unsigned int>>();
238  }
239  else if (fe_other.n_dofs_per_face() == 0)
240  {
241  // if the other element has no elements on faces at all,
242  // then it would be impossible to enforce any kind of
243  // continuity even if we knew exactly what kind of element
244  // we have -- simply because the other element declares
245  // that it is discontinuous because it has no DoFs on
246  // its faces. in that case, just state that we have no
247  // constraints to declare
248  return std::vector<std::pair<unsigned int, unsigned int>>();
249  }
250  else
251  {
252  Assert(false, ExcNotImplemented());
253  return std::vector<std::pair<unsigned int, unsigned int>>();
254  }
255 }
256 
257 
258 template <int dim, int spacedim>
259 std::vector<std::pair<unsigned int, unsigned int>>
261  const FiniteElement<dim, spacedim> &) const
262 {
263  // Since this fe is not interpolatory but on the vertices, we can
264  // not identify dofs on lines and on quads even if there are dofs
265  // on lines and on quads.
266  //
267  // we also have nothing to say about interpolation to other finite
268  // elements. consequently, we never have anything to say at all
269  return std::vector<std::pair<unsigned int, unsigned int>>();
270 }
271 
272 
273 template <int dim, int spacedim>
274 std::vector<std::pair<unsigned int, unsigned int>>
277  const unsigned int) const
278 {
279  // Since this fe is not interpolatory but on the vertices, we can
280  // not identify dofs on lines and on quads even if there are dofs
281  // on lines and on quads.
282  //
283  // we also have nothing to say about interpolation to other finite
284  // elements. consequently, we never have anything to say at all
285  return std::vector<std::pair<unsigned int, unsigned int>>();
286 }
287 
288 
289 template <int dim, int spacedim>
292  const FiniteElement<dim, spacedim> &fe_other,
293  const unsigned int codim) const
294 {
295  Assert(codim <= dim, ExcImpossibleInDim(dim));
296 
297  // vertex/line/face domination
298  // (if fe_other is derived from FE_DGQ)
299  // ------------------------------------
300  if (codim > 0)
301  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
302  // there are no requirements between continuous and discontinuous elements
304 
305  // vertex/line/face domination
306  // (if fe_other is not derived from FE_DGQ)
307  // & cell domination
308  // ----------------------------------------
309  if (const FE_Bernstein<dim, spacedim> *fe_b_other =
310  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
311  {
312  if (this->degree < fe_b_other->degree)
314  else if (this->degree == fe_b_other->degree)
316  else
318  }
319  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
320  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
321  {
322  if (fe_nothing->is_dominating())
324  else
325  // the FE_Nothing has no degrees of freedom and it is typically used
326  // in a context where we don't require any continuity along the
327  // interface
329  }
330 
331  Assert(false, ExcNotImplemented());
333 }
334 
335 
336 template <int dim, int spacedim>
337 std::string
339 {
340  // note that the FETools::get_fe_by_name function depends on the
341  // particular format of the string this function returns, so they have to be
342  // kept in synch
343 
344  std::ostringstream namebuf;
345  namebuf << "FE_Bernstein<" << dim << ">(" << this->degree << ")";
346  return namebuf.str();
347 }
348 
349 
350 template <int dim, int spacedim>
351 std::unique_ptr<FiniteElement<dim, spacedim>>
353 {
354  return std::make_unique<FE_Bernstein<dim, spacedim>>(*this);
355 }
356 
357 
361 template <int dim, int spacedim>
362 std::vector<unsigned int>
364 {
365  AssertThrow(deg > 0, ExcMessage("FE_Bernstein needs to be of degree > 0."));
366  std::vector<unsigned int> dpo(dim + 1, 1U);
367  for (unsigned int i = 1; i < dpo.size(); ++i)
368  dpo[i] = dpo[i - 1] * (deg - 1);
369  return dpo;
370 }
371 
372 
373 template <int dim, int spacedim>
376 {
378  ::generate_complete_bernstein_basis<double>(deg));
379  tpp.set_numbering(FETools::hierarchic_to_lexicographic_numbering<dim>(deg));
380  return tpp;
381 }
382 
383 
384 // explicit instantiations
385 #include "fe_bernstein.inst"
386 
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_bernstein.cc:52
size_type m() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2401
FE_Bernstein(const unsigned int p)
Definition: fe_bernstein.cc:38
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
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
const unsigned int degree
Definition: fe_base.h:431
void set_numbering(const std::vector< unsigned int > &renumber)
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
STL namespace.
static const char U
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
static ::ExceptionBase & ExcInterpolationNotImplemented()
ReferenceCell::Type reference_cell_type() const
virtual double shape_value(const unsigned int i, const Point< dim > &p) const override
static ::ExceptionBase & ExcMessage(std::string arg1)
Definition: fe_q.h:548
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2415
T sum(const T &t, const MPI_Comm &mpi_communicator)
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:1071
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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 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 override
Definition: fe_q_base.cc:1048
Expression fabs(const Expression &x)
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
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
TensorProductPolynomials< dim > renumber_bases(const unsigned int degree)
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
static void project_to_face(const SubQuadrature &quadrature, const unsigned int face_no, std::vector< Point< dim >> &q_points)
virtual bool hp_constraints_are_implemented() const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
virtual std::string get_name() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
static ::ExceptionBase & ExcNotImplemented()
T max(const T &t, const MPI_Comm &mpi_communicator)
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
static ::ExceptionBase & ExcInternalError()