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