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