Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_dgq.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2001 - 2019 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/quadrature.h>
18 #include <deal.II/base/quadrature_lib.h>
19 #include <deal.II/base/std_cxx14/memory.h>
20 
21 #include <deal.II/fe/fe.h>
22 #include <deal.II/fe/fe_bernstein.h>
23 #include <deal.II/fe/fe_dgq.h>
24 #include <deal.II/fe/fe_nothing.h>
25 #include <deal.II/fe/fe_q.h>
26 #include <deal.II/fe/fe_q_bubbles.h>
27 #include <deal.II/fe/fe_q_dg0.h>
28 #include <deal.II/fe/fe_q_hierarchical.h>
29 #include <deal.II/fe/fe_q_iso_q1.h>
30 #include <deal.II/fe/fe_tools.h>
31 
32 #include <deal.II/lac/vector.h>
33 
34 #include <iostream>
35 #include <sstream>
36 
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 
41 namespace internal
42 {
43  namespace FE_DGQ
44  {
45  namespace
46  {
47  std::vector<Point<1>>
48  get_QGaussLobatto_points(const unsigned int degree)
49  {
50  if (degree > 0)
51  return QGaussLobatto<1>(degree + 1).get_points();
52  else
53  return std::vector<Point<1>>(1, Point<1>(0.5));
54  }
55  } // namespace
56  } // namespace FE_DGQ
57 } // namespace internal
58 
59 
60 
61 template <int dim, int spacedim>
62 FE_DGQ<dim, spacedim>::FE_DGQ(const unsigned int degree)
63  : FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>(
65  Polynomials::generate_complete_Lagrange_basis(
66  internal::FE_DGQ::get_QGaussLobatto_points(degree))),
67  FiniteElementData<dim>(get_dpo_vector(degree),
68  1,
69  degree,
70  FiniteElementData<dim>::L2),
71  std::vector<bool>(
72  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
73  true),
74  std::vector<ComponentMask>(
75  FiniteElementData<dim>(get_dpo_vector(degree), 1, degree).dofs_per_cell,
76  std::vector<bool>(1, true)))
77 {
78  // Compute support points, which are the tensor product of the Lagrange
79  // interpolation points in the constructor.
80  this->unit_support_points =
81  Quadrature<dim>(internal::FE_DGQ::get_QGaussLobatto_points(degree))
82  .get_points();
83 
84  // do not initialize embedding and restriction here. these matrices are
85  // initialized on demand in get_restriction_matrix and
86  // get_prolongation_matrix
87 
88  // note: no face support points for DG elements
89 }
90 
91 
92 
93 template <int dim, int spacedim>
95  const std::vector<Polynomials::Polynomial<double>> &polynomials)
96  : FE_Poly<TensorProductPolynomials<dim>, dim, spacedim>(
97  TensorProductPolynomials<dim>(polynomials),
98  FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
99  1,
100  polynomials.size() - 1,
101  FiniteElementData<dim>::L2),
102  std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(
103  polynomials.size() - 1),
104  1,
105  polynomials.size() - 1)
106  .dofs_per_cell,
107  true),
108  std::vector<ComponentMask>(
109  FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
110  1,
111  polynomials.size() - 1)
112  .dofs_per_cell,
113  std::vector<bool>(1, true)))
114 {
115  // No support points can be defined in general. Derived classes might define
116  // support points like the class FE_DGQArbitraryNodes
117 
118  // do not initialize embedding and restriction here. these matrices are
119  // initialized on demand in get_restriction_matrix and
120  // get_prolongation_matrix
121 }
122 
123 
124 
125 template <int dim, int spacedim>
126 std::string
128 {
129  // note that the FETools::get_fe_by_name function depends on the
130  // particular format of the string this function returns, so they have to be
131  // kept in sync
132 
133  std::ostringstream namebuf;
134  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
135  << this->degree << ")";
136  return namebuf.str();
137 }
138 
139 
140 
141 template <int dim, int spacedim>
142 void
144  const std::vector<Vector<double>> &support_point_values,
145  std::vector<double> & nodal_values) const
146 {
147  AssertDimension(support_point_values.size(),
148  this->get_unit_support_points().size());
149  AssertDimension(support_point_values.size(), nodal_values.size());
150  AssertDimension(this->dofs_per_cell, nodal_values.size());
151 
152  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
153  {
154  AssertDimension(support_point_values[i].size(), 1);
155 
156  nodal_values[i] = support_point_values[i](0);
157  }
158 }
159 
160 
161 
162 template <int dim, int spacedim>
163 std::unique_ptr<FiniteElement<dim, spacedim>>
165 {
166  return std_cxx14::make_unique<FE_DGQ<dim, spacedim>>(*this);
167 }
168 
169 
170 //---------------------------------------------------------------------------
171 // Auxiliary functions
172 //---------------------------------------------------------------------------
173 
174 
175 template <int dim, int spacedim>
176 std::vector<unsigned int>
178 {
179  std::vector<unsigned int> dpo(dim + 1, 0U);
180  dpo[dim] = deg + 1;
181  for (unsigned int i = 1; i < dim; ++i)
182  dpo[dim] *= deg + 1;
183  return dpo;
184 }
185 
186 
187 
188 template <int dim, int spacedim>
189 void
191  const char direction) const
192 {
193  const unsigned int n = this->degree + 1;
194  unsigned int s = n;
195  for (unsigned int i = 1; i < dim; ++i)
196  s *= n;
197  numbers.resize(s);
198 
199  unsigned int l = 0;
200 
201  if (dim == 1)
202  {
203  // Mirror around midpoint
204  for (unsigned int i = n; i > 0;)
205  numbers[l++] = --i;
206  }
207  else
208  {
209  switch (direction)
210  {
211  // Rotate xy-plane
212  // counter-clockwise
213  case 'z':
214  for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
215  for (unsigned int j = 0; j < n; ++j)
216  for (unsigned int i = 0; i < n; ++i)
217  {
218  unsigned int k = n * i - j + n - 1 + n * n * iz;
219  numbers[l++] = k;
220  }
221  break;
222  // Rotate xy-plane
223  // clockwise
224  case 'Z':
225  for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
226  for (unsigned int iy = 0; iy < n; ++iy)
227  for (unsigned int ix = 0; ix < n; ++ix)
228  {
229  unsigned int k = n * ix - iy + n - 1 + n * n * iz;
230  numbers[k] = l++;
231  }
232  break;
233  // Rotate yz-plane
234  // counter-clockwise
235  case 'x':
236  Assert(dim > 2, ExcDimensionMismatch(dim, 3));
237  for (unsigned int iz = 0; iz < n; ++iz)
238  for (unsigned int iy = 0; iy < n; ++iy)
239  for (unsigned int ix = 0; ix < n; ++ix)
240  {
241  unsigned int k = n * (n * iy - iz + n - 1) + ix;
242  numbers[l++] = k;
243  }
244  break;
245  // Rotate yz-plane
246  // clockwise
247  case 'X':
248  Assert(dim > 2, ExcDimensionMismatch(dim, 3));
249  for (unsigned int iz = 0; iz < n; ++iz)
250  for (unsigned int iy = 0; iy < n; ++iy)
251  for (unsigned int ix = 0; ix < n; ++ix)
252  {
253  unsigned int k = n * (n * iy - iz + n - 1) + ix;
254  numbers[k] = l++;
255  }
256  break;
257  default:
258  Assert(false, ExcNotImplemented());
259  }
260  }
261 }
262 
263 
264 
265 template <int dim, int spacedim>
266 void
268  const FiniteElement<dim, spacedim> &x_source_fe,
269  FullMatrix<double> & interpolation_matrix) const
270 {
271  // this is only implemented, if the
272  // source FE is also a
273  // DGQ element
274  using FE = FiniteElement<dim, spacedim>;
275  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
276  nullptr),
277  typename FE::ExcInterpolationNotImplemented());
278 
279  // ok, source is a Q element, so
280  // we will be able to do the work
281  const FE_DGQ<dim, spacedim> &source_fe =
282  dynamic_cast<const FE_DGQ<dim, spacedim> &>(x_source_fe);
283 
284  Assert(interpolation_matrix.m() == this->dofs_per_cell,
285  ExcDimensionMismatch(interpolation_matrix.m(), this->dofs_per_cell));
286  Assert(interpolation_matrix.n() == source_fe.dofs_per_cell,
287  ExcDimensionMismatch(interpolation_matrix.n(),
288  source_fe.dofs_per_cell));
289 
290 
291  // compute the interpolation
292  // matrices in much the same way as
293  // we do for the embedding matrices
294  // from mother to child.
295  FullMatrix<double> cell_interpolation(this->dofs_per_cell,
296  this->dofs_per_cell);
297  FullMatrix<double> source_interpolation(this->dofs_per_cell,
298  source_fe.dofs_per_cell);
299  FullMatrix<double> tmp(this->dofs_per_cell, source_fe.dofs_per_cell);
300  for (unsigned int j = 0; j < this->dofs_per_cell; ++j)
301  {
302  // generate a point on this
303  // cell and evaluate the
304  // shape functions there
305  const Point<dim> p = this->unit_support_points[j];
306  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
307  cell_interpolation(j, i) = this->poly_space.compute_value(i, p);
308 
309  for (unsigned int i = 0; i < source_fe.dofs_per_cell; ++i)
310  source_interpolation(j, i) = source_fe.poly_space.compute_value(i, p);
311  }
312 
313  // then compute the
314  // interpolation matrix matrix
315  // for this coordinate
316  // direction
317  cell_interpolation.gauss_jordan();
318  cell_interpolation.mmult(interpolation_matrix, source_interpolation);
319 
320  // cut off very small values
321  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
322  for (unsigned int j = 0; j < source_fe.dofs_per_cell; ++j)
323  if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
324  interpolation_matrix(i, j) = 0.;
325 
326  // make sure that the row sum of
327  // each of the matrices is 1 at
328  // this point. this must be so
329  // since the shape functions sum up
330  // to 1
331  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
332  {
333  double sum = 0.;
334  for (unsigned int j = 0; j < source_fe.dofs_per_cell; ++j)
335  sum += interpolation_matrix(i, j);
336 
337  Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
338  ExcInternalError());
339  }
340 }
341 
342 
343 
344 template <int dim, int spacedim>
345 void
347  const FiniteElement<dim, spacedim> &x_source_fe,
348  FullMatrix<double> & interpolation_matrix) const
349 {
350  // this is only implemented, if the source
351  // FE is also a DGQ element. in that case,
352  // both elements have no dofs on their
353  // faces and the face interpolation matrix
354  // is necessarily empty -- i.e. there isn't
355  // much we need to do here.
356  (void)interpolation_matrix;
357  using FE = FiniteElement<dim, spacedim>;
358  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
359  nullptr),
360  typename FE::ExcInterpolationNotImplemented());
361 
362  Assert(interpolation_matrix.m() == 0,
363  ExcDimensionMismatch(interpolation_matrix.m(), 0));
364  Assert(interpolation_matrix.n() == 0,
365  ExcDimensionMismatch(interpolation_matrix.m(), 0));
366 }
367 
368 
369 
370 template <int dim, int spacedim>
371 void
373  const FiniteElement<dim, spacedim> &x_source_fe,
374  const unsigned int,
375  FullMatrix<double> &interpolation_matrix) const
376 {
377  // this is only implemented, if the source
378  // FE is also a DGQ element. in that case,
379  // both elements have no dofs on their
380  // faces and the face interpolation matrix
381  // is necessarily empty -- i.e. there isn't
382  // much we need to do here.
383  (void)interpolation_matrix;
384  using FE = FiniteElement<dim, spacedim>;
385  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
386  nullptr),
387  typename FE::ExcInterpolationNotImplemented());
388 
389  Assert(interpolation_matrix.m() == 0,
390  ExcDimensionMismatch(interpolation_matrix.m(), 0));
391  Assert(interpolation_matrix.n() == 0,
392  ExcDimensionMismatch(interpolation_matrix.m(), 0));
393 }
394 
395 
396 
397 template <int dim, int spacedim>
398 const FullMatrix<double> &
400  const unsigned int child,
401  const RefinementCase<dim> &refinement_case) const
402 {
404  ExcIndexRange(refinement_case,
405  0,
407  Assert(refinement_case != RefinementCase<dim>::no_refinement,
408  ExcMessage(
409  "Prolongation matrices are only available for refined cells!"));
410  Assert(child < GeometryInfo<dim>::n_children(refinement_case),
411  ExcIndexRange(child,
412  0,
413  GeometryInfo<dim>::n_children(refinement_case)));
414 
415  // initialization upon first request
416  if (this->prolongation[refinement_case - 1][child].n() == 0)
417  {
418  std::lock_guard<std::mutex> lock(this->mutex);
419 
420  // if matrix got updated while waiting for the lock
421  if (this->prolongation[refinement_case - 1][child].n() ==
422  this->dofs_per_cell)
423  return this->prolongation[refinement_case - 1][child];
424 
425  // now do the work. need to get a non-const version of data in order to
426  // be able to modify them inside a const function
427  FE_DGQ<dim, spacedim> &this_nonconst =
428  const_cast<FE_DGQ<dim, spacedim> &>(*this);
429  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
430  {
431  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
433  isotropic_matrices.back().resize(
435  FullMatrix<double>(this->dofs_per_cell, this->dofs_per_cell));
436  if (dim == spacedim)
438  isotropic_matrices,
439  true);
440  else
442  isotropic_matrices,
443  true);
444  this_nonconst.prolongation[refinement_case - 1].swap(
445  isotropic_matrices.back());
446  }
447  else
448  {
449  // must compute both restriction and prolongation matrices because
450  // we only check for their size and the reinit call initializes them
451  // all
453  if (dim == spacedim)
454  {
456  this_nonconst.prolongation);
458  this_nonconst.restriction);
459  }
460  else
461  {
462  FE_DGQ<dim> tmp(this->degree);
464  this_nonconst.prolongation);
466  this_nonconst.restriction);
467  }
468  }
469  }
470 
471  // finally return the matrix
472  return this->prolongation[refinement_case - 1][child];
473 }
474 
475 
476 
477 template <int dim, int spacedim>
478 const FullMatrix<double> &
480  const unsigned int child,
481  const RefinementCase<dim> &refinement_case) const
482 {
484  ExcIndexRange(refinement_case,
485  0,
487  Assert(refinement_case != RefinementCase<dim>::no_refinement,
488  ExcMessage(
489  "Restriction matrices are only available for refined cells!"));
490  Assert(child < GeometryInfo<dim>::n_children(refinement_case),
491  ExcIndexRange(child,
492  0,
493  GeometryInfo<dim>::n_children(refinement_case)));
494 
495  // initialization upon first request
496  if (this->restriction[refinement_case - 1][child].n() == 0)
497  {
498  std::lock_guard<std::mutex> lock(this->mutex);
499 
500  // if matrix got updated while waiting for the lock...
501  if (this->restriction[refinement_case - 1][child].n() ==
502  this->dofs_per_cell)
503  return this->restriction[refinement_case - 1][child];
504 
505  // now do the work. need to get a non-const version of data in order to
506  // be able to modify them inside a const function
507  FE_DGQ<dim, spacedim> &this_nonconst =
508  const_cast<FE_DGQ<dim, spacedim> &>(*this);
509  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
510  {
511  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
513  isotropic_matrices.back().resize(
515  FullMatrix<double>(this->dofs_per_cell, this->dofs_per_cell));
516  if (dim == spacedim)
518  isotropic_matrices,
519  true);
520  else
522  isotropic_matrices,
523  true);
524  this_nonconst.restriction[refinement_case - 1].swap(
525  isotropic_matrices.back());
526  }
527  else
528  {
529  // must compute both restriction and prolongation matrices because
530  // we only check for their size and the reinit call initializes them
531  // all
533  if (dim == spacedim)
534  {
536  this_nonconst.prolongation);
538  this_nonconst.restriction);
539  }
540  else
541  {
542  FE_DGQ<dim> tmp(this->degree);
544  this_nonconst.prolongation);
546  this_nonconst.restriction);
547  }
548  }
549  }
550 
551  // finally return the matrix
552  return this->restriction[refinement_case - 1][child];
553 }
554 
555 
556 
557 template <int dim, int spacedim>
558 bool
560 {
561  return true;
562 }
563 
564 
565 
566 template <int dim, int spacedim>
567 std::vector<std::pair<unsigned int, unsigned int>>
569  const FiniteElement<dim, spacedim> & /*fe_other*/) const
570 {
571  // this element is discontinuous, so by definition there can
572  // be no identities between its dofs and those of any neighbor
573  // (of whichever type the neighbor may be -- after all, we have
574  // no face dofs on this side to begin with)
575  return std::vector<std::pair<unsigned int, unsigned int>>();
576 }
577 
578 
579 
580 template <int dim, int spacedim>
581 std::vector<std::pair<unsigned int, unsigned int>>
583  const FiniteElement<dim, spacedim> & /*fe_other*/) const
584 {
585  // this element is discontinuous, so by definition there can
586  // be no identities between its dofs and those of any neighbor
587  // (of whichever type the neighbor may be -- after all, we have
588  // no face dofs on this side to begin with)
589  return std::vector<std::pair<unsigned int, unsigned int>>();
590 }
591 
592 
593 
594 template <int dim, int spacedim>
595 std::vector<std::pair<unsigned int, unsigned int>>
597  const FiniteElement<dim, spacedim> & /*fe_other*/) const
598 {
599  // this element is discontinuous, so by definition there can
600  // be no identities between its dofs and those of any neighbor
601  // (of whichever type the neighbor may be -- after all, we have
602  // no face dofs on this side to begin with)
603  return std::vector<std::pair<unsigned int, unsigned int>>();
604 }
605 
606 
607 
608 template <int dim, int spacedim>
611  const FiniteElement<dim, spacedim> &fe_other,
612  const unsigned int codim) const
613 {
614  Assert(codim <= dim, ExcImpossibleInDim(dim));
615 
616  // vertex/line/face domination
617  // ---------------------------
618  if (codim > 0)
619  // this is a discontinuous element, so by definition there will
620  // be no constraints wherever this element comes together with
621  // any other kind of element
623 
624  // cell domination
625  // ---------------
626  // The following block of conditionals is rather ugly, but there is currently
627  // no other way how to deal with a robust comparison of FE_DGQ elements with
628  // relevant others in the current implementation.
629  if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
630  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
631  {
632  if (this->degree < fe_dgq_other->degree)
634  else if (this->degree == fe_dgq_other->degree)
636  else
638  }
639  else if (const FE_Q<dim, spacedim> *fe_q_other =
640  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
641  {
642  if (this->degree < fe_q_other->degree)
644  else if (this->degree == fe_q_other->degree)
646  else
648  }
649  else if (const FE_Bernstein<dim, spacedim> *fe_bernstein_other =
650  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
651  {
652  if (this->degree < fe_bernstein_other->degree)
654  else if (this->degree == fe_bernstein_other->degree)
656  else
658  }
659  else if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
660  dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
661  {
662  if (this->degree < fe_bubbles_other->degree)
664  else if (this->degree == fe_bubbles_other->degree)
666  else
668  }
669  else if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
670  dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
671  {
672  if (this->degree < fe_dg0_other->degree)
674  else if (this->degree == fe_dg0_other->degree)
676  else
678  }
679  else if (const FE_Q_iso_Q1<dim, spacedim> *fe_q_iso_q1_other =
680  dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe_other))
681  {
682  if (this->degree < fe_q_iso_q1_other->degree)
684  else if (this->degree == fe_q_iso_q1_other->degree)
686  else
688  }
689  else if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
690  dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
691  {
692  if (this->degree < fe_hierarchical_other->degree)
694  else if (this->degree == fe_hierarchical_other->degree)
696  else
698  }
699  else if (const FE_Nothing<dim> *fe_nothing =
700  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
701  {
702  if (fe_nothing->is_dominating())
704  else
705  // the FE_Nothing has no degrees of freedom and it is typically used
706  // in a context where we don't require any continuity along the
707  // interface
709  }
710 
711  Assert(false, ExcNotImplemented());
713 }
714 
715 
716 
717 template <int dim, int spacedim>
718 bool
719 FE_DGQ<dim, spacedim>::has_support_on_face(const unsigned int shape_index,
720  const unsigned int face_index) const
721 {
722  Assert(shape_index < this->dofs_per_cell,
723  ExcIndexRange(shape_index, 0, this->dofs_per_cell));
726 
727  unsigned int n = this->degree + 1;
728 
729  // For DGQ elements that do not define support points, we cannot define
730  // whether they have support at the boundary easily, so return true in any
731  // case
732  if (this->unit_support_points.empty())
733  return true;
734 
735  // for DGQ(0) elements or arbitrary node DGQ with support points not located
736  // at the element boundary, the single shape functions is constant and
737  // therefore lives on the boundary
738  bool support_points_on_boundary = true;
739  for (unsigned int d = 0; d < dim; ++d)
740  if (std::abs(this->unit_support_points[0][d]) > 1e-13)
741  support_points_on_boundary = false;
742  for (unsigned int d = 0; d < dim; ++d)
743  if (std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
744  support_points_on_boundary = false;
745  if (support_points_on_boundary == false)
746  return true;
747 
748  unsigned int n2 = n * n;
749 
750  switch (dim)
751  {
752  case 1:
753  {
754  // in 1d, things are simple. since
755  // there is only one degree of
756  // freedom per vertex in this
757  // class, the first is on vertex 0
758  // (==face 0 in some sense), the
759  // second on face 1:
760  return (((shape_index == 0) && (face_index == 0)) ||
761  ((shape_index == this->degree) && (face_index == 1)));
762  }
763 
764  case 2:
765  {
766  if (face_index == 0 && (shape_index % n) == 0)
767  return true;
768  if (face_index == 1 && (shape_index % n) == this->degree)
769  return true;
770  if (face_index == 2 && shape_index < n)
771  return true;
772  if (face_index == 3 && shape_index >= this->dofs_per_cell - n)
773  return true;
774  return false;
775  }
776 
777  case 3:
778  {
779  const unsigned int in2 = shape_index % n2;
780 
781  // x=0
782  if (face_index == 0 && (shape_index % n) == 0)
783  return true;
784  // x=1
785  if (face_index == 1 && (shape_index % n) == n - 1)
786  return true;
787  // y=0
788  if (face_index == 2 && in2 < n)
789  return true;
790  // y=1
791  if (face_index == 3 && in2 >= n2 - n)
792  return true;
793  // z=0
794  if (face_index == 4 && shape_index < n2)
795  return true;
796  // z=1
797  if (face_index == 5 && shape_index >= this->dofs_per_cell - n2)
798  return true;
799  return false;
800  }
801 
802  default:
803  Assert(false, ExcNotImplemented());
804  }
805  return true;
806 }
807 
808 
809 
810 template <int dim, int spacedim>
811 std::pair<Table<2, bool>, std::vector<unsigned int>>
813 {
814  Table<2, bool> constant_modes(1, this->dofs_per_cell);
815  constant_modes.fill(true);
816  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
817  constant_modes, std::vector<unsigned int>(1, 0));
818 }
819 
820 
821 
822 template <int dim, int spacedim>
823 std::size_t
825 {
826  Assert(false, ExcNotImplemented());
827  return 0;
828 }
829 
830 
831 
832 // ------------------------------ FE_DGQArbitraryNodes -----------------------
833 
834 template <int dim, int spacedim>
836  const Quadrature<1> &points)
837  : FE_DGQ<dim, spacedim>(
838  Polynomials::generate_complete_Lagrange_basis(points.get_points()))
839 {
840  Assert(points.size() > 0,
842  this->unit_support_points = Quadrature<dim>(points).get_points();
843 }
844 
845 
846 
847 template <int dim, int spacedim>
848 std::string
850 {
851  // note that the FETools::get_fe_by_name function does not work for
852  // FE_DGQArbitraryNodes since there is no initialization by a degree value.
853  std::ostringstream namebuf;
854  bool equidistant = true;
855  std::vector<double> points(this->degree + 1);
856 
857  std::vector<unsigned int> lexicographic =
858  this->poly_space.get_numbering_inverse();
859  for (unsigned int j = 0; j <= this->degree; j++)
860  points[j] = this->unit_support_points[lexicographic[j]][0];
861 
862  // Check whether the support points are equidistant.
863  for (unsigned int j = 0; j <= this->degree; j++)
864  if (std::abs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
865  {
866  equidistant = false;
867  break;
868  }
869  if (this->degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
870  equidistant = true;
871 
872  if (equidistant == true)
873  {
874  if (this->degree > 2)
875  namebuf << "FE_DGQArbitraryNodes<"
876  << Utilities::dim_string(dim, spacedim)
877  << ">(QIterated(QTrapez()," << this->degree << "))";
878  else
879  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
880  << this->degree << ")";
881  return namebuf.str();
882  }
883 
884  // Check whether the support points come from QGaussLobatto.
885  const QGaussLobatto<1> points_gl(this->degree + 1);
886  bool gauss_lobatto = true;
887  for (unsigned int j = 0; j <= this->degree; j++)
888  if (points[j] != points_gl.point(j)(0))
889  {
890  gauss_lobatto = false;
891  break;
892  }
893 
894  if (gauss_lobatto == true)
895  {
896  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
897  << this->degree << ")";
898  return namebuf.str();
899  }
900 
901  // Check whether the support points come from QGauss.
902  const QGauss<1> points_g(this->degree + 1);
903  bool gauss = true;
904  for (unsigned int j = 0; j <= this->degree; j++)
905  if (points[j] != points_g.point(j)(0))
906  {
907  gauss = false;
908  break;
909  }
910 
911  if (gauss == true)
912  {
913  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
914  << ">(QGauss(" << this->degree + 1 << "))";
915  return namebuf.str();
916  }
917 
918  // Check whether the support points come from QGauss.
919  const QGaussLog<1> points_glog(this->degree + 1);
920  bool gauss_log = true;
921  for (unsigned int j = 0; j <= this->degree; j++)
922  if (points[j] != points_glog.point(j)(0))
923  {
924  gauss_log = false;
925  break;
926  }
927 
928  if (gauss_log == true)
929  {
930  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
931  << ">(QGaussLog(" << this->degree + 1 << "))";
932  return namebuf.str();
933  }
934 
935  // All guesses exhausted
936  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
937  << ">(QUnknownNodes(" << this->degree + 1 << "))";
938  return namebuf.str();
939 }
940 
941 
942 
943 template <int dim, int spacedim>
944 void
947  const std::vector<Vector<double>> &support_point_values,
948  std::vector<double> & nodal_values) const
949 {
950  AssertDimension(support_point_values.size(),
951  this->get_unit_support_points().size());
952  AssertDimension(support_point_values.size(), nodal_values.size());
953  AssertDimension(this->dofs_per_cell, nodal_values.size());
954 
955  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
956  {
957  AssertDimension(support_point_values[i].size(), 1);
958 
959  nodal_values[i] = support_point_values[i](0);
960  }
961 }
962 
963 
964 
965 template <int dim, int spacedim>
966 std::unique_ptr<FiniteElement<dim, spacedim>>
968 {
969  // Construct a dummy quadrature formula containing the FE's nodes:
970  std::vector<Point<1>> qpoints(this->degree + 1);
971  std::vector<unsigned int> lexicographic =
972  this->poly_space.get_numbering_inverse();
973  for (unsigned int i = 0; i <= this->degree; ++i)
974  qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
975  Quadrature<1> pquadrature(qpoints);
976 
977  return std_cxx14::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(
978  pquadrature);
979 }
980 
981 
982 
983 // ---------------------------------- FE_DGQLegendre -------------------------
984 
985 template <int dim, int spacedim>
987  : FE_DGQ<dim, spacedim>(
988  Polynomials::Legendre::generate_complete_basis(degree))
989 {}
990 
991 
992 
993 template <int dim, int spacedim>
994 std::pair<Table<2, bool>, std::vector<unsigned int>>
996 {
997  // Legendre represents a constant function by one in the first basis
998  // function and zero in all others
999  Table<2, bool> constant_modes(1, this->dofs_per_cell);
1000  constant_modes(0, 0) = true;
1001  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1002  constant_modes, std::vector<unsigned int>(1, 0));
1003 }
1004 
1005 
1006 
1007 template <int dim, int spacedim>
1008 std::string
1010 {
1011  return "FE_DGQLegendre<" + Utilities::dim_string(dim, spacedim) + ">(" +
1012  Utilities::int_to_string(this->degree) + ")";
1013 }
1014 
1015 
1016 
1017 template <int dim, int spacedim>
1018 std::unique_ptr<FiniteElement<dim, spacedim>>
1020 {
1021  return std_cxx14::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1022 }
1023 
1024 
1025 
1026 // ---------------------------------- FE_DGQHermite --------------------------
1027 
1028 template <int dim, int spacedim>
1030  : FE_DGQ<dim, spacedim>(
1031  Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1032 {}
1033 
1034 
1035 
1036 template <int dim, int spacedim>
1037 std::string
1039 {
1040  return "FE_DGQHermite<" + Utilities::dim_string(dim, spacedim) + ">(" +
1041  Utilities::int_to_string(this->degree) + ")";
1042 }
1043 
1044 
1045 
1046 template <int dim, int spacedim>
1047 std::unique_ptr<FiniteElement<dim, spacedim>>
1049 {
1050  return std_cxx14::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1051 }
1052 
1053 
1054 
1055 // explicit instantiations
1056 #include "fe_dgq.inst"
1057 
1058 
1059 DEAL_II_NAMESPACE_CLOSE
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_dgq.cc:995
static ::ExceptionBase & ExcFEHasNoSupportPoints()
size_type m() const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_dgq.cc:596
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_dgq.cc:399
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2442
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
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)
virtual bool hp_constraints_are_implemented() const override
Definition: fe_dgq.cc:559
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:1048
void gauss_jordan()
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_dgq.cc:568
FE_DGQLegendre(const unsigned int degree)
Definition: fe_dgq.cc:986
Definition: fe_dgq.h:109
const unsigned int degree
Definition: fe_base.h:298
virtual std::string get_name() const override
Definition: fe_dgq.cc:849
const Point< dim > & point(const unsigned int i) const
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_dgq.cc:346
STL namespace.
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Definition: point.h:110
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2480
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_dgq.cc:610
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:1019
virtual std::size_t memory_consumption() const override
Definition: fe_dgq.cc:824
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
static ::ExceptionBase & ExcMessage(std::string arg1)
Definition: fe_q.h:554
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
size_type n() const
virtual std::string get_name() const override
Definition: fe_dgq.cc:1009
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2456
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_dgq.cc:812
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_dgq.cc:177
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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
Definition: fe_dgq.cc:143
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix) const override
Definition: fe_dgq.cc:372
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:164
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
Definition: fe_dgq.cc:190
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:383
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:967
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:457
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
unsigned int size() const
const unsigned int dofs_per_cell
Definition: fe_base.h:282
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_dgq.cc:719
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
Definition: fe_dgq.cc:835
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_dgq.cc:267
FE_DGQHermite(const unsigned int degree)
Definition: fe_dgq.cc:1029
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
Definition: fe_dgq.cc:946
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:276
virtual std::string get_name() const override
Definition: fe_dgq.cc:1038
virtual std::string get_name() const override
Definition: fe_dgq.cc:127
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition: fe_dgq.cc:479
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:37
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition: fe_dgq.cc:582
double compute_value(const unsigned int i, const Point< dim > &p) const
void fill(InputIterator entries, const bool C_style_indexing=true)
friend class FE_DGQ
Definition: fe_dgq.h:384
static ::ExceptionBase & ExcInternalError()