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