Reference documentation for deal.II version 9.2.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 - 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 
20 
21 #include <deal.II/fe/fe.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>
27 #include <deal.II/fe/fe_q_dg0.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 
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>(
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),
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 {
403  AssertIndexRange(refinement_case,
405  Assert(refinement_case != RefinementCase<dim>::no_refinement,
406  ExcMessage(
407  "Prolongation matrices are only available for refined cells!"));
408  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
409 
410  // initialization upon first request
411  if (this->prolongation[refinement_case - 1][child].n() == 0)
412  {
413  std::lock_guard<std::mutex> lock(this->mutex);
414 
415  // if matrix got updated while waiting for the lock
416  if (this->prolongation[refinement_case - 1][child].n() ==
417  this->dofs_per_cell)
418  return this->prolongation[refinement_case - 1][child];
419 
420  // now do the work. need to get a non-const version of data in order to
421  // be able to modify them inside a const function
422  FE_DGQ<dim, spacedim> &this_nonconst =
423  const_cast<FE_DGQ<dim, spacedim> &>(*this);
424  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
425  {
426  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
428  isotropic_matrices.back().resize(
431  if (dim == spacedim)
433  isotropic_matrices,
434  true);
435  else
437  isotropic_matrices,
438  true);
439  this_nonconst.prolongation[refinement_case - 1].swap(
440  isotropic_matrices.back());
441  }
442  else
443  {
444  // must compute both restriction and prolongation matrices because
445  // we only check for their size and the reinit call initializes them
446  // all
448  if (dim == spacedim)
449  {
451  this_nonconst.prolongation);
453  this_nonconst.restriction);
454  }
455  else
456  {
457  FE_DGQ<dim> tmp(this->degree);
459  this_nonconst.prolongation);
461  this_nonconst.restriction);
462  }
463  }
464  }
465 
466  // finally return the matrix
467  return this->prolongation[refinement_case - 1][child];
468 }
469 
470 
471 
472 template <int dim, int spacedim>
473 const FullMatrix<double> &
475  const unsigned int child,
476  const RefinementCase<dim> &refinement_case) const
477 {
478  AssertIndexRange(refinement_case,
480  Assert(refinement_case != RefinementCase<dim>::no_refinement,
481  ExcMessage(
482  "Restriction matrices are only available for refined cells!"));
483  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
484 
485  // initialization upon first request
486  if (this->restriction[refinement_case - 1][child].n() == 0)
487  {
488  std::lock_guard<std::mutex> lock(this->mutex);
489 
490  // if matrix got updated while waiting for the lock...
491  if (this->restriction[refinement_case - 1][child].n() ==
492  this->dofs_per_cell)
493  return this->restriction[refinement_case - 1][child];
494 
495  // now do the work. need to get a non-const version of data in order to
496  // be able to modify them inside a const function
497  FE_DGQ<dim, spacedim> &this_nonconst =
498  const_cast<FE_DGQ<dim, spacedim> &>(*this);
499  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
500  {
501  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
503  isotropic_matrices.back().resize(
506  if (dim == spacedim)
508  isotropic_matrices,
509  true);
510  else
512  isotropic_matrices,
513  true);
514  this_nonconst.restriction[refinement_case - 1].swap(
515  isotropic_matrices.back());
516  }
517  else
518  {
519  // must compute both restriction and prolongation matrices because
520  // we only check for their size and the reinit call initializes them
521  // all
523  if (dim == spacedim)
524  {
526  this_nonconst.prolongation);
528  this_nonconst.restriction);
529  }
530  else
531  {
532  FE_DGQ<dim> tmp(this->degree);
534  this_nonconst.prolongation);
536  this_nonconst.restriction);
537  }
538  }
539  }
540 
541  // finally return the matrix
542  return this->restriction[refinement_case - 1][child];
543 }
544 
545 
546 
547 template <int dim, int spacedim>
548 bool
550 {
551  return true;
552 }
553 
554 
555 
556 template <int dim, int spacedim>
557 std::vector<std::pair<unsigned int, unsigned int>>
559  const FiniteElement<dim, spacedim> & /*fe_other*/) const
560 {
561  // this element is discontinuous, so by definition there can
562  // be no identities between its dofs and those of any neighbor
563  // (of whichever type the neighbor may be -- after all, we have
564  // no face dofs on this side to begin with)
565  return std::vector<std::pair<unsigned int, unsigned int>>();
566 }
567 
568 
569 
570 template <int dim, int spacedim>
571 std::vector<std::pair<unsigned int, unsigned int>>
573  const FiniteElement<dim, spacedim> & /*fe_other*/) const
574 {
575  // this element is discontinuous, so by definition there can
576  // be no identities between its dofs and those of any neighbor
577  // (of whichever type the neighbor may be -- after all, we have
578  // no face dofs on this side to begin with)
579  return std::vector<std::pair<unsigned int, unsigned int>>();
580 }
581 
582 
583 
584 template <int dim, int spacedim>
585 std::vector<std::pair<unsigned int, unsigned int>>
587  const FiniteElement<dim, spacedim> & /*fe_other*/) const
588 {
589  // this element is discontinuous, so by definition there can
590  // be no identities between its dofs and those of any neighbor
591  // (of whichever type the neighbor may be -- after all, we have
592  // no face dofs on this side to begin with)
593  return std::vector<std::pair<unsigned int, unsigned int>>();
594 }
595 
596 
597 
598 template <int dim, int spacedim>
601  const FiniteElement<dim, spacedim> &fe_other,
602  const unsigned int codim) const
603 {
604  Assert(codim <= dim, ExcImpossibleInDim(dim));
605 
606  // vertex/line/face domination
607  // ---------------------------
608  if (codim > 0)
609  // this is a discontinuous element, so by definition there will
610  // be no constraints wherever this element comes together with
611  // any other kind of element
613 
614  // cell domination
615  // ---------------
616  // The following block of conditionals is rather ugly, but there is currently
617  // no other way how to deal with a robust comparison of FE_DGQ elements with
618  // relevant others in the current implementation.
619  if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
620  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
621  {
622  if (this->degree < fe_dgq_other->degree)
624  else if (this->degree == fe_dgq_other->degree)
626  else
628  }
629  else if (const FE_Q<dim, spacedim> *fe_q_other =
630  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
631  {
632  if (this->degree < fe_q_other->degree)
634  else if (this->degree == fe_q_other->degree)
636  else
638  }
639  else if (const FE_Bernstein<dim, spacedim> *fe_bernstein_other =
640  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
641  {
642  if (this->degree < fe_bernstein_other->degree)
644  else if (this->degree == fe_bernstein_other->degree)
646  else
648  }
649  else if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
650  dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
651  {
652  if (this->degree < fe_bubbles_other->degree)
654  else if (this->degree == fe_bubbles_other->degree)
656  else
658  }
659  else if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
660  dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
661  {
662  if (this->degree < fe_dg0_other->degree)
664  else if (this->degree == fe_dg0_other->degree)
666  else
668  }
669  else if (const FE_Q_iso_Q1<dim, spacedim> *fe_q_iso_q1_other =
670  dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe_other))
671  {
672  if (this->degree < fe_q_iso_q1_other->degree)
674  else if (this->degree == fe_q_iso_q1_other->degree)
676  else
678  }
679  else if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
680  dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
681  {
682  if (this->degree < fe_hierarchical_other->degree)
684  else if (this->degree == fe_hierarchical_other->degree)
686  else
688  }
689  else if (const FE_Nothing<dim> *fe_nothing =
690  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
691  {
692  if (fe_nothing->is_dominating())
694  else
695  // the FE_Nothing has no degrees of freedom and it is typically used
696  // in a context where we don't require any continuity along the
697  // interface
699  }
700 
701  Assert(false, ExcNotImplemented());
703 }
704 
705 
706 
707 template <int dim, int spacedim>
708 bool
709 FE_DGQ<dim, spacedim>::has_support_on_face(const unsigned int shape_index,
710  const unsigned int face_index) const
711 {
712  AssertIndexRange(shape_index, this->dofs_per_cell);
714 
715  unsigned int n = this->degree + 1;
716 
717  // For DGQ elements that do not define support points, we cannot define
718  // whether they have support at the boundary easily, so return true in any
719  // case
720  if (this->unit_support_points.empty())
721  return true;
722 
723  // for DGQ(0) elements or arbitrary node DGQ with support points not located
724  // at the element boundary, the single shape functions is constant and
725  // therefore lives on the boundary
726  bool support_points_on_boundary = true;
727  for (unsigned int d = 0; d < dim; ++d)
728  if (std::abs(this->unit_support_points[0][d]) > 1e-13)
729  support_points_on_boundary = false;
730  for (unsigned int d = 0; d < dim; ++d)
731  if (std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
732  support_points_on_boundary = false;
733  if (support_points_on_boundary == false)
734  return true;
735 
736  unsigned int n2 = n * n;
737 
738  switch (dim)
739  {
740  case 1:
741  {
742  // in 1d, things are simple. since
743  // there is only one degree of
744  // freedom per vertex in this
745  // class, the first is on vertex 0
746  // (==face 0 in some sense), the
747  // second on face 1:
748  return (((shape_index == 0) && (face_index == 0)) ||
749  ((shape_index == this->degree) && (face_index == 1)));
750  }
751 
752  case 2:
753  {
754  if (face_index == 0 && (shape_index % n) == 0)
755  return true;
756  if (face_index == 1 && (shape_index % n) == this->degree)
757  return true;
758  if (face_index == 2 && shape_index < n)
759  return true;
760  if (face_index == 3 && shape_index >= this->dofs_per_cell - n)
761  return true;
762  return false;
763  }
764 
765  case 3:
766  {
767  const unsigned int in2 = shape_index % n2;
768 
769  // x=0
770  if (face_index == 0 && (shape_index % n) == 0)
771  return true;
772  // x=1
773  if (face_index == 1 && (shape_index % n) == n - 1)
774  return true;
775  // y=0
776  if (face_index == 2 && in2 < n)
777  return true;
778  // y=1
779  if (face_index == 3 && in2 >= n2 - n)
780  return true;
781  // z=0
782  if (face_index == 4 && shape_index < n2)
783  return true;
784  // z=1
785  if (face_index == 5 && shape_index >= this->dofs_per_cell - n2)
786  return true;
787  return false;
788  }
789 
790  default:
791  Assert(false, ExcNotImplemented());
792  }
793  return true;
794 }
795 
796 
797 
798 template <int dim, int spacedim>
799 std::pair<Table<2, bool>, std::vector<unsigned int>>
801 {
802  Table<2, bool> constant_modes(1, this->dofs_per_cell);
803  constant_modes.fill(true);
804  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
805  constant_modes, std::vector<unsigned int>(1, 0));
806 }
807 
808 
809 
810 // ------------------------------ FE_DGQArbitraryNodes -----------------------
811 
812 template <int dim, int spacedim>
814  const Quadrature<1> &points)
815  : FE_DGQ<dim, spacedim>(
816  Polynomials::generate_complete_Lagrange_basis(points.get_points()))
817 {
818  Assert(points.size() > 0,
820  this->unit_support_points = Quadrature<dim>(points).get_points();
821 }
822 
823 
824 
825 template <int dim, int spacedim>
826 std::string
828 {
829  // note that the FETools::get_fe_by_name function does not work for
830  // FE_DGQArbitraryNodes since there is no initialization by a degree value.
831  std::ostringstream namebuf;
832  bool equidistant = true;
833  std::vector<double> points(this->degree + 1);
834 
835  std::vector<unsigned int> lexicographic =
837  for (unsigned int j = 0; j <= this->degree; j++)
838  points[j] = this->unit_support_points[lexicographic[j]][0];
839 
840  // Check whether the support points are equidistant.
841  for (unsigned int j = 0; j <= this->degree; j++)
842  if (std::abs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
843  {
844  equidistant = false;
845  break;
846  }
847  if (this->degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
848  equidistant = true;
849 
850  if (equidistant == true)
851  {
852  if (this->degree > 2)
853  namebuf << "FE_DGQArbitraryNodes<"
854  << Utilities::dim_string(dim, spacedim)
855  << ">(QIterated(QTrapez()," << this->degree << "))";
856  else
857  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
858  << this->degree << ")";
859  return namebuf.str();
860  }
861 
862  // Check whether the support points come from QGaussLobatto.
863  const QGaussLobatto<1> points_gl(this->degree + 1);
864  bool gauss_lobatto = true;
865  for (unsigned int j = 0; j <= this->degree; j++)
866  if (points[j] != points_gl.point(j)(0))
867  {
868  gauss_lobatto = false;
869  break;
870  }
871 
872  if (gauss_lobatto == true)
873  {
874  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
875  << this->degree << ")";
876  return namebuf.str();
877  }
878 
879  // Check whether the support points come from QGauss.
880  const QGauss<1> points_g(this->degree + 1);
881  bool gauss = true;
882  for (unsigned int j = 0; j <= this->degree; j++)
883  if (points[j] != points_g.point(j)(0))
884  {
885  gauss = false;
886  break;
887  }
888 
889  if (gauss == true)
890  {
891  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
892  << ">(QGauss(" << this->degree + 1 << "))";
893  return namebuf.str();
894  }
895 
896  // Check whether the support points come from QGauss.
897  const QGaussLog<1> points_glog(this->degree + 1);
898  bool gauss_log = true;
899  for (unsigned int j = 0; j <= this->degree; j++)
900  if (points[j] != points_glog.point(j)(0))
901  {
902  gauss_log = false;
903  break;
904  }
905 
906  if (gauss_log == true)
907  {
908  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
909  << ">(QGaussLog(" << this->degree + 1 << "))";
910  return namebuf.str();
911  }
912 
913  // All guesses exhausted
914  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
915  << ">(QUnknownNodes(" << this->degree + 1 << "))";
916  return namebuf.str();
917 }
918 
919 
920 
921 template <int dim, int spacedim>
922 void
925  const std::vector<Vector<double>> &support_point_values,
926  std::vector<double> & nodal_values) const
927 {
928  AssertDimension(support_point_values.size(),
929  this->get_unit_support_points().size());
930  AssertDimension(support_point_values.size(), nodal_values.size());
931  AssertDimension(this->dofs_per_cell, nodal_values.size());
932 
933  for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
934  {
935  AssertDimension(support_point_values[i].size(), 1);
936 
937  nodal_values[i] = support_point_values[i](0);
938  }
939 }
940 
941 
942 
943 template <int dim, int spacedim>
944 std::unique_ptr<FiniteElement<dim, spacedim>>
946 {
947  // Construct a dummy quadrature formula containing the FE's nodes:
948  std::vector<Point<1>> qpoints(this->degree + 1);
949  std::vector<unsigned int> lexicographic =
951  for (unsigned int i = 0; i <= this->degree; ++i)
952  qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
953  Quadrature<1> pquadrature(qpoints);
954 
955  return std_cxx14::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(
956  pquadrature);
957 }
958 
959 
960 
961 // ---------------------------------- FE_DGQLegendre -------------------------
962 
963 template <int dim, int spacedim>
965  : FE_DGQ<dim, spacedim>(
966  Polynomials::Legendre::generate_complete_basis(degree))
967 {}
968 
969 
970 
971 template <int dim, int spacedim>
972 std::pair<Table<2, bool>, std::vector<unsigned int>>
974 {
975  // Legendre represents a constant function by one in the first basis
976  // function and zero in all others
977  Table<2, bool> constant_modes(1, this->dofs_per_cell);
978  constant_modes(0, 0) = true;
979  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
980  constant_modes, std::vector<unsigned int>(1, 0));
981 }
982 
983 
984 
985 template <int dim, int spacedim>
986 std::string
988 {
989  return "FE_DGQLegendre<" + Utilities::dim_string(dim, spacedim) + ">(" +
990  Utilities::int_to_string(this->degree) + ")";
991 }
992 
993 
994 
995 template <int dim, int spacedim>
996 std::unique_ptr<FiniteElement<dim, spacedim>>
998 {
999  return std_cxx14::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1000 }
1001 
1002 
1003 
1004 // ---------------------------------- FE_DGQHermite --------------------------
1005 
1006 template <int dim, int spacedim>
1008  : FE_DGQ<dim, spacedim>(
1009  Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1010 {}
1011 
1012 
1013 
1014 template <int dim, int spacedim>
1015 std::string
1017 {
1018  return "FE_DGQHermite<" + Utilities::dim_string(dim, spacedim) + ">(" +
1019  Utilities::int_to_string(this->degree) + ")";
1020 }
1021 
1022 
1023 
1024 template <int dim, int spacedim>
1025 std::unique_ptr<FiniteElement<dim, spacedim>>
1027 {
1028  return std_cxx14::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1029 }
1030 
1031 
1032 
1033 // explicit instantiations
1034 #include "fe_dgq.inst"
1035 
1036 
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_dgq.cc:973
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:586
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
Threads::Mutex mutex
Definition: fe_dgq.h:370
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2414
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
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:549
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:1026
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:558
FE_DGQLegendre(const unsigned int degree)
Definition: fe_dgq.cc:964
Definition: fe_dgq.h:112
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const unsigned int degree
Definition: fe_base.h:298
virtual std::string get_name() const override
Definition: fe_dgq.cc:827
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.
static const char U
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
Definition: point.h:111
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2452
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_dgq.cc:600
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:997
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:987
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2428
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:800
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_dgq.cc:177
#define Assert(cond, exc)
Definition: exceptions.h:1419
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
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
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
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:945
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:559
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:709
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
Definition: fe_dgq.cc:813
const std::vector< Point< dim > > & get_unit_support_points() const
Definition: fe.cc:991
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:1007
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
const std::vector< unsigned int > & get_numbering_inverse() const
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:924
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
Definition: fe.cc:275
virtual std::string get_name() const override
Definition: fe_dgq.cc:1016
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:474
static ::ExceptionBase & ExcNotImplemented()
Definition: table.h:699
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:572
double compute_value(const unsigned int i, const Point< dim > &p) const
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:171
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:374
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()