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_bdm.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2003 - 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
20
22
23#include <deal.II/fe/fe.h>
24#include <deal.II/fe/fe_bdm.h>
25#include <deal.II/fe/fe_tools.h>
27#include <deal.II/fe/mapping.h>
28
29#include <deal.II/grid/tria.h>
31
32#include <iostream>
33#include <memory>
34#include <sstream>
35
36
38
39// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
40// adjust_line_dof_index_for_line_orientation_table fields, and write tests
41// similar to bits/face_orientation_and_fe_q_*
42
43template <int dim>
44FE_BDM<dim>::FE_BDM(const unsigned int deg)
45 : FE_PolyTensor<dim>(
46 PolynomialsBDM<dim>(deg),
47 FiniteElementData<dim>(get_dpo_vector(deg),
48 dim,
49 deg + 1,
50 FiniteElementData<dim>::Hdiv),
51 get_ria_vector(deg),
52 std::vector<ComponentMask>(PolynomialsBDM<dim>::n_polynomials(deg),
53 std::vector<bool>(dim, true)))
54{
55 Assert(dim >= 2, ExcImpossibleInDim(dim));
56 Assert(
57 deg > 0,
59 "Lowest order BDM element are degree 1, but you asked for degree 0"));
60
61 const unsigned int n_dofs = this->n_dofs_per_cell();
62
63 this->mapping_kind = {mapping_bdm};
64 // These must be done first, since
65 // they change the evaluation of
66 // basis functions
67
68 // Set up the generalized support
69 // points
71
72 // Now compute the inverse node matrix, generating the correct
73 // basis functions from the raw ones. For a discussion of what
74 // exactly happens here, see FETools::compute_node_matrix.
76 this->inverse_node_matrix.reinit(n_dofs, n_dofs);
78 // From now on, the shape functions provided by FiniteElement::shape_value
79 // and similar functions will be the correct ones, not
80 // the raw shape functions from the polynomial space anymore.
81
82 // Embedding errors become pretty large, so we just replace the
83 // regular threshold in both "computing_..." functions by 1.
86
88 const unsigned int face_no = 0;
89
91 for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
92 face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
93 this->n_dofs_per_face(face_no));
94 FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.);
95 this->interface_constraints.reinit((1 << (dim - 1)) *
96 this->n_dofs_per_face(face_no),
97 this->n_dofs_per_face(face_no));
98 unsigned int target_row = 0;
99 for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
100 for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
101 {
102 for (unsigned int j = 0; j < face_embeddings[d].n(); ++j)
103 this->interface_constraints(target_row, j) = face_embeddings[d](i, j);
104 ++target_row;
105 }
106
107 // We need to initialize the dof permutation table and the one for the sign
108 // change.
110}
111
112
113template <int dim>
114void
116{
117 // for 1d and 2d, do nothing
118 if (dim < 3)
119 return;
120
121 // TODO: Implement this for this class
122 return;
123}
124
125
126template <int dim>
127std::string
129{
130 // note that the
131 // FETools::get_fe_by_name
132 // function depends on the
133 // particular format of the string
134 // this function returns, so they
135 // have to be kept in synch
136
137 // note that this->degree is the maximal
138 // polynomial degree and is thus one higher
139 // than the argument given to the
140 // constructor
141 std::ostringstream namebuf;
142 namebuf << "FE_BDM<" << dim << ">(" << this->degree - 1 << ")";
143
144 return namebuf.str();
145}
146
147
148template <int dim>
149std::unique_ptr<FiniteElement<dim, dim>>
151{
152 return std::make_unique<FE_BDM<dim>>(*this);
153}
154
155
156
157template <int dim>
158void
160 const std::vector<Vector<double>> &support_point_values,
161 std::vector<double> & nodal_values) const
162{
163 Assert(support_point_values.size() == this->generalized_support_points.size(),
164 ExcDimensionMismatch(support_point_values.size(),
165 this->generalized_support_points.size()));
166 AssertDimension(support_point_values[0].size(), dim);
167 Assert(nodal_values.size() == this->n_dofs_per_cell(),
168 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
169
170 // First do interpolation on faces. There, the component evaluated
171 // depends on the face direction and orientation.
172
173 // The index of the first dof on this face or the cell
174 unsigned int dbase = 0;
175 // The index of the first generalized support point on this face or the cell
176 unsigned int pbase = 0;
177 for (auto f : GeometryInfo<dim>::face_indices())
178 {
179 // Old version with no moments in 2d. See comment below in
180 // initialize_support_points()
181 if (test_values_face.size() == 0)
182 {
183 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
184 nodal_values[dbase + i] =
185 support_point_values[pbase + i]
187 pbase += this->n_dofs_per_face(f);
188 }
189 else
190 {
191 for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
192 {
193 double s = 0.;
194 for (unsigned int k = 0; k < test_values_face.size(); ++k)
195 s +=
196 support_point_values
198 test_values_face[k][i];
199 nodal_values[dbase + i] = s;
200 }
201 pbase += test_values_face.size();
202 }
203 dbase += this->n_dofs_per_face(f);
204 }
205
206 AssertDimension(this->n_unique_faces(), 1);
207 const unsigned int face_no = 0;
208 (void)face_no;
209
210 AssertDimension(dbase,
211 this->n_dofs_per_face(face_no) *
213 AssertDimension(pbase,
214 this->generalized_support_points.size() -
215 test_values_cell.size());
216
217 // Done for BDM1
218 if (dbase == this->n_dofs_per_cell())
219 return;
220
221 // What's missing are the interior
222 // degrees of freedom. In each
223 // point, we take all components of
224 // the solution.
225 Assert((this->n_dofs_per_cell() - dbase) % dim == 0, ExcInternalError());
226
227 for (unsigned int d = 0; d < dim; ++d, dbase += test_values_cell[0].size())
228 {
229 for (unsigned int i = 0; i < test_values_cell[0].size(); ++i)
230 {
231 double s = 0.;
232 for (unsigned int k = 0; k < test_values_cell.size(); ++k)
233 s += support_point_values[pbase + k][d] * test_values_cell[k][i];
234 nodal_values[dbase + i] = s;
235 }
236 }
237
238 Assert(dbase == this->n_dofs_per_cell(), ExcInternalError());
239}
240
241
242
243template <int dim>
244std::vector<unsigned int>
245FE_BDM<dim>::get_dpo_vector(const unsigned int deg)
246{
247 // compute the number of unknowns per cell interior/face/edge
248 //
249 // for the number of interior dofs, this is the number of
250 // polynomials up to degree deg-2 in dim dimensions.
251 //
252 // the element is face-based and we have as many degrees of freedom
253 // on the faces as there are polynomials of degree up to
254 // deg. Observe the odd convention of
255 // PolynomialSpace::n_polynomials()!
256
257 std::vector<unsigned int> dpo(dim + 1, 0u);
258 dpo[dim] =
259 (deg > 1 ? dim * PolynomialSpace<dim>::n_polynomials(deg - 1) : 0u);
260 dpo[dim - 1] = PolynomialSpace<dim - 1>::n_polynomials(deg + 1);
261
262 return dpo;
263}
264
265
266
267template <int dim>
268std::vector<bool>
269FE_BDM<dim>::get_ria_vector(const unsigned int deg)
270{
271 if (dim == 1)
272 {
273 Assert(false, ExcImpossibleInDim(1));
274 return std::vector<bool>();
275 }
276
277 const unsigned int dofs_per_cell = PolynomialsBDM<dim>::n_polynomials(deg);
278 const unsigned int dofs_per_face =
280
281 Assert(GeometryInfo<dim>::faces_per_cell * dofs_per_face <= dofs_per_cell,
283
284 // all dofs need to be
285 // non-additive, since they have
286 // continuity requirements.
287 // however, the interior dofs are
288 // made additive
289 std::vector<bool> ret_val(dofs_per_cell, false);
290 for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
291 i < dofs_per_cell;
292 ++i)
293 ret_val[i] = true;
294
295 return ret_val;
296}
297
298
299namespace internal
300{
301 namespace FE_BDM
302 {
303 namespace
304 {
305 // This function sets up the values of the polynomials we want to
306 // take moments with in the quadrature points. In fact, we multiply
307 // those by the weights, such that the sum of function values and
308 // test_values over quadrature points yields the interpolated degree
309 // of freedom.
310 template <int dim>
311 void
312 initialize_test_values(std::vector<std::vector<double>> &test_values,
313 const Quadrature<dim> & quadrature,
314 const unsigned int deg)
315 {
316 PolynomialsP<dim> poly(deg);
317 std::vector<Tensor<1, dim>> dummy1;
318 std::vector<Tensor<2, dim>> dummy2;
319 std::vector<Tensor<3, dim>> dummy3;
320 std::vector<Tensor<4, dim>> dummy4;
321
322 test_values.resize(quadrature.size());
323
324 for (unsigned int k = 0; k < quadrature.size(); ++k)
325 {
326 test_values[k].resize(poly.n());
327 poly.evaluate(quadrature.point(k),
328 test_values[k],
329 dummy1,
330 dummy2,
331 dummy3,
332 dummy4);
333 for (unsigned int i = 0; i < poly.n(); ++i)
334 {
335 test_values[k][i] *= quadrature.weight(k);
336 }
337 }
338 }
339
340 // This specialization only serves to avoid error messages. Nothing
341 // useful can be computed in dimension zero and thus the vector
342 // length stays zero.
343 template <>
344 void
345 initialize_test_values(std::vector<std::vector<double>> &,
346 const Quadrature<0> &,
347 const unsigned int)
348 {}
349 } // namespace
350 } // namespace FE_BDM
351} // namespace internal
352
353
354template <int dim>
355void
357{
358 // Our support points are quadrature points on faces and inside the
359 // cell. First on the faces, we have to test polynomials of degree
360 // up to deg, which means we need dg+1 points in each direction. The
361 // fact that we do not have tensor product polynomials will be
362 // considered later. In 2d, we can use point values.
363 QGauss<dim - 1> face_points(deg + 1);
364
365 // TODO: the implementation makes the assumption that all faces have the
366 // same number of dofs
367 AssertDimension(this->n_unique_faces(), 1);
368 const unsigned int face_no = 0;
369
370 // Copy the quadrature formula to the face points.
371 this->generalized_face_support_points[face_no].resize(face_points.size());
372 for (unsigned int k = 0; k < face_points.size(); ++k)
373 this->generalized_face_support_points[face_no][k] = face_points.point(k);
374
375 // In the interior, we only test with polynomials of degree up to
376 // deg-2, thus we use deg points. Note that deg>=1 and the lowest
377 // order element has no points in the cell, such that we have to
378 // distinguish this case.
379 QGauss<dim> cell_points(deg == 1 ? 0 : deg);
380
381 // Compute the size of the whole support point set
382 const unsigned int npoints =
383 cell_points.size() + GeometryInfo<dim>::faces_per_cell * face_points.size();
384
385 this->generalized_support_points.resize(npoints);
386
387 Quadrature<dim> faces =
388 QProjector<dim>::project_to_all_faces(this->reference_cell(), face_points);
389 for (unsigned int k = 0;
390 k < face_points.size() * GeometryInfo<dim>::faces_per_cell;
391 ++k)
392 this->generalized_support_points[k] = faces.point(
393 k +
394 QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
395 0,
396 true,
397 false,
398 false,
399 this->n_dofs_per_face(face_no)));
400
401 // Currently, for backward compatibility, we do not use moments, but
402 // point values on faces in 2d. In 3d, this is impossible, since the
403 // moments are only taken with respect to PolynomialsP.
404 if (dim > 2)
405 internal::FE_BDM::initialize_test_values(test_values_face,
406 face_points,
407 deg);
408
409 if (deg <= 1)
410 return;
411
412 // Remember where interior points start
413 const unsigned int ibase =
414 face_points.size() * GeometryInfo<dim>::faces_per_cell;
415 for (unsigned int k = 0; k < cell_points.size(); ++k)
416 {
417 this->generalized_support_points[ibase + k] = cell_points.point(k);
418 }
419 // Finally, compute the values of
420 // the test functions in the
421 // interior quadrature points
422
423 internal::FE_BDM::initialize_test_values(test_values_cell,
424 cell_points,
425 deg - 2);
426}
427
428
429
430/*-------------- Explicit Instantiations -------------------------------*/
431#include "fe_bdm.inst"
432
void initialize_support_points(const unsigned int bdm_degree)
Definition fe_bdm.cc:356
FE_BDM(const unsigned int p)
Definition fe_bdm.cc:44
static std::vector< bool > get_ria_vector(const unsigned int degree)
Definition fe_bdm.cc:269
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
Definition fe_bdm.cc:150
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
Definition fe_bdm.cc:159
void initialize_quad_dof_index_permutation_and_sign_change()
Definition fe_bdm.cc:115
virtual std::string get_name() const override
Definition fe_bdm.cc:128
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition fe_bdm.cc:245
FullMatrix< double > inverse_node_matrix
std::vector< MappingKind > mapping_kind
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_unique_faces() const
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
FullMatrix< double > interface_constraints
Definition fe.h:2428
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2416
size_type n() const
void invert(const FullMatrix< number2 > &M)
size_type m() const
static unsigned int n_polynomials(const unsigned int n)
void evaluate(const Point< dim > &unit_point, std::vector< double > &values, std::vector< Tensor< 1, dim > > &grads, std::vector< Tensor< 2, dim > > &grad_grads, std::vector< Tensor< 3, dim > > &third_derivatives, std::vector< Tensor< 4, dim > > &fourth_derivatives) const override
static unsigned int n_polynomials(const unsigned int degree)
static Quadrature< dim > project_to_all_faces(const ReferenceCell &reference_cell, const hp::QCollection< dim - 1 > &quadrature)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
@ mapping_bdm
Definition mapping.h:138
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
void compute_face_embedding_matrices(const FiniteElement< dim, spacedim > &fe, FullMatrix< number >(&matrices)[GeometryInfo< dim >::max_children_per_face], const unsigned int face_coarse, const unsigned int face_fine, const double threshold=1.e-12)
FullMatrix< double > compute_node_matrix(const FiniteElement< dim, spacedim > &fe)
STL namespace.
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()