Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_bernstein.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2018 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 
17 #include <deal.II/base/polynomials_bernstein.h>
18 #include <deal.II/base/qprojector.h>
19 #include <deal.II/base/quadrature.h>
20 #include <deal.II/base/quadrature_lib.h>
21 #include <deal.II/base/std_cxx14/memory.h>
22 
23 #include <deal.II/fe/fe_bernstein.h>
24 #include <deal.II/fe/fe_dgq.h>
25 #include <deal.II/fe/fe_nothing.h>
26 #include <deal.II/fe/fe_q.h>
27 #include <deal.II/fe/fe_tools.h>
28 
29 #include <sstream>
30 #include <vector>
31 
32 
33 DEAL_II_NAMESPACE_OPEN
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) const
99 {
100  Assert(dim > 1, ExcImpossibleInDim(1));
101  get_subface_interpolation_matrix(source_fe,
103  interpolation_matrix);
104 }
105 
106 
107 template <int dim, int spacedim>
108 void
110  const FiniteElement<dim, spacedim> &x_source_fe,
111  const unsigned int subface,
112  FullMatrix<double> & interpolation_matrix) const
113 {
114  Assert(interpolation_matrix.m() == x_source_fe.dofs_per_face,
115  ExcDimensionMismatch(interpolation_matrix.m(),
116  x_source_fe.dofs_per_face));
117 
118  // see if source is a Bernstein element
119  if (const FE_Bernstein<dim, spacedim> *source_fe =
120  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&x_source_fe))
121  {
122  // have this test in here since a table of size 2x0 reports its size as
123  // 0x0
124  Assert(interpolation_matrix.n() == this->dofs_per_face,
125  ExcDimensionMismatch(interpolation_matrix.n(),
126  this->dofs_per_face));
127 
128  // Make sure that the element for which the DoFs should be constrained
129  // is the one with the higher polynomial degree. Actually the procedure
130  // will work also if this assertion is not satisfied. But the matrices
131  // produced in that case might lead to problems in the hp procedures,
132  // which use this method.
133  Assert(
134  this->dofs_per_face <= source_fe->dofs_per_face,
135  (typename FiniteElement<dim,
136  spacedim>::ExcInterpolationNotImplemented()));
137 
138  const Quadrature<dim - 1> quad_face_support(
139  FE_Q<dim, spacedim>(QIterated<1>(QTrapez<1>(), source_fe->degree))
141 
142  // Rule of thumb for FP accuracy, that can be expected for a given
143  // polynomial degree. This value is used to cut off values close to
144  // zero.
145  double eps =
146  2e-13 * std::max(this->degree, source_fe->degree) * (dim - 1);
147 
148  // compute the interpolation matrix by simply taking the value at the
149  // support points.
150  // TODO: Verify that all faces are the same with respect to
151  // these support points. Furthermore, check if something has to
152  // be done for the face orientation flag in 3D.
153  const Quadrature<dim> subface_quadrature =
154  subface == numbers::invalid_unsigned_int ?
155  QProjector<dim>::project_to_face(quad_face_support, 0) :
156  QProjector<dim>::project_to_subface(quad_face_support, 0, subface);
157 
158  for (unsigned int i = 0; i < source_fe->dofs_per_face; ++i)
159  {
160  const Point<dim> &p = subface_quadrature.point(i);
161  for (unsigned int j = 0; j < this->dofs_per_face; ++j)
162  {
163  double matrix_entry =
164  this->shape_value(this->face_to_cell_index(j, 0), p);
165 
166  // Correct the interpolated value. I.e. if it is close to 1 or
167  // 0, make it exactly 1 or 0. Unfortunately, this is required to
168  // avoid problems with higher order elements.
169  if (std::fabs(matrix_entry - 1.0) < eps)
170  matrix_entry = 1.0;
171  if (std::fabs(matrix_entry) < eps)
172  matrix_entry = 0.0;
173 
174  interpolation_matrix(i, j) = matrix_entry;
175  }
176  }
177 
178  // make sure that the row sum of each of the matrices is 1 at this
179  // point. this must be so since the shape functions sum up to 1
180  for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
181  {
182  double sum = 0.;
183 
184  for (unsigned int i = 0; i < this->dofs_per_face; ++i)
185  sum += interpolation_matrix(j, i);
186 
187  Assert(std::fabs(sum - 1) < eps, ExcInternalError());
188  }
189  }
190  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
191  {
192  // nothing to do here, the FE_Nothing has no degrees of freedom anyway
193  }
194  else
195  AssertThrow(
196  false,
197  (typename FiniteElement<dim,
198  spacedim>::ExcInterpolationNotImplemented()));
199 }
200 
201 
202 
203 template <int dim, int spacedim>
204 bool
206 {
207  return true;
208 }
209 
210 
211 template <int dim, int spacedim>
212 std::vector<std::pair<unsigned int, unsigned int>>
214  const FiniteElement<dim, spacedim> &fe_other) const
215 {
216  // we can presently only compute these identities if both FEs are
217  // FE_Bernsteins or if the other one is an FE_Nothing. in the first case,
218  // there should be exactly one single DoF of each FE at a vertex, and they
219  // should have identical value
220  if (dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other) != nullptr)
221  {
222  return std::vector<std::pair<unsigned int, unsigned int>>(
223  1, std::make_pair(0U, 0U));
224  }
225  else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
226  {
227  // the FE_Nothing has no degrees of freedom, so there are no
228  // equivalencies to be recorded
229  return std::vector<std::pair<unsigned int, unsigned int>>();
230  }
231  else if (fe_other.dofs_per_face == 0)
232  {
233  // if the other element has no elements on faces at all,
234  // then it would be impossible to enforce any kind of
235  // continuity even if we knew exactly what kind of element
236  // we have -- simply because the other element declares
237  // that it is discontinuous because it has no DoFs on
238  // its faces. in that case, just state that we have no
239  // constraints to declare
240  return std::vector<std::pair<unsigned int, unsigned int>>();
241  }
242  else
243  {
244  Assert(false, ExcNotImplemented());
245  return std::vector<std::pair<unsigned int, unsigned int>>();
246  }
247 }
248 
249 
250 template <int dim, int spacedim>
251 std::vector<std::pair<unsigned int, unsigned int>>
253  const FiniteElement<dim, spacedim> &) const
254 {
255  // Since this fe is not interpolatory but on the vertices, we can
256  // not identify dofs on lines and on quads even if there are dofs
257  // on lines and on quads.
258  //
259  // we also have nothing to say about interpolation to other finite
260  // elements. consequently, we never have anything to say at all
261  return std::vector<std::pair<unsigned int, unsigned int>>();
262 }
263 
264 
265 template <int dim, int spacedim>
266 std::vector<std::pair<unsigned int, unsigned int>>
268  const FiniteElement<dim, spacedim> &) const
269 {
270  // Since this fe is not interpolatory but on the vertices, we can
271  // not identify dofs on lines and on quads even if there are dofs
272  // on lines and on quads.
273  //
274  // we also have nothing to say about interpolation to other finite
275  // elements. consequently, we never have anything to say at all
276  return std::vector<std::pair<unsigned int, unsigned int>>();
277 }
278 
279 
280 template <int dim, int spacedim>
283  const FiniteElement<dim, spacedim> &fe_other,
284  const unsigned int codim) const
285 {
286  Assert(codim <= dim, ExcImpossibleInDim(dim));
287 
288  // vertex/line/face domination
289  // (if fe_other is derived from FE_DGQ)
290  // ------------------------------------
291  if (codim > 0)
292  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
293  // there are no requirements between continuous and discontinuous elements
295 
296  // vertex/line/face domination
297  // (if fe_other is not derived from FE_DGQ)
298  // & cell domination
299  // ----------------------------------------
300  if (const FE_Bernstein<dim, spacedim> *fe_b_other =
301  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
302  {
303  if (this->degree < fe_b_other->degree)
305  else if (this->degree == fe_b_other->degree)
307  else
309  }
310  else if (const FE_Nothing<dim, spacedim> *fe_nothing =
311  dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
312  {
313  if (fe_nothing->is_dominating())
315  else
316  // the FE_Nothing has no degrees of freedom and it is typically used
317  // in a context where we don't require any continuity along the
318  // interface
320  }
321 
322  Assert(false, ExcNotImplemented());
324 }
325 
326 
327 template <int dim, int spacedim>
328 std::string
330 {
331  // note that the FETools::get_fe_by_name function depends on the
332  // particular format of the string this function returns, so they have to be
333  // kept in synch
334 
335  std::ostringstream namebuf;
336  namebuf << "FE_Bernstein<" << dim << ">(" << this->degree << ")";
337  return namebuf.str();
338 }
339 
340 
341 template <int dim, int spacedim>
342 std::unique_ptr<FiniteElement<dim, spacedim>>
344 {
345  return std_cxx14::make_unique<FE_Bernstein<dim, spacedim>>(*this);
346 }
347 
348 
352 template <int dim, int spacedim>
353 std::vector<unsigned int>
355 {
356  AssertThrow(deg > 0, ExcMessage("FE_Bernstein needs to be of degree > 0."));
357  std::vector<unsigned int> dpo(dim + 1, 1U);
358  for (unsigned int i = 1; i < dpo.size(); ++i)
359  dpo[i] = dpo[i - 1] * (deg - 1);
360  return dpo;
361 }
362 
363 
364 template <int dim, int spacedim>
367 {
369  ::generate_complete_bernstein_basis<double>(deg));
370  std::vector<unsigned int> renumber(Utilities::fixed_power<dim>(deg + 1));
371  const FiniteElementData<dim> fe(this->get_dpo_vector(deg), 1, deg);
373  tpp.set_numbering(renumber);
374  return tpp;
375 }
376 
377 
378 // explicit instantiations
379 #include "fe_bernstein.inst"
380 
381 DEAL_II_NAMESPACE_CLOSE
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:173
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)
const std::vector< Point< dim - 1 > > & get_unit_face_support_points() const
Definition: fe.cc:1069
Definition: fe_dgq.h:109
void hierarchic_to_lexicographic_numbering(unsigned int degree, std::vector< unsigned int > &h2l)
void set_numbering(const std::vector< unsigned int > &renumber)
STL namespace.
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcMessage(std::string arg1)
Definition: fe_q.h:554
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
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
virtual std::string get_name() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
const unsigned int dofs_per_face
Definition: fe_base.h:275
static ::ExceptionBase & ExcNotImplemented()
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_bernstein.cc:96
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
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 override
static ::ExceptionBase & ExcInternalError()