Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20:02+00:00
fe_q_bubbles.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2012 - 2021 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
23
24 #include <deal.II/fe/fe_dgq.h>
25 #include <deal.II/fe/fe_nothing.h>
27 #include <deal.II/fe/fe_tools.h>
28 #include <deal.II/fe/fe_values.h>
29 #include <deal.II/fe/mapping_q1.h>
30
32 #include <deal.II/grid/tria.h>
33
34 #include <memory>
35 #include <sstream>
36 #include <vector>
37
39
40
41 namespace internal
42 {
43  namespace FE_Q_Bubbles
44  {
45  namespace
46  {
47  template <int dim, int spacedim>
48  inline void
50  const ::FE_Q_Bubbles<dim, spacedim> & fe,
51  std::vector<std::vector<FullMatrix<double>>> &matrices,
52  const bool isotropic_only)
53  {
54  const unsigned int dpc = fe.n_dofs_per_cell();
55  const unsigned int degree = fe.degree;
56
57  // Initialize quadrature formula on fine cells
60  std::vector<double>(1, 1.));
61  switch (dim)
62  {
63  case 1:
64  if (spacedim == 1)
65  q_fine = std::make_unique<QGauss<dim>>(degree + 1);
66  else if (spacedim == 2)
67  q_fine =
68  std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
69  q_dummy);
70  else
71  q_fine =
72  std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
73  q_dummy,
74  q_dummy);
75  break;
76  case 2:
77  if (spacedim == 2)
78  q_fine = std::make_unique<QGauss<dim>>(degree + 1);
79  else
80  q_fine =
81  std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
82  QGauss<1>(degree + 1),
83  q_dummy);
84  break;
85  case 3:
86  q_fine = std::make_unique<QGauss<dim>>(degree + 1);
87  break;
88  default:
89  Assert(false, ExcInternalError());
90  }
91
92  Assert(q_fine.get() != nullptr, ExcInternalError());
93  const unsigned int nq = q_fine->size();
94
95  // loop over all possible refinement cases
96  unsigned int ref_case = (isotropic_only) ?
99  for (; ref_case <= RefinementCase<dim>::isotropic_refinement;
100  ++ref_case)
101  {
102  const unsigned int nc =
104
105  for (unsigned int i = 0; i < nc; ++i)
106  {
107  Assert(matrices[ref_case - 1][i].n() == dpc,
108  ExcDimensionMismatch(matrices[ref_case - 1][i].n(),
109  dpc));
110  Assert(matrices[ref_case - 1][i].m() == dpc,
111  ExcDimensionMismatch(matrices[ref_case - 1][i].m(),
112  dpc));
113  }
114
115  // create a respective refinement on the triangulation
117  GridGenerator::hyper_cube(tr, 0, 1);
118  tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
120
122  dh.distribute_dofs(fe);
123
125  fe,
126  *q_fine,
129  update_values);
130
131  const unsigned int n_dofs = dh.n_dofs();
132
133  FullMatrix<double> fine_mass(n_dofs);
134  FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
135
136  std::vector<std::vector<types::global_dof_index>> child_ldi(
137  nc, std::vector<types::global_dof_index>(fe.n_dofs_per_cell()));
138
139  // now create the mass matrix and all the right_hand sides
140  unsigned int child_no = 0;
141  typename ::DoFHandler<dim>::active_cell_iterator cell =
142  dh.begin_active();
143  for (; cell != dh.end(); ++cell, ++child_no)
144  {
145  fine.reinit(cell);
146  cell->get_dof_indices(child_ldi[child_no]);
147
148  for (unsigned int q = 0; q < nq; ++q)
149  for (unsigned int i = 0; i < dpc; ++i)
150  for (unsigned int j = 0; j < dpc; ++j)
151  {
152  const unsigned int gdi = child_ldi[child_no][i];
153  const unsigned int gdj = child_ldi[child_no][j];
154  fine_mass(gdi, gdj) += fine.shape_value(i, q) *
155  fine.shape_value(j, q) *
156  fine.JxW(q);
158  for (unsigned int k = 0; k < dim; ++k)
160  coarse_rhs_matrix(gdi, j) +=
161  fine.shape_value(i, q) * fe.shape_value(j, quad_tmp) *
162  fine.JxW(q);
163  }
164  }
165
166  // now solve for all right-hand sides simultaneously
167  ::FullMatrix<double> solution(n_dofs, dpc);
168  fine_mass.gauss_jordan();
169  fine_mass.mmult(solution, coarse_rhs_matrix);
170
171  // and distribute to the fine cell matrices
172  for (unsigned int child_no = 0; child_no < nc; ++child_no)
173  for (unsigned int i = 0; i < dpc; ++i)
174  for (unsigned int j = 0; j < dpc; ++j)
175  {
176  const unsigned int gdi = child_ldi[child_no][i];
177  // remove small entries
178  if (std::fabs(solution(gdi, j)) > 1.e-12)
179  matrices[ref_case - 1][child_no](i, j) = solution(gdi, j);
180  }
181  }
182  }
183  } // namespace
184  } // namespace FE_Q_Bubbles
185 } // namespace internal
186
187
188 template <int dim, int spacedim>
189 FE_Q_Bubbles<dim, spacedim>::FE_Q_Bubbles(const unsigned int q_degree)
190  : FE_Q_Base<dim, spacedim>(TensorProductPolynomialsBubbles<dim>(
192  QGaussLobatto<1>(q_degree + 1).get_points())),
193  FiniteElementData<dim>(get_dpo_vector(q_degree),
194  1,
195  q_degree + 1,
196  FiniteElementData<dim>::H1),
197  get_riaf_vector(q_degree))
198  , n_bubbles((q_degree <= 1) ? 1 : dim)
199 {
200  Assert(q_degree > 0,
201  ExcMessage("This element can only be used for polynomial degrees "
202  "greater than zero"));
203
204  this->initialize(QGaussLobatto<1>(q_degree + 1).get_points());
205
206  // adjust unit support point for discontinuous node
208  for (unsigned int d = 0; d < dim; ++d)
209  point[d] = 0.5;
210  for (unsigned int i = 0; i < n_bubbles; ++i)
211  this->unit_support_points.push_back(point);
212  AssertDimension(this->n_dofs_per_cell(), this->unit_support_points.size());
213
215  if (dim == spacedim)
216  {
218  this->prolongation,
219  false);
220  // Fill restriction matrices with L2-projection
222  }
223 }
224
225
226
227 template <int dim, int spacedim>
229  : FE_Q_Base<dim, spacedim>(
231  Polynomials::generate_complete_Lagrange_basis(points.get_points())),
232  FiniteElementData<dim>(get_dpo_vector(points.size() - 1),
233  1,
234  points.size(),
235  FiniteElementData<dim>::H1),
236  get_riaf_vector(points.size() - 1))
237  , n_bubbles((points.size() - 1 <= 1) ? 1 : dim)
238 {
239  Assert(points.size() > 1,
240  ExcMessage("This element can only be used for polynomial degrees "
241  "at least one"));
242
243  this->initialize(points.get_points());
244
245  // adjust unit support point for discontinuous node
247  for (unsigned int d = 0; d < dim; ++d)
248  point[d] = 0.5;
249  for (unsigned int i = 0; i < n_bubbles; ++i)
250  this->unit_support_points.push_back(point);
251  AssertDimension(this->n_dofs_per_cell(), this->unit_support_points.size());
252
254  if (dim == spacedim)
255  {
257  this->prolongation,
258  false);
259  // Fill restriction matrices with L2-projection
261  }
262 }
263
264
265
266 template <int dim, int spacedim>
267 std::string
269 {
270  // note that the FETools::get_fe_by_name function depends on the
271  // particular format of the string this function returns, so they have to be
272  // kept in synch
273
274  std::ostringstream namebuf;
275  bool type = true;
276  const unsigned int n_points = this->degree;
277  std::vector<double> points(n_points);
278  const unsigned int dofs_per_cell = this->n_dofs_per_cell();
279  const std::vector<Point<dim>> &unit_support_points =
280  this->unit_support_points;
281  unsigned int index = 0;
282
283  // Decode the support points in one coordinate direction.
284  for (unsigned int j = 0; j < dofs_per_cell; ++j)
285  {
286  if ((dim > 1) ? (unit_support_points[j](1) == 0 &&
287  ((dim > 2) ? unit_support_points[j](2) == 0 : true)) :
288  true)
289  {
290  if (index == 0)
291  points[index] = unit_support_points[j](0);
292  else if (index == 1)
293  points[n_points - 1] = unit_support_points[j](0);
294  else
295  points[index - 1] = unit_support_points[j](0);
296
297  index++;
298  }
299  }
300  // Do not consider the discontinuous node for dimension 1
301  Assert(index == n_points || (dim == 1 && index == n_points + n_bubbles),
302  ExcMessage(
303  "Could not decode support points in one coordinate direction."));
304
305  // Check whether the support points are equidistant.
306  for (unsigned int j = 0; j < n_points; ++j)
307  if (std::fabs(points[j] - static_cast<double>(j) / (this->degree - 1)) >
308  1e-15)
309  {
310  type = false;
311  break;
312  }
313
314  if (type == true)
315  {
316  if (this->degree > 3)
317  namebuf << "FE_Q_Bubbles<" << Utilities::dim_string(dim, spacedim)
318  << ">(QIterated(QTrapezoid()," << this->degree - 1 << "))";
319  else
320  namebuf << "FE_Q_Bubbles<" << Utilities::dim_string(dim, spacedim)
321  << ">(" << this->degree - 1 << ")";
322  }
323  else
324  {
325  // Check whether the support points come from QGaussLobatto.
326  const QGaussLobatto<1> points_gl(n_points);
327  type = true;
328  for (unsigned int j = 0; j < n_points; ++j)
329  if (points[j] != points_gl.point(j)(0))
330  {
331  type = false;
332  break;
333  }
334  if (type == true)
335  namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree - 1 << ")";
336  else
337  namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree
338  << "))";
339  }
340  return namebuf.str();
341 }
342
343
344
345 template <int dim, int spacedim>
346 std::unique_ptr<FiniteElement<dim, spacedim>>
348 {
349  return std::make_unique<FE_Q_Bubbles<dim, spacedim>>(*this);
350 }
351
352
353
354 template <int dim, int spacedim>
355 void
358  const std::vector<Vector<double>> &support_point_values,
359  std::vector<double> & nodal_values) const
360 {
361  Assert(support_point_values.size() == this->unit_support_points.size(),
362  ExcDimensionMismatch(support_point_values.size(),
363  this->unit_support_points.size()));
364  Assert(nodal_values.size() == this->n_dofs_per_cell(),
365  ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
366  Assert(support_point_values[0].size() == this->n_components(),
367  ExcDimensionMismatch(support_point_values[0].size(),
368  this->n_components()));
369
370  for (unsigned int i = 0; i < this->n_dofs_per_cell() - 1; ++i)
371  {
372  const std::pair<unsigned int, unsigned int> index =
373  this->system_to_component_index(i);
374  nodal_values[i] = support_point_values[i](index.first);
375  }
376
377  // We don't use the bubble functions for local interpolation
378  for (unsigned int i = 0; i < n_bubbles; ++i)
379  nodal_values[nodal_values.size() - i - 1] = 0.;
380 }
381
382
383
384 template <int dim, int spacedim>
385 void
387  const FiniteElement<dim, spacedim> &x_source_fe,
388  FullMatrix<double> & interpolation_matrix) const
389 {
390  // We don't know how to do this properly, yet.
391  // However, for SolutionTransfer to work we need to provide an implementation
392  // for the case that the x_source_fe is identical to this FE
393  using FEQBUBBLES = FE_Q_Bubbles<dim, spacedim>;
394
395  AssertThrow(
396  (x_source_fe.get_name().find("FE_Q_Bubbles<") == 0) ||
397  (dynamic_cast<const FEQBUBBLES *>(&x_source_fe) != nullptr),
399  Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
400  ExcDimensionMismatch(interpolation_matrix.m(),
401  this->n_dofs_per_cell()));
402  Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
403  ExcDimensionMismatch(interpolation_matrix.m(),
404  x_source_fe.n_dofs_per_cell()));
405
406  // Provide a short cut in case we are just inquiring the identity
407  auto casted_fe = dynamic_cast<const FEQBUBBLES *>(&x_source_fe);
408  if (casted_fe != nullptr && casted_fe->degree == this->degree)
409  for (unsigned int i = 0; i < interpolation_matrix.m(); ++i)
410  interpolation_matrix.set(i, i, 1.);
411  // else we need to do more...
412  else
413  Assert(
414  false,
415  (typename FiniteElement<dim,
416  spacedim>::ExcInterpolationNotImplemented()));
417 }
418
419
420
421 template <int dim, int spacedim>
422 std::vector<bool>
424 {
425  const unsigned int n_cont_dofs = Utilities::fixed_power<dim>(q_deg + 1);
426  const unsigned int n_bubbles = (q_deg <= 1 ? 1 : dim);
427  return std::vector<bool>(n_cont_dofs + n_bubbles, true);
428 }
429
430
431
432 template <int dim, int spacedim>
433 std::vector<unsigned int>
435 {
436  std::vector<unsigned int> dpo(dim + 1, 1U);
437  for (unsigned int i = 1; i < dpo.size(); ++i)
438  dpo[i] = dpo[i - 1] * (q_deg - 1);
439
440  // Then add the bubble functions; they are all associated with the
441  // cell interior
442  dpo[dim] += (q_deg <= 1 ? 1 : dim);
443  return dpo;
444 }
445
446
447
448 template <int dim, int spacedim>
449 bool
451  const unsigned int shape_index,
452  const unsigned int face_index) const
453 {
454  // discontinuous functions have no support on faces
455  if (shape_index >= this->n_dofs_per_cell() - n_bubbles)
456  return false;
457  else
459  face_index);
460 }
461
462
463
464 template <int dim, int spacedim>
465 const FullMatrix<double> &
467  const unsigned int child,
468  const RefinementCase<dim> &refinement_case) const
469 {
470  AssertIndexRange(refinement_case,
472  Assert(refinement_case != RefinementCase<dim>::no_refinement,
473  ExcMessage(
474  "Prolongation matrices are only available for refined cells!"));
475  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
476
477  Assert(this->prolongation[refinement_case - 1][child].n() != 0,
478  ExcMessage("This prolongation matrix has not been computed yet!"));
479  // finally return the matrix
480  return this->prolongation[refinement_case - 1][child];
481 }
482
483
484
485 template <int dim, int spacedim>
486 const FullMatrix<double> &
488  const unsigned int child,
489  const RefinementCase<dim> &refinement_case) const
490 {
491  AssertIndexRange(refinement_case,
493  Assert(refinement_case != RefinementCase<dim>::no_refinement,
494  ExcMessage(
495  "Restriction matrices are only available for refined cells!"));
496  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
497
498  Assert(this->restriction[refinement_case - 1][child].n() != 0,
499  ExcMessage("This restriction matrix has not been computed yet!"));
500
501  // finally return the matrix
502  return this->restriction[refinement_case - 1][child];
503 }
504
505
506
507 template <int dim, int spacedim>
510  const FiniteElement<dim, spacedim> &fe_other,
511  const unsigned int codim) const
512 {
513  Assert(codim <= dim, ExcImpossibleInDim(dim));
514
515  // vertex/line/face domination
516  // (if fe_other is derived from FE_DGQ)
517  // ------------------------------------
518  if (codim > 0)
519  if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
520  // there are no requirements between continuous and discontinuous elements
522
523  // vertex/line/face domination
524  // (if fe_other is not derived from FE_DGQ)
525  // & cell domination
526  // ----------------------------------------
527  if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
528  dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
529  {
530  if (this->degree < fe_bubbles_other->degree)
532  else if (this->degree == fe_bubbles_other->degree)
534  else
536  }
537  else if (const FE_Nothing<dim> *fe_nothing =
538  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
539  {
540  if (fe_nothing->is_dominating())
542  else
543  // the FE_Nothing has no degrees of freedom and it is typically used
544  // in a context where we don't require any continuity along the
545  // interface
547  }
548
549  Assert(false, ExcNotImplemented());
551 }
552
553
554 // explicit instantiations
555 #include "fe_q_bubbles.inst"
556
cell_iterator end() const
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
active_cell_iterator begin_active(const unsigned int level=0) const
types::global_dof_index n_dofs() const
const double & shape_value(const unsigned int function_no, const unsigned int point_no) const
const Point< spacedim > & quadrature_point(const unsigned int q) const
double JxW(const unsigned int quadrature_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell)
Definition: fe_dgq.h:113
void initialize(const std::vector< Point< 1 >> &support_points_1d)
Definition: fe_q_base.cc:436
const unsigned int q_degree
Definition: fe_q_base.h:346
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_q_base.cc:1590
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
static std::vector< bool > get_riaf_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) const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case) const override
const unsigned int n_bubbles
Definition: fe_q_bubbles.h:195
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
FE_Q_Bubbles(const unsigned int p)
virtual std::string get_name() const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
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
unsigned int n_dofs_per_cell() const
virtual std::string get_name() const =0
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2386
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2424
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2400
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void set(const size_type i, const size_type j, const number value)
size_type n() const
void gauss_jordan()
size_type m() const
const std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
virtual void execute_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
@ update_values
Shape function values.
@ update_JxW_values
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition: mapping.cc:260
Expression fabs(const Expression &x)
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_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:702
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:558
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)
static unsigned int n_children(const RefinementCase< dim > &refinement_case)