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