Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2628-g5436af8a2d 2025-02-12 20:00:00+00:00
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_q_bubbles.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
19
22
23#include <deal.II/fe/fe_dgq.h>
26#include <deal.II/fe/fe_tools.h>
28
30#include <deal.II/grid/tria.h>
31
32#include <memory>
33#include <sstream>
34#include <vector>
35
37
38
39namespace internal
40{
41 namespace FE_Q_Bubbles
42 {
43 namespace
44 {
45 template <int dim, int spacedim>
46 inline void
47 compute_embedding_matrices(
48 const ::FE_Q_Bubbles<dim, spacedim> &fe,
49 std::vector<std::vector<FullMatrix<double>>> &matrices,
50 const bool isotropic_only)
51 {
52 const unsigned int dpc = fe.n_dofs_per_cell();
53 const unsigned int degree = fe.degree;
54
55 // Initialize quadrature formula on fine cells
56 std::unique_ptr<Quadrature<dim>> q_fine;
57 Quadrature<1> q_dummy(std::vector<Point<1>>(1),
58 std::vector<double>(1, 1.));
59 switch (dim)
60 {
61 case 1:
62 if (spacedim == 1)
63 q_fine = std::make_unique<QGauss<dim>>(degree + 1);
64 else if (spacedim == 2)
65 q_fine =
66 std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
67 q_dummy);
68 else
69 q_fine =
70 std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
71 q_dummy,
72 q_dummy);
73 break;
74 case 2:
75 if (spacedim == 2)
76 q_fine = std::make_unique<QGauss<dim>>(degree + 1);
77 else
78 q_fine =
79 std::make_unique<QAnisotropic<dim>>(QGauss<1>(degree + 1),
80 QGauss<1>(degree + 1),
81 q_dummy);
82 break;
83 case 3:
84 q_fine = std::make_unique<QGauss<dim>>(degree + 1);
85 break;
86 default:
88 }
89
90 Assert(q_fine.get() != nullptr, ExcInternalError());
91 const unsigned int nq = q_fine->size();
92
93 // loop over all possible refinement cases
94 for (unsigned int ref_case =
97 ref_case <= RefinementCase<dim>::isotropic_refinement;
98 ++ref_case)
99 {
100 const unsigned int nc =
101 fe.reference_cell().template n_children<dim>(
102 RefinementCase<dim>(ref_case));
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<" << Utilities::dim_string(dim, spacedim)
335 << ">(" << this->degree - 1 << ")";
336 else
337 namebuf << "FE_Q_Bubbles<" << Utilities::dim_string(dim, spacedim)
338 << ">(QUnknownNodes(" << this->degree << "))";
339 }
340 return namebuf.str();
341}
342
343
344
345template <int dim, int spacedim>
346std::unique_ptr<FiniteElement<dim, spacedim>>
348{
349 return std::make_unique<FE_Q_Bubbles<dim, spacedim>>(*this);
350}
351
352
353
354template <int dim, int spacedim>
355void
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
384template <int dim, int spacedim>
385void
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
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
421template <int dim, int spacedim>
422std::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
432template <int dim, int spacedim>
433std::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
448template <int dim, int spacedim>
449bool
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
464template <int dim, int spacedim>
465const FullMatrix<double> &
467 const unsigned int child,
468 const RefinementCase<dim> &refinement_case) const
469{
470 AssertIndexRange(refinement_case,
474 "Prolongation matrices are only available for refined cells!"));
476 child, this->reference_cell().template n_children<dim>(refinement_case));
477
478 Assert(this->prolongation[refinement_case - 1][child].n() != 0,
479 ExcMessage("This prolongation matrix has not been computed yet!"));
480 // finally return the matrix
481 return this->prolongation[refinement_case - 1][child];
482}
483
484
485
486template <int dim, int spacedim>
487const FullMatrix<double> &
489 const unsigned int child,
490 const RefinementCase<dim> &refinement_case) const
491{
492 AssertIndexRange(refinement_case,
496 "Restriction matrices are only available for refined cells!"));
498 child, this->reference_cell().template n_children<dim>(refinement_case));
499
500 Assert(this->restriction[refinement_case - 1][child].n() != 0,
501 ExcMessage("This restriction matrix has not been computed yet!"));
502
503 // finally return the matrix
504 return this->restriction[refinement_case - 1][child];
505}
506
507
508
509template <int dim, int spacedim>
512 const FiniteElement<dim, spacedim> &fe_other,
513 const unsigned int codim) const
514{
515 Assert(codim <= dim, ExcImpossibleInDim(dim));
516
517 // vertex/line/face domination
518 // (if fe_other is derived from FE_DGQ)
519 // ------------------------------------
520 if (codim > 0)
521 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
522 // there are no requirements between continuous and discontinuous elements
524
525 // vertex/line/face domination
526 // (if fe_other is not derived from FE_DGQ)
527 // & cell domination
528 // ----------------------------------------
529 if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
530 dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
531 {
532 if (this->degree < fe_bubbles_other->degree)
534 else if (this->degree == fe_bubbles_other->degree)
536 else
538 }
539 else if (const FE_Nothing<dim> *fe_nothing =
540 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
541 {
542 if (fe_nothing->is_dominating())
544 else
545 // the FE_Nothing has no degrees of freedom and it is typically used
546 // in a context where we don't require any continuity along the
547 // interface
549 }
550
553}
554
555
556// explicit instantiations
557#include "fe/fe_q_bubbles.inst"
558
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:340
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:2523
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:2561
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2537
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:113
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:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#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.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:316
std::size_t size
Definition mpi.cc:739
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:555