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