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