Reference documentation for deal.II version GIT relicensing-214-g6e74dec06b 2024-03-27 18:10:01+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\}}\)
Loading...
Searching...
No Matches
fe_bernstein.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
20
22#include <deal.II/fe/fe_dgq.h>
24#include <deal.II/fe/fe_q.h>
25#include <deal.II/fe/fe_tools.h>
26
27#include <memory>
28#include <sstream>
29#include <vector>
30
31
33
34
35
36template <int dim, int spacedim>
38 : FE_Q_Base<dim, spacedim>(this->renumber_bases(degree),
39 FiniteElementData<dim>(this->get_dpo_vector(
40 degree),
41 1,
42 degree,
43 FiniteElementData<dim>::H1),
44 std::vector<bool>(1, false))
45{}
46
47
48
49template <int dim, int spacedim>
50void
53 FullMatrix<double> &) const
54{
55 // no interpolation possible. throw exception, as documentation says
57 false,
59}
60
61
62
63template <int dim, int spacedim>
66 const unsigned int,
67 const RefinementCase<dim> &) const
68{
69 AssertThrow(false,
71 // return dummy, nothing will happen because the base class FE_Q_Base
72 // implements lazy evaluation of those matrices
73 return this->restriction[0][0];
74}
75
76
77
78template <int dim, int spacedim>
81 const unsigned int,
82 const RefinementCase<dim> &) const
83{
84 AssertThrow(false,
86 // return dummy, nothing will happen because the base class FE_Q_Base
87 // implements lazy evaluation of those matrices
88 return this->prolongation[0][0];
89}
90
91
92
93template <int dim, int spacedim>
94void
96 const FiniteElement<dim, spacedim> &source_fe,
97 FullMatrix<double> &interpolation_matrix,
98 const unsigned int face_no) const
99{
100 get_subface_interpolation_matrix(source_fe,
102 interpolation_matrix,
103 face_no);
104}
105
106
107template <int dim, int spacedim>
108void
110 const FiniteElement<dim, spacedim> &x_source_fe,
111 const unsigned int subface,
112 FullMatrix<double> &interpolation_matrix,
113 const unsigned int face_no) const
114{
115 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
116 ExcDimensionMismatch(interpolation_matrix.m(),
117 x_source_fe.n_dofs_per_face(face_no)));
118
119 // see if source is a Bernstein element
120 if (const FE_Bernstein<dim, spacedim> *source_fe =
121 dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&x_source_fe))
122 {
123 // have this test in here since a table of size 2x0 reports its size as
124 // 0x0
125 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
126 ExcDimensionMismatch(interpolation_matrix.n(),
127 this->n_dofs_per_face(face_no)));
128
129 // Make sure that the element for which the DoFs should be constrained
130 // is the one with the higher polynomial degree. Actually the procedure
131 // will work also if this assertion is not satisfied. But the matrices
132 // produced in that case might lead to problems in the hp-procedures,
133 // which use this method.
134 Assert(
135 this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
136 (typename FiniteElement<dim,
137 spacedim>::ExcInterpolationNotImplemented()));
138
139 const Quadrature<dim - 1> quad_face_support(
140 FE_Q<dim, spacedim>(QIterated<1>(QTrapezoid<1>(), source_fe->degree))
142
143 // Rule of thumb for FP accuracy, that can be expected for a given
144 // polynomial degree. This value is used to cut off values close to
145 // zero.
146 const double eps = 2e-13 * std::max(this->degree, source_fe->degree) *
147 std::max(dim - 1, 1);
148
149 // compute the interpolation matrix by simply taking the value at the
150 // support points.
151 // TODO: Verify that all faces are the same with respect to
152 // these support points. Furthermore, check if something has to
153 // be done for the face orientation flag in 3d.
154 const Quadrature<dim> subface_quadrature =
156 QProjector<dim>::project_to_face(this->reference_cell(),
157 quad_face_support,
158 0) :
159 QProjector<dim>::project_to_subface(this->reference_cell(),
160 quad_face_support,
161 0,
162 subface);
163
164 for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
165 {
166 const Point<dim> &p = subface_quadrature.point(i);
167 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
168 {
169 double matrix_entry =
170 this->shape_value(this->face_to_cell_index(j, 0), p);
171
172 // Correct the interpolated value. I.e. if it is close to 1 or
173 // 0, make it exactly 1 or 0. Unfortunately, this is required to
174 // avoid problems with higher order elements.
175 if (std::fabs(matrix_entry - 1.0) < eps)
176 matrix_entry = 1.0;
177 if (std::fabs(matrix_entry) < eps)
178 matrix_entry = 0.0;
179
180 interpolation_matrix(i, j) = matrix_entry;
181 }
182 }
183
184#ifdef DEBUG
185 // make sure that the row sum of each of the matrices is 1 at this
186 // point. this must be so since the shape functions sum up to 1
187 for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
188 {
189 double sum = 0.;
190
191 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
192 sum += interpolation_matrix(j, i);
193
194 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
195 }
196#endif
197 }
198 else
199 {
200 // When the incoming element is not FE_Bernstein we can just delegate to
201 // the base class to create the interpolation matrix.
203 x_source_fe, subface, interpolation_matrix, face_no);
204 }
205}
206
207
208
209template <int dim, int spacedim>
210bool
215
216
217template <int dim, int spacedim>
218std::vector<std::pair<unsigned int, unsigned int>>
220 const FiniteElement<dim, spacedim> &fe_other) const
221{
222 // we can presently only compute these identities if both FEs are
223 // FE_Bernsteins or if the other one is an FE_Nothing. in the first case,
224 // there should be exactly one single DoF of each FE at a vertex, and they
225 // should have identical value
226 if (dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other) != nullptr)
227 {
228 return std::vector<std::pair<unsigned int, unsigned int>>(
229 1, std::make_pair(0U, 0U));
230 }
231 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
232 {
233 // the FE_Nothing has no degrees of freedom, so there are no
234 // equivalencies to be recorded
235 return std::vector<std::pair<unsigned int, unsigned int>>();
236 }
237 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
238 {
239 // if the other element has no elements on faces at all,
240 // then it would be impossible to enforce any kind of
241 // continuity even if we knew exactly what kind of element
242 // we have -- simply because the other element declares
243 // that it is discontinuous because it has no DoFs on
244 // its faces. in that case, just state that we have no
245 // constraints to declare
246 return std::vector<std::pair<unsigned int, unsigned int>>();
247 }
248 else
249 {
251 return std::vector<std::pair<unsigned int, unsigned int>>();
252 }
253}
254
255
256template <int dim, int spacedim>
257std::vector<std::pair<unsigned int, unsigned int>>
259 const FiniteElement<dim, spacedim> &) const
260{
261 // Since this FE is not interpolatory but on the vertices, we can
262 // not identify dofs on lines and on quads even if there are dofs
263 // on lines and on quads.
264 //
265 // we also have nothing to say about interpolation to other finite
266 // elements. consequently, we never have anything to say at all
267 return std::vector<std::pair<unsigned int, unsigned int>>();
268}
269
270
271template <int dim, int spacedim>
272std::vector<std::pair<unsigned int, unsigned int>>
275 const unsigned int) const
276{
277 // Since this FE is not interpolatory but on the vertices, we can
278 // not identify dofs on lines and on quads even if there are dofs
279 // on lines and on quads.
280 //
281 // we also have nothing to say about interpolation to other finite
282 // elements. consequently, we never have anything to say at all
283 return std::vector<std::pair<unsigned int, unsigned int>>();
284}
285
286
287template <int dim, int spacedim>
290 const FiniteElement<dim, spacedim> &fe_other,
291 const unsigned int codim) const
292{
293 Assert(codim <= dim, ExcImpossibleInDim(dim));
294
295 // vertex/line/face domination
296 // (if fe_other is derived from FE_DGQ)
297 // ------------------------------------
298 if (codim > 0)
299 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
300 // there are no requirements between continuous and discontinuous elements
302
303 // vertex/line/face domination
304 // (if fe_other is not derived from FE_DGQ)
305 // & cell domination
306 // ----------------------------------------
307 if (const FE_Bernstein<dim, spacedim> *fe_b_other =
308 dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
309 {
310 if (this->degree < fe_b_other->degree)
312 else if (this->degree == fe_b_other->degree)
314 else
316 }
317 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
318 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
319 {
320 if (fe_nothing->is_dominating())
322 else
323 // the FE_Nothing has no degrees of freedom and it is typically used
324 // in a context where we don't require any continuity along the
325 // interface
327 }
328
331}
332
333
334template <int dim, int spacedim>
335std::string
337{
338 // note that the FETools::get_fe_by_name function depends on the
339 // particular format of the string this function returns, so they have to be
340 // kept in synch
341
342 std::ostringstream namebuf;
343 namebuf << "FE_Bernstein<" << dim << ">(" << this->degree << ")";
344 return namebuf.str();
345}
346
347
348template <int dim, int spacedim>
349std::unique_ptr<FiniteElement<dim, spacedim>>
351{
352 return std::make_unique<FE_Bernstein<dim, spacedim>>(*this);
353}
354
355
359template <int dim, int spacedim>
360std::vector<unsigned int>
362{
363 AssertThrow(deg > 0, ExcMessage("FE_Bernstein needs to be of degree > 0."));
364 std::vector<unsigned int> dpo(dim + 1, 1U);
365 for (unsigned int i = 1; i < dpo.size(); ++i)
366 dpo[i] = dpo[i - 1] * (deg - 1);
367 return dpo;
368}
369
370
371template <int dim, int spacedim>
374{
376 ::generate_complete_bernstein_basis<double>(deg));
377 tpp.set_numbering(FETools::hierarchic_to_lexicographic_numbering<dim>(deg));
378 return tpp;
379}
380
381
382// explicit instantiations
383#include "fe_bernstein.inst"
384
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:550
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: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)
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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)
#define DEAL_II_NOT_IMPLEMENTED()
static const unsigned int invalid_unsigned_int
Definition types.h:220
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)