Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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>
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
37template <int dim, int spacedim>
39 : FE_Q_Base<dim, spacedim>(this->renumber_bases(degree),
40 FiniteElementData<dim>(this->get_dpo_vector(
41 degree),
42 1,
43 degree,
44 FiniteElementData<dim>::H1),
45 std::vector<bool>(1, false))
46{}
47
48
49
50template <int dim, int spacedim>
51void
54 FullMatrix<double> &) const
55{
56 // no interpolation possible. throw exception, as documentation says
58 false,
60}
61
62
63
64template <int dim, int spacedim>
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
79template <int dim, int spacedim>
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
94template <int dim, int spacedim>
95void
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
108template <int dim, int spacedim>
109void
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))
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 =
157 QProjector<dim>::project_to_face(this->reference_cell(),
158 quad_face_support,
159 0) :
160 QProjector<dim>::project_to_subface(this->reference_cell(),
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
195 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
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
210template <int dim, int spacedim>
211bool
213{
214 return true;
215}
216
217
218template <int dim, int spacedim>
219std::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
257template <int dim, int spacedim>
258std::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
272template <int dim, int spacedim>
273std::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
288template <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
335template <int dim, int spacedim>
336std::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
349template <int dim, int spacedim>
350std::unique_ptr<FiniteElement<dim, spacedim>>
352{
353 return std::make_unique<FE_Bernstein<dim, spacedim>>(*this);
354}
355
356
360template <int dim, int spacedim>
361std::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
372template <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)
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
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
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
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
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
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.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
void set_numbering(const std::vector< unsigned int > &renumber)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
static const unsigned int invalid_unsigned_int
Definition types.h:213
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)