Reference documentation for deal.II version 9.5.0
\(\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_q_bubbles.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2012 - 2023 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>
27#include <deal.II/fe/fe_tools.h>
29
31#include <deal.II/grid/tria.h>
32
33#include <memory>
34#include <sstream>
35#include <vector>
36
38
39
40namespace internal
41{
42 namespace FE_Q_Bubbles
43 {
44 namespace
45 {
46 template <int dim, int spacedim>
47 inline void
48 compute_embedding_matrices(
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
57 std::unique_ptr<Quadrature<dim>> q_fine;
58 Quadrature<1> q_dummy(std::vector<Point<1>>(1),
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
117 tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
119
121 dh.distribute_dofs(fe);
122
124 fe,
125 *q_fine,
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;
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);
156 Point<dim> quad_tmp;
157 for (unsigned int k = 0; k < dim; ++k)
158 quad_tmp(k) = fine.quadrature_point(q)(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
187template <int dim, int spacedim>
189 : FE_Q_Base<dim, spacedim>(TensorProductPolynomialsBubbles<dim>(
190 Polynomials::generate_complete_Lagrange_basis(
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
206 Point<dim> point;
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);
212
214 if (dim == spacedim)
215 {
216 internal::FE_Q_Bubbles::compute_embedding_matrices(*this,
217 this->prolongation,
218 false);
219 // Fill restriction matrices with L2-projection
221 }
222}
223
224
225
226template <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
245 Point<dim> point;
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);
251
253 if (dim == spacedim)
254 {
255 internal::FE_Q_Bubbles::compute_embedding_matrices(*this,
256 this->prolongation,
257 false);
258 // Fill restriction matrices with L2-projection
260 }
261}
262
263
264
265template <int dim, int spacedim>
266std::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),
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
344template <int dim, int spacedim>
345std::unique_ptr<FiniteElement<dim, spacedim>>
347{
348 return std::make_unique<FE_Q_Bubbles<dim, spacedim>>(*this);
349}
350
351
352
353template <int dim, int spacedim>
354void
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
383template <int dim, int spacedim>
384void
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
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
420template <int dim, int spacedim>
421std::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
431template <int dim, int spacedim>
432std::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
447template <int dim, int spacedim>
448bool
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
463template <int dim, int spacedim>
464const FullMatrix<double> &
466 const unsigned int child,
467 const RefinementCase<dim> &refinement_case) const
468{
469 AssertIndexRange(refinement_case,
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
484template <int dim, int spacedim>
485const FullMatrix<double> &
487 const unsigned int child,
488 const RefinementCase<dim> &refinement_case) const
489{
490 AssertIndexRange(refinement_case,
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
506template <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 Point< spacedim > & quadrature_point(const unsigned int q_point) const
double JxW(const unsigned int q_point) const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
const unsigned int q_degree
Definition fe_q_base.h:346
void initialize(const std::vector< Point< 1 > > &support_points_1d)
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
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
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
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)
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
Definition point.h:112
const Point< dim > & point(const unsigned int i) const
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
virtual void execute_coarsening_and_refinement()
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:285
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:556
static unsigned int n_children(const RefinementCase< dim > &refinement_case)