Reference documentation for deal.II version GIT relicensing-233-g802318d791 2024-03-28 20:20:02+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\}}\)
Loading...
Searching...
No Matches
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 =
102
103 for (unsigned int i = 0; i < nc; ++i)
104 {
105 Assert(matrices[ref_case - 1][i].n() == dpc,
106 ExcDimensionMismatch(matrices[ref_case - 1][i].n(),
107 dpc));
108 Assert(matrices[ref_case - 1][i].m() == dpc,
109 ExcDimensionMismatch(matrices[ref_case - 1][i].m(),
110 dpc));
111 }
112
113 // create a respective refinement on the triangulation
116 tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
118
120 dh.distribute_dofs(fe);
121
123 fe,
124 *q_fine,
128
129 const unsigned int n_dofs = dh.n_dofs();
130
131 FullMatrix<double> fine_mass(n_dofs);
132 FullMatrix<double> coarse_rhs_matrix(n_dofs, dpc);
133
134 std::vector<std::vector<types::global_dof_index>> child_ldi(
135 nc, std::vector<types::global_dof_index>(fe.n_dofs_per_cell()));
136
137 // now create the mass matrix and all the right_hand sides
138 unsigned int child_no = 0;
140 dh.begin_active();
141 for (; cell != dh.end(); ++cell, ++child_no)
142 {
143 fine.reinit(cell);
144 cell->get_dof_indices(child_ldi[child_no]);
145
146 for (unsigned int q = 0; q < nq; ++q)
147 for (unsigned int i = 0; i < dpc; ++i)
148 for (unsigned int j = 0; j < dpc; ++j)
149 {
150 const unsigned int gdi = child_ldi[child_no][i];
151 const unsigned int gdj = child_ldi[child_no][j];
152 fine_mass(gdi, gdj) += fine.shape_value(i, q) *
153 fine.shape_value(j, q) *
154 fine.JxW(q);
155 Point<dim> quad_tmp;
156 for (unsigned int k = 0; k < dim; ++k)
157 quad_tmp[k] = fine.quadrature_point(q)[k];
158 coarse_rhs_matrix(gdi, j) +=
159 fine.shape_value(i, q) * fe.shape_value(j, quad_tmp) *
160 fine.JxW(q);
161 }
162 }
163
164 // now solve for all right-hand sides simultaneously
165 ::FullMatrix<double> solution(n_dofs, dpc);
166 fine_mass.gauss_jordan();
167 fine_mass.mmult(solution, coarse_rhs_matrix);
168
169 // and distribute to the fine cell matrices
170 for (unsigned int child_no = 0; child_no < nc; ++child_no)
171 for (unsigned int i = 0; i < dpc; ++i)
172 for (unsigned int j = 0; j < dpc; ++j)
173 {
174 const unsigned int gdi = child_ldi[child_no][i];
175 // remove small entries
176 if (std::fabs(solution(gdi, j)) > 1.e-12)
177 matrices[ref_case - 1][child_no](i, j) = solution(gdi, j);
178 }
179 }
180 }
181 } // namespace
182 } // namespace FE_Q_Bubbles
183} // namespace internal
184
185
186template <int dim, int spacedim>
188 : FE_Q_Base<dim, spacedim>(TensorProductPolynomialsBubbles<dim>(
189 Polynomials::generate_complete_Lagrange_basis(
190 QGaussLobatto<1>(q_degree + 1).get_points())),
191 FiniteElementData<dim>(get_dpo_vector(q_degree),
192 1,
193 q_degree + 1,
194 FiniteElementData<dim>::H1),
195 get_riaf_vector(q_degree))
196 , n_bubbles((q_degree <= 1) ? 1 : dim)
197{
198 Assert(q_degree > 0,
199 ExcMessage("This element can only be used for polynomial degrees "
200 "greater than zero"));
201
202 this->initialize(QGaussLobatto<1>(q_degree + 1).get_points());
203
204 // adjust unit support point for discontinuous node
205 Point<dim> point;
206 for (unsigned int d = 0; d < dim; ++d)
207 point[d] = 0.5;
208 for (unsigned int i = 0; i < n_bubbles; ++i)
209 this->unit_support_points.push_back(point);
211
213 if (dim == spacedim)
214 {
215 internal::FE_Q_Bubbles::compute_embedding_matrices(*this,
216 this->prolongation,
217 false);
218 // Fill restriction matrices with L2-projection
220 }
221}
222
223
224
225template <int dim, int spacedim>
227 : FE_Q_Base<dim, spacedim>(
229 Polynomials::generate_complete_Lagrange_basis(points.get_points())),
230 FiniteElementData<dim>(get_dpo_vector(points.size() - 1),
231 1,
232 points.size(),
233 FiniteElementData<dim>::H1),
234 get_riaf_vector(points.size() - 1))
235 , n_bubbles((points.size() - 1 <= 1) ? 1 : dim)
236{
237 Assert(points.size() > 1,
238 ExcMessage("This element can only be used for polynomial degrees "
239 "at least one"));
240
241 this->initialize(points.get_points());
242
243 // adjust unit support point for discontinuous node
244 Point<dim> point;
245 for (unsigned int d = 0; d < dim; ++d)
246 point[d] = 0.5;
247 for (unsigned int i = 0; i < n_bubbles; ++i)
248 this->unit_support_points.push_back(point);
250
252 if (dim == spacedim)
253 {
254 internal::FE_Q_Bubbles::compute_embedding_matrices(*this,
255 this->prolongation,
256 false);
257 // Fill restriction matrices with L2-projection
259 }
260}
261
262
263
264template <int dim, int spacedim>
265std::string
267{
268 // note that the FETools::get_fe_by_name function depends on the
269 // particular format of the string this function returns, so they have to be
270 // kept in synch
271
272 std::ostringstream namebuf;
273 bool type = true;
274 const unsigned int n_points = this->degree;
275 std::vector<double> points(n_points);
276 const unsigned int dofs_per_cell = this->n_dofs_per_cell();
277 const std::vector<Point<dim>> &unit_support_points =
278 this->unit_support_points;
279 unsigned int index = 0;
280
281 // Decode the support points in one coordinate direction.
282 for (unsigned int j = 0; j < dofs_per_cell; ++j)
283 {
284 if ((dim > 1) ? (unit_support_points[j][1] == 0 &&
285 ((dim > 2) ? unit_support_points[j][2] == 0 : true)) :
286 true)
287 {
288 if (index == 0)
289 points[index] = unit_support_points[j][0];
290 else if (index == 1)
291 points[n_points - 1] = unit_support_points[j][0];
292 else
293 points[index - 1] = unit_support_points[j][0];
294
295 ++index;
296 }
297 }
298 // Do not consider the discontinuous node for dimension 1
299 Assert(index == n_points || (dim == 1 && index == n_points + n_bubbles),
301 "Could not decode support points in one coordinate direction."));
302
303 // Check whether the support points are equidistant.
304 for (unsigned int j = 0; j < n_points; ++j)
305 if (std::fabs(points[j] - static_cast<double>(j) / (this->degree - 1)) >
306 1e-15)
307 {
308 type = false;
309 break;
310 }
311
312 if (type == true)
313 {
314 if (this->degree > 3)
315 namebuf << "FE_Q_Bubbles<" << Utilities::dim_string(dim, spacedim)
316 << ">(QIterated(QTrapezoid()," << this->degree - 1 << "))";
317 else
318 namebuf << "FE_Q_Bubbles<" << Utilities::dim_string(dim, spacedim)
319 << ">(" << this->degree - 1 << ")";
320 }
321 else
322 {
323 // Check whether the support points come from QGaussLobatto.
324 const QGaussLobatto<1> points_gl(n_points);
325 type = true;
326 for (unsigned int j = 0; j < n_points; ++j)
327 if (points[j] != points_gl.point(j)[0])
328 {
329 type = false;
330 break;
331 }
332 if (type == true)
333 namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree - 1 << ")";
334 else
335 namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree
336 << "))";
337 }
338 return namebuf.str();
339}
340
341
342
343template <int dim, int spacedim>
344std::unique_ptr<FiniteElement<dim, spacedim>>
346{
347 return std::make_unique<FE_Q_Bubbles<dim, spacedim>>(*this);
348}
349
350
351
352template <int dim, int spacedim>
353void
356 const std::vector<Vector<double>> &support_point_values,
357 std::vector<double> &nodal_values) const
358{
359 Assert(support_point_values.size() == this->unit_support_points.size(),
360 ExcDimensionMismatch(support_point_values.size(),
361 this->unit_support_points.size()));
362 Assert(nodal_values.size() == this->n_dofs_per_cell(),
363 ExcDimensionMismatch(nodal_values.size(), this->n_dofs_per_cell()));
364 Assert(support_point_values[0].size() == this->n_components(),
365 ExcDimensionMismatch(support_point_values[0].size(),
366 this->n_components()));
367
368 for (unsigned int i = 0; i < this->n_dofs_per_cell() - 1; ++i)
369 {
370 const std::pair<unsigned int, unsigned int> index =
371 this->system_to_component_index(i);
372 nodal_values[i] = support_point_values[i](index.first);
373 }
374
375 // We don't use the bubble functions for local interpolation
376 for (unsigned int i = 0; i < n_bubbles; ++i)
377 nodal_values[nodal_values.size() - i - 1] = 0.;
378}
379
380
381
382template <int dim, int spacedim>
383void
385 const FiniteElement<dim, spacedim> &x_source_fe,
386 FullMatrix<double> &interpolation_matrix) const
387{
388 // We don't know how to do this properly, yet.
389 // However, for SolutionTransfer to work we need to provide an implementation
390 // for the case that the x_source_fe is identical to this FE
391 using FEQBUBBLES = FE_Q_Bubbles<dim, spacedim>;
392
394 (x_source_fe.get_name().find("FE_Q_Bubbles<") == 0) ||
395 (dynamic_cast<const FEQBUBBLES *>(&x_source_fe) != nullptr),
397 Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
398 ExcDimensionMismatch(interpolation_matrix.m(),
399 this->n_dofs_per_cell()));
400 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
401 ExcDimensionMismatch(interpolation_matrix.m(),
402 x_source_fe.n_dofs_per_cell()));
403
404 // Provide a short cut in case we are just inquiring the identity
405 auto casted_fe = dynamic_cast<const FEQBUBBLES *>(&x_source_fe);
406 if (casted_fe != nullptr && casted_fe->degree == this->degree)
407 for (unsigned int i = 0; i < interpolation_matrix.m(); ++i)
408 interpolation_matrix.set(i, i, 1.);
409 // else we need to do more...
410 else
411 Assert(
412 false,
413 (typename FiniteElement<dim,
414 spacedim>::ExcInterpolationNotImplemented()));
415}
416
417
418
419template <int dim, int spacedim>
420std::vector<bool>
422{
423 const unsigned int n_cont_dofs = Utilities::fixed_power<dim>(q_deg + 1);
424 const unsigned int n_bubbles = (q_deg <= 1 ? 1 : dim);
425 return std::vector<bool>(n_cont_dofs + n_bubbles, true);
426}
427
428
429
430template <int dim, int spacedim>
431std::vector<unsigned int>
433{
434 std::vector<unsigned int> dpo(dim + 1, 1U);
435 for (unsigned int i = 1; i < dpo.size(); ++i)
436 dpo[i] = dpo[i - 1] * (q_deg - 1);
437
438 // Then add the bubble functions; they are all associated with the
439 // cell interior
440 dpo[dim] += (q_deg <= 1 ? 1 : dim);
441 return dpo;
442}
443
444
445
446template <int dim, int spacedim>
447bool
449 const unsigned int shape_index,
450 const unsigned int face_index) const
451{
452 // discontinuous functions have no support on faces
453 if (shape_index >= this->n_dofs_per_cell() - n_bubbles)
454 return false;
455 else
457 face_index);
458}
459
460
461
462template <int dim, int spacedim>
463const FullMatrix<double> &
465 const unsigned int child,
466 const RefinementCase<dim> &refinement_case) const
467{
468 AssertIndexRange(refinement_case,
472 "Prolongation matrices are only available for refined cells!"));
473 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
474
475 Assert(this->prolongation[refinement_case - 1][child].n() != 0,
476 ExcMessage("This prolongation matrix has not been computed yet!"));
477 // finally return the matrix
478 return this->prolongation[refinement_case - 1][child];
479}
480
481
482
483template <int dim, int spacedim>
484const FullMatrix<double> &
486 const unsigned int child,
487 const RefinementCase<dim> &refinement_case) const
488{
489 AssertIndexRange(refinement_case,
493 "Restriction matrices are only available for refined cells!"));
494 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
495
496 Assert(this->restriction[refinement_case - 1][child].n() != 0,
497 ExcMessage("This restriction matrix has not been computed yet!"));
498
499 // finally return the matrix
500 return this->restriction[refinement_case - 1][child];
501}
502
503
504
505template <int dim, int spacedim>
508 const FiniteElement<dim, spacedim> &fe_other,
509 const unsigned int codim) const
510{
511 Assert(codim <= dim, ExcImpossibleInDim(dim));
512
513 // vertex/line/face domination
514 // (if fe_other is derived from FE_DGQ)
515 // ------------------------------------
516 if (codim > 0)
517 if (dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other) != nullptr)
518 // there are no requirements between continuous and discontinuous elements
520
521 // vertex/line/face domination
522 // (if fe_other is not derived from FE_DGQ)
523 // & cell domination
524 // ----------------------------------------
525 if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
526 dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
527 {
528 if (this->degree < fe_bubbles_other->degree)
530 else if (this->degree == fe_bubbles_other->degree)
532 else
534 }
535 else if (const FE_Nothing<dim> *fe_nothing =
536 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
537 {
538 if (fe_nothing->is_dominating())
540 else
541 // the FE_Nothing has no degrees of freedom and it is typically used
542 // in a context where we don't require any continuity along the
543 // interface
545 }
546
549}
550
551
552// explicit instantiations
553#include "fe_q_bubbles.inst"
554
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:2397
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:2435
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2411
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:111
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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:284
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
static unsigned int n_children(const RefinementCase< dim > &refinement_case)