Reference documentation for deal.II version GIT relicensing-31-gf283a0fc39 2024-02-29 06:50: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_dgp_monomial.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2004 - 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
16
19#include <deal.II/fe/fe_tools.h>
20
21#include <memory>
22#include <sstream>
23
24
26
27namespace internal
28{
30 {
31 namespace
32 {
33 // storage of hand-chosen support
34 // points
35 //
36 // For dim=2, dofs_per_cell of
37 // FE_DGPMonomial(k) is given by
38 // 0.5(k+1)(k+2), i.e.
39 //
40 // k 0 1 2 3 4 5 6 7
41 // dofs 1 3 6 10 15 21 28 36
42 //
43 // indirect access of unit points:
44 // the points for degree k are
45 // located at
46 //
47 // points[start_index[k]..start_index[k+1]-1]
48 const unsigned int start_index2d[6] = {0, 1, 4, 10, 20, 35};
49 const double points2d[35][2] = {
50 {0, 0}, {0, 0}, {1, 0}, {0, 1}, {0, 0},
51 {1, 0}, {0, 1}, {1, 1}, {0.5, 0}, {0, 0.5},
52 {0, 0}, {1, 0}, {0, 1}, {1, 1}, {1. / 3., 0},
53 {2. / 3., 0}, {0, 1. / 3.}, {0, 2. / 3.}, {0.5, 1}, {1, 0.5},
54 {0, 0}, {1, 0}, {0, 1}, {1, 1}, {0.25, 0},
55 {0.5, 0}, {0.75, 0}, {0, 0.25}, {0, 0.5}, {0, 0.75},
56 {1. / 3., 1}, {2. / 3., 1}, {1, 1. / 3.}, {1, 2. / 3.}, {0.5, 0.5}};
57
58 // For dim=3, dofs_per_cell of
59 // FE_DGPMonomial(k) is given by
60 // 1./6.(k+1)(k+2)(k+3), i.e.
61 //
62 // k 0 1 2 3 4 5 6 7
63 // dofs 1 4 10 20 35 56 84 120
64 const unsigned int start_index3d[6] = {0, 1, 5, 15 /*,35*/};
65 const double points3d[35][3] = {{0, 0, 0},
66 {0, 0, 0},
67 {1, 0, 0},
68 {0, 1, 0},
69 {0, 0, 1},
70 {0, 0, 0},
71 {1, 0, 0},
72 {0, 1, 0},
73 {0, 0, 1},
74 {0.5, 0, 0},
75 {0, 0.5, 0},
76 {0, 0, 0.5},
77 {1, 1, 0},
78 {1, 0, 1},
79 {0, 1, 1}};
80
81
82 template <int dim>
83 void
84 generate_unit_points(const unsigned int, std::vector<Point<dim>> &);
85
86 template <>
87 void
88 generate_unit_points(const unsigned int k, std::vector<Point<1>> &p)
89 {
90 Assert(p.size() == k + 1, ExcDimensionMismatch(p.size(), k + 1));
91 const double h = 1. / k;
92 for (unsigned int i = 0; i < p.size(); ++i)
93 p[i][0] = i * h;
94 }
95
96 template <>
97 void
98 generate_unit_points(const unsigned int k, std::vector<Point<2>> &p)
99 {
100 Assert(k <= 4, ExcNotImplemented());
101 Assert(p.size() == start_index2d[k + 1] - start_index2d[k],
103 for (unsigned int i = 0; i < p.size(); ++i)
104 {
105 p[i][0] = points2d[start_index2d[k] + i][0];
106 p[i][1] = points2d[start_index2d[k] + i][1];
107 }
108 }
109
110 template <>
111 void
112 generate_unit_points(const unsigned int k, std::vector<Point<3>> &p)
113 {
114 Assert(k <= 2, ExcNotImplemented());
115 Assert(p.size() == start_index3d[k + 1] - start_index3d[k],
117 for (unsigned int i = 0; i < p.size(); ++i)
118 {
119 p[i][0] = points3d[start_index3d[k] + i][0];
120 p[i][1] = points3d[start_index3d[k] + i][1];
121 p[i][2] = points3d[start_index3d[k] + i][2];
122 }
123 }
124 } // namespace
125 } // namespace FE_DGPMonomial
126} // namespace internal
127
128
129
130template <int dim>
131FE_DGPMonomial<dim>::FE_DGPMonomial(const unsigned int degree)
132 : FE_Poly<dim>(PolynomialsP<dim>(degree),
133 FiniteElementData<dim>(get_dpo_vector(degree),
134 1,
135 degree,
136 FiniteElementData<dim>::L2),
137 std::vector<bool>(
138 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
139 .n_dofs_per_cell(),
140 true),
141 std::vector<ComponentMask>(
142 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
143 .n_dofs_per_cell(),
144 ComponentMask(std::vector<bool>(1, true))))
145{
146 Assert(this->poly_space->n() == this->n_dofs_per_cell(), ExcInternalError());
147 Assert(this->poly_space->degree() == this->degree, ExcInternalError());
148
149 // DG doesn't have constraints, so
150 // leave them empty
151
152 // Reinit the vectors of
153 // restriction and prolongation
154 // matrices to the right sizes
156 // Fill prolongation matrices with embedding operators
158 // Fill restriction matrices with L2-projection
160}
161
162
163
164template <int dim>
165std::string
167{
168 // note that the
169 // FETools::get_fe_by_name
170 // function depends on the
171 // particular format of the string
172 // this function returns, so they
173 // have to be kept in synch
174
175 std::ostringstream namebuf;
176 namebuf << "FE_DGPMonomial<" << dim << ">(" << this->degree << ")";
177
178 return namebuf.str();
179}
180
181
182
183template <int dim>
184std::unique_ptr<FiniteElement<dim, dim>>
186{
187 return std::make_unique<FE_DGPMonomial<dim>>(*this);
188}
189
190
191
192// TODO: Remove this function and use the one in FETools, if needed
193template <int dim>
194void
196 const FiniteElement<dim> &source_fe,
197 FullMatrix<double> &interpolation_matrix) const
198{
199 const FE_DGPMonomial<dim> *source_dgp_monomial =
200 dynamic_cast<const FE_DGPMonomial<dim> *>(&source_fe);
201
202 if (source_dgp_monomial)
203 {
204 // ok, source_fe is a DGP_Monomial
205 // element. Then, the interpolation
206 // matrix is simple
207 const unsigned int m = interpolation_matrix.m();
208 const unsigned int n = interpolation_matrix.n();
209 (void)m;
210 (void)n;
211 Assert(m == this->n_dofs_per_cell(),
212 ExcDimensionMismatch(m, this->n_dofs_per_cell()));
213 Assert(n == source_dgp_monomial->n_dofs_per_cell(),
214 ExcDimensionMismatch(n, source_dgp_monomial->n_dofs_per_cell()));
215
216 const unsigned int min_mn =
217 interpolation_matrix.m() < interpolation_matrix.n() ?
218 interpolation_matrix.m() :
219 interpolation_matrix.n();
220
221 for (unsigned int i = 0; i < min_mn; ++i)
222 interpolation_matrix(i, i) = 1.;
223 }
224 else
225 {
226 std::vector<Point<dim>> unit_points(this->n_dofs_per_cell());
227 internal::FE_DGPMonomial::generate_unit_points(this->degree, unit_points);
228
229 FullMatrix<double> source_fe_matrix(unit_points.size(),
230 source_fe.n_dofs_per_cell());
231 for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
232 for (unsigned int k = 0; k < unit_points.size(); ++k)
233 source_fe_matrix(k, j) = source_fe.shape_value(j, unit_points[k]);
234
235 FullMatrix<double> this_matrix(this->n_dofs_per_cell(),
236 this->n_dofs_per_cell());
237 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
238 for (unsigned int k = 0; k < unit_points.size(); ++k)
239 this_matrix(k, j) =
240 this->poly_space->compute_value(j, unit_points[k]);
241
242 this_matrix.gauss_jordan();
243
244 this_matrix.mmult(interpolation_matrix, source_fe_matrix);
245 }
246}
247
248
249
250template <int dim>
251void
256
257
258//---------------------------------------------------------------------------
259// Auxiliary functions
260//---------------------------------------------------------------------------
261
262
263template <int dim>
264std::vector<unsigned int>
266{
267 std::vector<unsigned int> dpo(dim + 1, 0U);
268 dpo[dim] = deg + 1;
269 for (unsigned int i = 1; i < dim; ++i)
270 {
271 dpo[dim] *= deg + 1 + i;
272 dpo[dim] /= i + 1;
273 }
274 return dpo;
275}
276
277
278template <int dim>
279void
281 const FiniteElement<dim> &x_source_fe,
282 FullMatrix<double> &interpolation_matrix,
283 const unsigned int) const
284{
285 // this is only implemented, if the source
286 // FE is also a DGPMonomial element. in that case,
287 // both elements have no dofs on their
288 // faces and the face interpolation matrix
289 // is necessarily empty -- i.e. there isn't
290 // much we need to do here.
291 (void)interpolation_matrix;
292 AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
293 (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
294 nullptr),
296
297 Assert(interpolation_matrix.m() == 0,
298 ExcDimensionMismatch(interpolation_matrix.m(), 0));
299 Assert(interpolation_matrix.n() == 0,
300 ExcDimensionMismatch(interpolation_matrix.n(), 0));
301}
302
303
304
305template <int dim>
306void
308 const FiniteElement<dim> &x_source_fe,
309 const unsigned int,
310 FullMatrix<double> &interpolation_matrix,
311 const unsigned int) const
312{
313 // this is only implemented, if the source
314 // FE is also a DGPMonomial element. in that case,
315 // both elements have no dofs on their
316 // faces and the face interpolation matrix
317 // is necessarily empty -- i.e. there isn't
318 // much we need to do here.
319 (void)interpolation_matrix;
320 AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
321 (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
322 nullptr),
324
325 Assert(interpolation_matrix.m() == 0,
326 ExcDimensionMismatch(interpolation_matrix.m(), 0));
327 Assert(interpolation_matrix.n() == 0,
328 ExcDimensionMismatch(interpolation_matrix.n(), 0));
329}
330
331
332
333template <int dim>
334bool
339
340
341
342template <int dim>
343std::vector<std::pair<unsigned int, unsigned int>>
345 const FiniteElement<dim> &fe_other) const
346{
347 // there are no such constraints for DGPMonomial
348 // elements at all
349 if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
350 return std::vector<std::pair<unsigned int, unsigned int>>();
351 else
352 {
354 return std::vector<std::pair<unsigned int, unsigned int>>();
355 }
356}
357
358
359
360template <int dim>
361std::vector<std::pair<unsigned int, unsigned int>>
363 const FiniteElement<dim> &fe_other) const
364{
365 // there are no such constraints for DGPMonomial
366 // elements at all
367 if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
368 return std::vector<std::pair<unsigned int, unsigned int>>();
369 else
370 {
372 return std::vector<std::pair<unsigned int, unsigned int>>();
373 }
374}
375
376
377
378template <int dim>
379std::vector<std::pair<unsigned int, unsigned int>>
381 const unsigned int) const
382{
383 // there are no such constraints for DGPMonomial
384 // elements at all
385 if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
386 return std::vector<std::pair<unsigned int, unsigned int>>();
387 else
388 {
390 return std::vector<std::pair<unsigned int, unsigned int>>();
391 }
392}
393
394
395
396template <int dim>
399 const unsigned int codim) const
400{
401 Assert(codim <= dim, ExcImpossibleInDim(dim));
402
403 // vertex/line/face domination
404 // ---------------------------
405 if (codim > 0)
406 // this is a discontinuous element, so by definition there will
407 // be no constraints wherever this element comes together with
408 // any other kind of element
410
411 // cell domination
412 // ---------------
413 if (const FE_DGPMonomial<dim> *fe_monomial_other =
414 dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other))
415 {
416 if (this->degree < fe_monomial_other->degree)
418 else if (this->degree == fe_monomial_other->degree)
420 else
422 }
423 else if (const FE_Nothing<dim> *fe_nothing =
424 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
425 {
426 if (fe_nothing->is_dominating())
428 else
429 // the FE_Nothing has no degrees of freedom and it is typically used
430 // in a context where we don't require any continuity along the
431 // interface
433 }
434
437}
438
439
440
441template <>
442bool
444 const unsigned int face_index) const
445{
446 return face_index == 1 || (face_index == 0 && this->degree == 0);
447}
448
449
450
451template <>
452bool
453FE_DGPMonomial<2>::has_support_on_face(const unsigned int shape_index,
454 const unsigned int face_index) const
455{
456 bool support_on_face = false;
457 if (face_index == 1 || face_index == 2)
458 support_on_face = true;
459 else
460 {
461 auto *const polynomial_space_p =
462 dynamic_cast<PolynomialsP<2> *>(this->poly_space.get());
463 Assert(polynomial_space_p != nullptr, ExcInternalError());
464 const std::array<unsigned int, 2> degrees =
465 polynomial_space_p->directional_degrees(shape_index);
466
467 if ((face_index == 0 && degrees[1] == 0) ||
468 (face_index == 3 && degrees[0] == 0))
469 support_on_face = true;
470 }
471 return support_on_face;
472}
473
474
475
476template <>
477bool
478FE_DGPMonomial<3>::has_support_on_face(const unsigned int shape_index,
479 const unsigned int face_index) const
480{
481 bool support_on_face = false;
482 if (face_index == 1 || face_index == 3 || face_index == 4)
483 support_on_face = true;
484 else
485 {
486 auto *const polynomial_space_p =
487 dynamic_cast<PolynomialsP<3> *>(this->poly_space.get());
488 Assert(polynomial_space_p != nullptr, ExcInternalError());
489 const std::array<unsigned int, 3> degrees =
490 polynomial_space_p->directional_degrees(shape_index);
491
492 if ((face_index == 0 && degrees[1] == 0) ||
493 (face_index == 2 && degrees[2] == 0) ||
494 (face_index == 5 && degrees[0] == 0))
495 support_on_face = true;
496 }
497 return support_on_face;
498}
499
500
501
502template <int dim>
503std::size_t
509
510
511
512// explicit instantiations
513#include "fe_dgp_monomial.inst"
514
515
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
FE_DGPMonomial(const unsigned int p)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual void get_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::size_t memory_consumption() const override
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
virtual bool hp_constraints_are_implemented() const override
virtual std::string get_name() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
void initialize_restriction()
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Definition fe_poly.h:532
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)
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
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
size_type n() const
void gauss_jordan()
size_type m() const
Definition point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
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)
STL namespace.