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