Reference documentation for deal.II version GIT f0f8c7fe18 2023-03-21 21:25: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\}}\)
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 
19 #include <deal.II/fe/fe_nothing.h>
20 #include <deal.II/fe/fe_tools.h>
21 
22 #include <memory>
23 #include <sstream>
24 
25 
27 
28 namespace internal
29 {
30  namespace FE_DGPMonomial
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],
103  ExcInternalError());
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],
117  ExcInternalError());
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 
131 template <int dim>
132 FE_DGPMonomial<dim>::FE_DGPMonomial(const unsigned int degree)
133  : FE_Poly<dim>(PolynomialsP<dim>(degree),
134  FiniteElementData<dim>(get_dpo_vector(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 
165 template <int dim>
166 std::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 
184 template <int dim>
185 std::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
194 template <int dim>
195 void
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 
251 template <int dim>
252 void
254 {
255  Assert(false, ExcNotImplemented());
256 }
257 
258 
259 //---------------------------------------------------------------------------
260 // Auxiliary functions
261 //---------------------------------------------------------------------------
262 
263 
264 template <int dim>
265 std::vector<unsigned int>
266 FE_DGPMonomial<dim>::get_dpo_vector(const unsigned int deg)
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 
279 template <int dim>
280 void
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 
306 template <int dim>
307 void
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 
334 template <int dim>
335 bool
337 {
338  return true;
339 }
340 
341 
342 
343 template <int dim>
344 std::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 
361 template <int dim>
362 std::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 
379 template <int dim>
380 std::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 
397 template <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 
442 template <>
443 bool
445  const unsigned int face_index) const
446 {
447  return face_index == 1 || (face_index == 0 && this->degree == 0);
448 }
449 
450 
451 
452 template <>
453 bool
454 FE_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 
477 template <>
478 bool
479 FE_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 
503 template <int dim>
504 std::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:533
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)
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
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
size_type n() const
void gauss_jordan()
size_type m() const
Definition: point.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
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