Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3094-gbcbfbe1d73 2025-04-21 14:50:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 =
157 this->reference_cell(),
158 quad_face_support,
159 0,
162 this->reference_cell(),
163 quad_face_support,
164 0,
165 subface,
168
169 for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
170 {
171 const Point<dim> &p = subface_quadrature.point(i);
172 for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
173 {
174 double matrix_entry =
175 this->shape_value(this->face_to_cell_index(j, 0), p);
176
177 // Correct the interpolated value. I.e. if it is close to 1 or
178 // 0, make it exactly 1 or 0. Unfortunately, this is required to
179 // avoid problems with higher order elements.
180 if (std::fabs(matrix_entry - 1.0) < eps)
181 matrix_entry = 1.0;
182 if (std::fabs(matrix_entry) < eps)
183 matrix_entry = 0.0;
184
185 interpolation_matrix(i, j) = matrix_entry;
186 }
187 }
188
189 if constexpr (running_in_debug_mode())
190 {
191 // make sure that the row sum of each of the matrices is 1 at this
192 // point. this must be so since the shape functions sum up to 1
193 for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
194 {
195 double sum = 0.;
196
197 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
198 sum += interpolation_matrix(j, i);
199
200 Assert(std::fabs(sum - 1) < eps, ExcInternalError());
201 }
202 }
203 }
204 else
205 {
206 // When the incoming element is not FE_Bernstein we can just delegate to
207 // the base class to create the interpolation matrix.
209 x_source_fe, subface, interpolation_matrix, face_no);
210 }
211}
212
213
214
215template <int dim, int spacedim>
216bool
221
222
223template <int dim, int spacedim>
224std::vector<std::pair<unsigned int, unsigned int>>
226 const FiniteElement<dim, spacedim> &fe_other) const
227{
228 // we can presently only compute these identities if both FEs are
229 // FE_Bernsteins or if the other one is an FE_Nothing. in the first case,
230 // there should be exactly one single DoF of each FE at a vertex, and they
231 // should have identical value
232 if (dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other) != nullptr)
233 {
234 return std::vector<std::pair<unsigned int, unsigned int>>(
235 1, std::make_pair(0U, 0U));
236 }
237 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
238 {
239 // the FE_Nothing has no degrees of freedom, so there are no
240 // equivalencies to be recorded
241 return std::vector<std::pair<unsigned int, unsigned int>>();
242 }
243 else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
244 {
245 // if the other element has no elements on faces at all,
246 // then it would be impossible to enforce any kind of
247 // continuity even if we knew exactly what kind of element
248 // we have -- simply because the other element declares
249 // that it is discontinuous because it has no DoFs on
250 // its faces. in that case, just state that we have no
251 // constraints to declare
252 return std::vector<std::pair<unsigned int, unsigned int>>();
253 }
254 else
255 {
257 return std::vector<std::pair<unsigned int, unsigned int>>();
258 }
259}
260
261
262template <int dim, int spacedim>
263std::vector<std::pair<unsigned int, unsigned int>>
265 const FiniteElement<dim, spacedim> &) const
266{
267 // Since this FE is not interpolatory but on the vertices, we can
268 // not identify dofs on lines and on quads even if there are dofs
269 // on lines and on quads.
270 //
271 // we also have nothing to say about interpolation to other finite
272 // elements. consequently, we never have anything to say at all
273 return std::vector<std::pair<unsigned int, unsigned int>>();
274}
275
276
277template <int dim, int spacedim>
278std::vector<std::pair<unsigned int, unsigned int>>
281 const unsigned int) const
282{
283 // Since this FE is not interpolatory but on the vertices, we can
284 // not identify dofs on lines and on quads even if there are dofs
285 // on lines and on quads.
286 //
287 // we also have nothing to say about interpolation to other finite
288 // elements. consequently, we never have anything to say at all
289 return std::vector<std::pair<unsigned int, unsigned int>>();
290}
291
292
293template <int dim, int spacedim>
296 const FiniteElement<dim, spacedim> &fe_other,
297 const unsigned int codim) const
298{
299 Assert(codim <= dim, ExcImpossibleInDim(dim));
300
301 // vertex/line/face domination
302 // (if fe_other is derived from FE_DGQ)
303 // ------------------------------------
304 if (codim > 0)
305 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
306 // there are no requirements between continuous and discontinuous elements
308
309 // vertex/line/face domination
310 // (if fe_other is not derived from FE_DGQ)
311 // & cell domination
312 // ----------------------------------------
313 if (const FE_Bernstein<dim, spacedim> *fe_b_other =
314 dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
315 {
316 if (this->degree < fe_b_other->degree)
318 else if (this->degree == fe_b_other->degree)
320 else
322 }
323 else if (const FE_Nothing<dim, spacedim> *fe_nothing =
324 dynamic_cast<const FE_Nothing<dim, spacedim> *>(&fe_other))
325 {
326 if (fe_nothing->is_dominating())
328 else
329 // the FE_Nothing has no degrees of freedom and it is typically used
330 // in a context where we don't require any continuity along the
331 // interface
333 }
334
337}
338
339
340template <int dim, int spacedim>
341std::string
343{
344 // note that the FETools::get_fe_by_name function depends on the
345 // particular format of the string this function returns, so they have to be
346 // kept in synch
347
348 std::ostringstream namebuf;
349 namebuf << "FE_Bernstein<" << Utilities::dim_string(dim, spacedim) << ">("
350 << this->degree << ")";
351 return namebuf.str();
352}
353
354
355template <int dim, int spacedim>
356std::unique_ptr<FiniteElement<dim, spacedim>>
358{
359 return std::make_unique<FE_Bernstein<dim, spacedim>>(*this);
360}
361
362
366template <int dim, int spacedim>
367std::vector<unsigned int>
369{
370 AssertThrow(deg > 0, ExcMessage("FE_Bernstein needs to be of degree > 0."));
371 std::vector<unsigned int> dpo(dim + 1, 1U);
372 for (unsigned int i = 1; i < dpo.size(); ++i)
373 dpo[i] = dpo[i - 1] * (deg - 1);
374 return dpo;
375}
376
377
378template <int dim, int spacedim>
381{
383 ::generate_complete_bernstein_basis<double>(deg));
384 tpp.set_numbering(FETools::hierarchic_to_lexicographic_numbering<dim>(deg));
385 return tpp;
386}
387
388
389// explicit instantiations
390#include "fe/fe_bernstein.inst"
391
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:554
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:113
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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_NOT_IMPLEMENTED()
#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)
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:551
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
constexpr types::geometric_orientation default_geometric_orientation
Definition types.h:352
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)