Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_dgp_monomial.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2004 - 2018 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 #include <deal.II/base/std_cxx14/memory.h>
18 
19 #include <deal.II/fe/fe_dgp_monomial.h>
20 #include <deal.II/fe/fe_nothing.h>
21 #include <deal.II/fe/fe_tools.h>
22 
23 #include <sstream>
24 
25 
26 DEAL_II_NAMESPACE_OPEN
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<PolynomialsP<dim>, dim>(
134  PolynomialsP<dim>(degree),
135  FiniteElementData<dim>(get_dpo_vector(degree),
136  1,
137  degree,
138  FiniteElementData<dim>::L2),
139  std::vector<bool>(
140  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
141  true),
142  std::vector<ComponentMask>(
143  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
144  std::vector<bool>(1, true)))
145 {
146  Assert(this->poly_space.n() == this->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 
164 template <int dim>
165 std::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 
183 template <int dim>
184 std::unique_ptr<FiniteElement<dim, dim>>
186 {
187  return std_cxx14::make_unique<FE_DGPMonomial<dim>>(*this);
188 }
189 
190 
191 
192 // TODO: Remove this function and use the one in FETools, if needed
193 template <int dim>
194 void
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->dofs_per_cell,
212  ExcDimensionMismatch(m, this->dofs_per_cell));
213  Assert(n == source_dgp_monomial->dofs_per_cell,
214  ExcDimensionMismatch(n, source_dgp_monomial->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->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.dofs_per_cell);
231  for (unsigned int j = 0; j < source_fe.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->dofs_per_cell, this->dofs_per_cell);
236  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
237  for (unsigned int k = 0; k < unit_points.size(); ++k)
238  this_matrix(k, j) = this->poly_space.compute_value(j, unit_points[k]);
239 
240  this_matrix.gauss_jordan();
241 
242  this_matrix.mmult(interpolation_matrix, source_fe_matrix);
243  }
244 }
245 
246 
247 
248 template <int dim>
249 void
251 {
252  Assert(false, ExcNotImplemented());
253 }
254 
255 
256 //---------------------------------------------------------------------------
257 // Auxiliary functions
258 //---------------------------------------------------------------------------
259 
260 
261 template <int dim>
262 std::vector<unsigned int>
263 FE_DGPMonomial<dim>::get_dpo_vector(const unsigned int deg)
264 {
265  std::vector<unsigned int> dpo(dim + 1, 0U);
266  dpo[dim] = deg + 1;
267  for (unsigned int i = 1; i < dim; ++i)
268  {
269  dpo[dim] *= deg + 1 + i;
270  dpo[dim] /= i + 1;
271  }
272  return dpo;
273 }
274 
275 
276 template <int dim>
277 void
279  const FiniteElement<dim> &x_source_fe,
280  FullMatrix<double> & interpolation_matrix) const
281 {
282  // this is only implemented, if the source
283  // FE is also a DGPMonomial element. in that case,
284  // both elements have no dofs on their
285  // faces and the face interpolation matrix
286  // is necessarily empty -- i.e. there isn't
287  // much we need to do here.
288  (void)interpolation_matrix;
289  AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
290  (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
291  nullptr),
293 
294  Assert(interpolation_matrix.m() == 0,
295  ExcDimensionMismatch(interpolation_matrix.m(), 0));
296  Assert(interpolation_matrix.n() == 0,
297  ExcDimensionMismatch(interpolation_matrix.n(), 0));
298 }
299 
300 
301 
302 template <int dim>
303 void
305  const FiniteElement<dim> &x_source_fe,
306  const unsigned int,
307  FullMatrix<double> &interpolation_matrix) const
308 {
309  // this is only implemented, if the source
310  // FE is also a DGPMonomial element. in that case,
311  // both elements have no dofs on their
312  // faces and the face interpolation matrix
313  // is necessarily empty -- i.e. there isn't
314  // much we need to do here.
315  (void)interpolation_matrix;
316  AssertThrow((x_source_fe.get_name().find("FE_DGPMonomial<") == 0) ||
317  (dynamic_cast<const FE_DGPMonomial<dim> *>(&x_source_fe) !=
318  nullptr),
320 
321  Assert(interpolation_matrix.m() == 0,
322  ExcDimensionMismatch(interpolation_matrix.m(), 0));
323  Assert(interpolation_matrix.n() == 0,
324  ExcDimensionMismatch(interpolation_matrix.n(), 0));
325 }
326 
327 
328 
329 template <int dim>
330 bool
332 {
333  return true;
334 }
335 
336 
337 
338 template <int dim>
339 std::vector<std::pair<unsigned int, unsigned int>>
341  const FiniteElement<dim> &fe_other) const
342 {
343  // there are no such constraints for DGPMonomial
344  // elements at all
345  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
346  return std::vector<std::pair<unsigned int, unsigned int>>();
347  else
348  {
349  Assert(false, ExcNotImplemented());
350  return std::vector<std::pair<unsigned int, unsigned int>>();
351  }
352 }
353 
354 
355 
356 template <int dim>
357 std::vector<std::pair<unsigned int, unsigned int>>
359  const FiniteElement<dim> &fe_other) const
360 {
361  // there are no such constraints for DGPMonomial
362  // elements at all
363  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
364  return std::vector<std::pair<unsigned int, unsigned int>>();
365  else
366  {
367  Assert(false, ExcNotImplemented());
368  return std::vector<std::pair<unsigned int, unsigned int>>();
369  }
370 }
371 
372 
373 
374 template <int dim>
375 std::vector<std::pair<unsigned int, unsigned int>>
377  const FiniteElement<dim> &fe_other) const
378 {
379  // there are no such constraints for DGPMonomial
380  // elements at all
381  if (dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other) != nullptr)
382  return std::vector<std::pair<unsigned int, unsigned int>>();
383  else
384  {
385  Assert(false, ExcNotImplemented());
386  return std::vector<std::pair<unsigned int, unsigned int>>();
387  }
388 }
389 
390 
391 
392 template <int dim>
395  const unsigned int codim) const
396 {
397  Assert(codim <= dim, ExcImpossibleInDim(dim));
398 
399  // vertex/line/face domination
400  // ---------------------------
401  if (codim > 0)
402  // this is a discontinuous element, so by definition there will
403  // be no constraints wherever this element comes together with
404  // any other kind of element
406 
407  // cell domination
408  // ---------------
409  if (const FE_DGPMonomial<dim> *fe_monomial_other =
410  dynamic_cast<const FE_DGPMonomial<dim> *>(&fe_other))
411  {
412  if (this->degree < fe_monomial_other->degree)
414  else if (this->degree == fe_monomial_other->degree)
416  else
418  }
419  else if (const FE_Nothing<dim> *fe_nothing =
420  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
421  {
422  if (fe_nothing->is_dominating())
424  else
425  // the FE_Nothing has no degrees of freedom and it is typically used
426  // in a context where we don't require any continuity along the
427  // interface
429  }
430 
431  Assert(false, ExcNotImplemented());
433 }
434 
435 
436 
437 template <>
438 bool
440  const unsigned int face_index) const
441 {
442  return face_index == 1 || (face_index == 0 && this->degree == 0);
443 }
444 
445 
446 
447 template <>
448 bool
449 FE_DGPMonomial<2>::has_support_on_face(const unsigned int shape_index,
450  const unsigned int face_index) const
451 {
452  bool support_on_face = false;
453  if (face_index == 1 || face_index == 2)
454  support_on_face = true;
455  else
456  {
457  const std::array<unsigned int, 2> degrees =
458  this->poly_space.directional_degrees(shape_index);
459 
460  if ((face_index == 0 && degrees[1] == 0) ||
461  (face_index == 3 && degrees[0] == 0))
462  support_on_face = true;
463  }
464  return support_on_face;
465 }
466 
467 
468 
469 template <>
470 bool
471 FE_DGPMonomial<3>::has_support_on_face(const unsigned int shape_index,
472  const unsigned int face_index) const
473 {
474  bool support_on_face = false;
475  if (face_index == 1 || face_index == 3 || face_index == 4)
476  support_on_face = true;
477  else
478  {
479  const std::array<unsigned int, 3> degrees =
480  this->poly_space.directional_degrees(shape_index);
481 
482  if ((face_index == 0 && degrees[1] == 0) ||
483  (face_index == 2 && degrees[2] == 0) ||
484  (face_index == 5 && degrees[0] == 0))
485  support_on_face = true;
486  }
487  return support_on_face;
488 }
489 
490 
491 
492 template <int dim>
493 std::size_t
495 {
496  Assert(false, ExcNotImplemented());
497  return 0;
498 }
499 
500 
501 
502 // explicit instantiations
503 #include "fe_dgp_monomial.inst"
504 
505 
506 DEAL_II_NAMESPACE_CLOSE
size_type m() const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other) const override
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2442
FE_DGPMonomial(const unsigned int p)
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 gauss_jordan()
const unsigned int degree
Definition: fe_base.h:298
virtual std::size_t memory_consumption() const override
STL namespace.
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual void get_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
void initialize_restriction()
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2456
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
virtual std::string get_name() const =0
const unsigned int dofs_per_cell
Definition: fe_base.h:282
virtual bool hp_constraints_are_implemented() const override
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
virtual double shape_value(const unsigned int i, const Point< dim > &p) const
Definition: fe.cc:159
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:276
static ::ExceptionBase & ExcNotImplemented()
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
PolynomialsP< dim > poly_space
Definition: fe_poly.h:482
virtual std::string get_name() const override
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
static ::ExceptionBase & ExcInternalError()