Reference documentation for deal.II version GIT 6bdf9a4b6c 2022-08-12 19:20:02+00:00
\(\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),
108  std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(
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 #ifdef DEBUG
334  // make sure that the row sum of
335  // each of the matrices is 1 at
336  // this point. this must be so
337  // since the shape functions sum up
338  // to 1
339  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
340  {
341  double sum = 0.;
342  for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
343  sum += interpolation_matrix(i, j);
344 
345  Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
346  ExcInternalError());
347  }
348 #endif
349 }
350 
351 
352 
353 template <int dim, int spacedim>
354 void
356  const FiniteElement<dim, spacedim> &x_source_fe,
357  FullMatrix<double> & interpolation_matrix,
358  const unsigned int) const
359 {
360  // this is only implemented, if the source
361  // FE is also a DGQ element. in that case,
362  // both elements have no dofs on their
363  // faces and the face interpolation matrix
364  // is necessarily empty -- i.e. there isn't
365  // much we need to do here.
366  (void)interpolation_matrix;
368  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
369  nullptr),
370  typename FE::ExcInterpolationNotImplemented());
371 
372  Assert(interpolation_matrix.m() == 0,
373  ExcDimensionMismatch(interpolation_matrix.m(), 0));
374  Assert(interpolation_matrix.n() == 0,
375  ExcDimensionMismatch(interpolation_matrix.m(), 0));
376 }
377 
378 
379 
380 template <int dim, int spacedim>
381 void
383  const FiniteElement<dim, spacedim> &x_source_fe,
384  const unsigned int,
385  FullMatrix<double> &interpolation_matrix,
386  const unsigned int) const
387 {
388  // this is only implemented, if the source
389  // FE is also a DGQ element. in that case,
390  // both elements have no dofs on their
391  // faces and the face interpolation matrix
392  // is necessarily empty -- i.e. there isn't
393  // much we need to do here.
394  (void)interpolation_matrix;
395  using FE = FiniteElement<dim, spacedim>;
396  AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
397  nullptr),
398  typename FE::ExcInterpolationNotImplemented());
399 
400  Assert(interpolation_matrix.m() == 0,
401  ExcDimensionMismatch(interpolation_matrix.m(), 0));
402  Assert(interpolation_matrix.n() == 0,
403  ExcDimensionMismatch(interpolation_matrix.m(), 0));
404 }
405 
406 
407 
408 template <int dim, int spacedim>
409 const FullMatrix<double> &
411  const unsigned int child,
412  const RefinementCase<dim> &refinement_case) const
413 {
414  AssertIndexRange(refinement_case,
416  Assert(refinement_case != RefinementCase<dim>::no_refinement,
417  ExcMessage(
418  "Prolongation matrices are only available for refined cells!"));
419  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
420 
421  // initialization upon first request
422  if (this->prolongation[refinement_case - 1][child].n() == 0)
423  {
424  std::lock_guard<std::mutex> lock(this->mutex);
425 
426  // if matrix got updated while waiting for the lock
427  if (this->prolongation[refinement_case - 1][child].n() ==
428  this->n_dofs_per_cell())
429  return this->prolongation[refinement_case - 1][child];
430 
431  // now do the work. need to get a non-const version of data in order to
432  // be able to modify them inside a const function
433  FE_DGQ<dim, spacedim> &this_nonconst =
434  const_cast<FE_DGQ<dim, spacedim> &>(*this);
435  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
436  {
437  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
439  isotropic_matrices.back().resize(
441  FullMatrix<double>(this->n_dofs_per_cell(),
442  this->n_dofs_per_cell()));
443  if (dim == spacedim)
445  isotropic_matrices,
446  true);
447  else
449  isotropic_matrices,
450  true);
451  this_nonconst.prolongation[refinement_case - 1].swap(
452  isotropic_matrices.back());
453  }
454  else
455  {
456  // must compute both restriction and prolongation matrices because
457  // we only check for their size and the reinit call initializes them
458  // all
460  if (dim == spacedim)
461  {
463  this_nonconst.prolongation);
465  this_nonconst.restriction);
466  }
467  else
468  {
469  FE_DGQ<dim> tmp(this->degree);
471  this_nonconst.prolongation);
473  this_nonconst.restriction);
474  }
475  }
476  }
477 
478  // finally return the matrix
479  return this->prolongation[refinement_case - 1][child];
480 }
481 
482 
483 
484 template <int dim, int spacedim>
485 const FullMatrix<double> &
487  const unsigned int child,
488  const RefinementCase<dim> &refinement_case) const
489 {
490  AssertIndexRange(refinement_case,
492  Assert(refinement_case != RefinementCase<dim>::no_refinement,
493  ExcMessage(
494  "Restriction matrices are only available for refined cells!"));
495  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
496 
497  // initialization upon first request
498  if (this->restriction[refinement_case - 1][child].n() == 0)
499  {
500  std::lock_guard<std::mutex> lock(this->mutex);
501 
502  // if matrix got updated while waiting for the lock...
503  if (this->restriction[refinement_case - 1][child].n() ==
504  this->n_dofs_per_cell())
505  return this->restriction[refinement_case - 1][child];
506 
507  // now do the work. need to get a non-const version of data in order to
508  // be able to modify them inside a const function
509  FE_DGQ<dim, spacedim> &this_nonconst =
510  const_cast<FE_DGQ<dim, spacedim> &>(*this);
511  if (refinement_case == RefinementCase<dim>::isotropic_refinement)
512  {
513  std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
515  isotropic_matrices.back().resize(
517  FullMatrix<double>(this->n_dofs_per_cell(),
518  this->n_dofs_per_cell()));
519  if (dim == spacedim)
521  isotropic_matrices,
522  true);
523  else
525  isotropic_matrices,
526  true);
527  this_nonconst.restriction[refinement_case - 1].swap(
528  isotropic_matrices.back());
529  }
530  else
531  {
532  // must compute both restriction and prolongation matrices because
533  // we only check for their size and the reinit call initializes them
534  // all
536  if (dim == spacedim)
537  {
539  this_nonconst.prolongation);
541  this_nonconst.restriction);
542  }
543  else
544  {
545  FE_DGQ<dim> tmp(this->degree);
547  this_nonconst.prolongation);
549  this_nonconst.restriction);
550  }
551  }
552  }
553 
554  // finally return the matrix
555  return this->restriction[refinement_case - 1][child];
556 }
557 
558 
559 
560 template <int dim, int spacedim>
561 bool
563 {
564  return true;
565 }
566 
567 
568 
569 template <int dim, int spacedim>
570 std::vector<std::pair<unsigned int, unsigned int>>
572  const FiniteElement<dim, spacedim> & /*fe_other*/) const
573 {
574  // this element is discontinuous, so by definition there can
575  // be no identities between its dofs and those of any neighbor
576  // (of whichever type the neighbor may be -- after all, we have
577  // no face dofs on this side to begin with)
578  return std::vector<std::pair<unsigned int, unsigned int>>();
579 }
580 
581 
582 
583 template <int dim, int spacedim>
584 std::vector<std::pair<unsigned int, unsigned int>>
586  const FiniteElement<dim, spacedim> & /*fe_other*/) const
587 {
588  // this element is discontinuous, so by definition there can
589  // be no identities between its dofs and those of any neighbor
590  // (of whichever type the neighbor may be -- after all, we have
591  // no face dofs on this side to begin with)
592  return std::vector<std::pair<unsigned int, unsigned int>>();
593 }
594 
595 
596 
597 template <int dim, int spacedim>
598 std::vector<std::pair<unsigned int, unsigned int>>
600  const FiniteElement<dim, spacedim> & /*fe_other*/,
601  const unsigned int) const
602 {
603  // this element is discontinuous, so by definition there can
604  // be no identities between its dofs and those of any neighbor
605  // (of whichever type the neighbor may be -- after all, we have
606  // no face dofs on this side to begin with)
607  return std::vector<std::pair<unsigned int, unsigned int>>();
608 }
609 
610 
611 
612 template <int dim, int spacedim>
615  const FiniteElement<dim, spacedim> &fe_other,
616  const unsigned int codim) const
617 {
618  Assert(codim <= dim, ExcImpossibleInDim(dim));
619 
620  // vertex/line/face domination
621  // ---------------------------
622  if (codim > 0)
623  // this is a discontinuous element, so by definition there will
624  // be no constraints wherever this element comes together with
625  // any other kind of element
627 
628  // cell domination
629  // ---------------
630  // The following block of conditionals is rather ugly, but there is currently
631  // no other way how to deal with a robust comparison of FE_DGQ elements with
632  // relevant others in the current implementation.
633  if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
634  dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
635  {
636  if (this->degree < fe_dgq_other->degree)
638  else if (this->degree == fe_dgq_other->degree)
640  else
642  }
643  else if (const FE_Q<dim, spacedim> *fe_q_other =
644  dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
645  {
646  if (this->degree < fe_q_other->degree)
648  else if (this->degree == fe_q_other->degree)
650  else
652  }
653  else if (const FE_Bernstein<dim, spacedim> *fe_bernstein_other =
654  dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
655  {
656  if (this->degree < fe_bernstein_other->degree)
658  else if (this->degree == fe_bernstein_other->degree)
660  else
662  }
663  else if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
664  dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
665  {
666  if (this->degree < fe_bubbles_other->degree)
668  else if (this->degree == fe_bubbles_other->degree)
670  else
672  }
673  else if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
674  dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
675  {
676  if (this->degree < fe_dg0_other->degree)
678  else if (this->degree == fe_dg0_other->degree)
680  else
682  }
683  else if (const FE_Q_iso_Q1<dim, spacedim> *fe_q_iso_q1_other =
684  dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe_other))
685  {
686  if (this->degree < fe_q_iso_q1_other->degree)
688  else if (this->degree == fe_q_iso_q1_other->degree)
690  else
692  }
693  else if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
694  dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
695  {
696  if (this->degree < fe_hierarchical_other->degree)
698  else if (this->degree == fe_hierarchical_other->degree)
700  else
702  }
703  else if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
704  dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
705  {
706  if (this->degree < fe_dgp_other->degree)
708  else if (this->degree == fe_dgp_other->degree)
710  else
712  }
713  else if (const FE_Nothing<dim> *fe_nothing =
714  dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
715  {
716  if (fe_nothing->is_dominating())
718  else
719  // the FE_Nothing has no degrees of freedom and it is typically used
720  // in a context where we don't require any continuity along the
721  // interface
723  }
724 
725  Assert(false, ExcNotImplemented());
727 }
728 
729 
730 
731 template <int dim, int spacedim>
732 bool
733 FE_DGQ<dim, spacedim>::has_support_on_face(const unsigned int shape_index,
734  const unsigned int face_index) const
735 {
736  AssertIndexRange(shape_index, this->n_dofs_per_cell());
738 
739  unsigned int n = this->degree + 1;
740 
741  // For DGQ elements that do not define support points, we cannot define
742  // whether they have support at the boundary easily, so return true in any
743  // case
744  if (this->unit_support_points.empty())
745  return true;
746 
747  // for DGQ(0) elements or arbitrary node DGQ with support points not located
748  // at the element boundary, the single shape functions is constant and
749  // therefore lives on the boundary
750  bool support_points_on_boundary = true;
751  for (unsigned int d = 0; d < dim; ++d)
752  if (std::abs(this->unit_support_points[0][d]) > 1e-13)
753  support_points_on_boundary = false;
754  for (unsigned int d = 0; d < dim; ++d)
755  if (std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
756  support_points_on_boundary = false;
757  if (support_points_on_boundary == false)
758  return true;
759 
760  unsigned int n2 = n * n;
761 
762  switch (dim)
763  {
764  case 1:
765  {
766  // in 1d, things are simple. since
767  // there is only one degree of
768  // freedom per vertex in this
769  // class, the first is on vertex 0
770  // (==face 0 in some sense), the
771  // second on face 1:
772  return (((shape_index == 0) && (face_index == 0)) ||
773  ((shape_index == this->degree) && (face_index == 1)));
774  }
775 
776  case 2:
777  {
778  if (face_index == 0 && (shape_index % n) == 0)
779  return true;
780  if (face_index == 1 && (shape_index % n) == this->degree)
781  return true;
782  if (face_index == 2 && shape_index < n)
783  return true;
784  if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
785  return true;
786  return false;
787  }
788 
789  case 3:
790  {
791  const unsigned int in2 = shape_index % n2;
792 
793  // x=0
794  if (face_index == 0 && (shape_index % n) == 0)
795  return true;
796  // x=1
797  if (face_index == 1 && (shape_index % n) == n - 1)
798  return true;
799  // y=0
800  if (face_index == 2 && in2 < n)
801  return true;
802  // y=1
803  if (face_index == 3 && in2 >= n2 - n)
804  return true;
805  // z=0
806  if (face_index == 4 && shape_index < n2)
807  return true;
808  // z=1
809  if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
810  return true;
811  return false;
812  }
813 
814  default:
815  Assert(false, ExcNotImplemented());
816  }
817  return true;
818 }
819 
820 
821 
822 template <int dim, int spacedim>
823 std::pair<Table<2, bool>, std::vector<unsigned int>>
825 {
826  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
827  constant_modes.fill(true);
828  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
829  constant_modes, std::vector<unsigned int>(1, 0));
830 }
831 
832 
833 
834 // ------------------------------ FE_DGQArbitraryNodes -----------------------
835 
836 template <int dim, int spacedim>
838  const Quadrature<1> &points)
839  : FE_DGQ<dim, spacedim>(
840  Polynomials::generate_complete_Lagrange_basis(points.get_points()))
841 {
842  Assert(points.size() > 0,
845 }
846 
847 
848 
849 template <int dim, int spacedim>
850 std::string
852 {
853  // note that the FETools::get_fe_by_name function does not work for
854  // FE_DGQArbitraryNodes since there is no initialization by a degree value.
855  std::ostringstream namebuf;
856  bool equidistant = true;
857  std::vector<double> points(this->degree + 1);
858 
859  auto *const polynomial_space =
860  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
861  Assert(polynomial_space != nullptr, ExcInternalError());
862  std::vector<unsigned int> lexicographic =
863  polynomial_space->get_numbering_inverse();
864  for (unsigned int j = 0; j <= this->degree; ++j)
865  points[j] = this->unit_support_points[lexicographic[j]][0];
866 
867  // Check whether the support points are equidistant.
868  for (unsigned int j = 0; j <= this->degree; ++j)
869  if (std::abs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
870  {
871  equidistant = false;
872  break;
873  }
874  if (this->degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
875  equidistant = true;
876 
877  if (equidistant == true)
878  {
879  if (this->degree > 2)
880  namebuf << "FE_DGQArbitraryNodes<"
881  << Utilities::dim_string(dim, spacedim)
882  << ">(QIterated(QTrapezoid()," << this->degree << "))";
883  else
884  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
885  << this->degree << ")";
886  return namebuf.str();
887  }
888 
889  // Check whether the support points come from QGaussLobatto.
890  const QGaussLobatto<1> points_gl(this->degree + 1);
891  bool gauss_lobatto = true;
892  for (unsigned int j = 0; j <= this->degree; ++j)
893  if (points[j] != points_gl.point(j)(0))
894  {
895  gauss_lobatto = false;
896  break;
897  }
898 
899  if (gauss_lobatto == true)
900  {
901  namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
902  << this->degree << ")";
903  return namebuf.str();
904  }
905 
906  // Check whether the support points come from QGauss.
907  const QGauss<1> points_g(this->degree + 1);
908  bool gauss = true;
909  for (unsigned int j = 0; j <= this->degree; ++j)
910  if (points[j] != points_g.point(j)(0))
911  {
912  gauss = false;
913  break;
914  }
915 
916  if (gauss == true)
917  {
918  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
919  << ">(QGauss(" << this->degree + 1 << "))";
920  return namebuf.str();
921  }
922 
923  // Check whether the support points come from QGauss.
924  const QGaussLog<1> points_glog(this->degree + 1);
925  bool gauss_log = true;
926  for (unsigned int j = 0; j <= this->degree; ++j)
927  if (points[j] != points_glog.point(j)(0))
928  {
929  gauss_log = false;
930  break;
931  }
932 
933  if (gauss_log == true)
934  {
935  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
936  << ">(QGaussLog(" << this->degree + 1 << "))";
937  return namebuf.str();
938  }
939 
940  // All guesses exhausted
941  namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
942  << ">(QUnknownNodes(" << this->degree + 1 << "))";
943  return namebuf.str();
944 }
945 
946 
947 
948 template <int dim, int spacedim>
949 void
952  const std::vector<Vector<double>> &support_point_values,
953  std::vector<double> & nodal_values) const
954 {
955  AssertDimension(support_point_values.size(),
956  this->get_unit_support_points().size());
957  AssertDimension(support_point_values.size(), nodal_values.size());
958  AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
959 
960  for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
961  {
962  AssertDimension(support_point_values[i].size(), 1);
963 
964  nodal_values[i] = support_point_values[i](0);
965  }
966 }
967 
968 
969 
970 template <int dim, int spacedim>
971 std::unique_ptr<FiniteElement<dim, spacedim>>
973 {
974  // Construct a dummy quadrature formula containing the FE's nodes:
975  std::vector<Point<1>> qpoints(this->degree + 1);
976  auto *const polynomial_space =
977  dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
978  Assert(polynomial_space != nullptr, ExcInternalError());
979  std::vector<unsigned int> lexicographic =
980  polynomial_space->get_numbering_inverse();
981  for (unsigned int i = 0; i <= this->degree; ++i)
982  qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
983  Quadrature<1> pquadrature(qpoints);
984 
985  return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
986 }
987 
988 
989 
990 // ---------------------------------- FE_DGQLegendre -------------------------
991 
992 template <int dim, int spacedim>
994  : FE_DGQ<dim, spacedim>(
995  Polynomials::Legendre::generate_complete_basis(degree))
996 {}
997 
998 
999 
1000 template <int dim, int spacedim>
1001 std::pair<Table<2, bool>, std::vector<unsigned int>>
1003 {
1004  // Legendre represents a constant function by one in the first basis
1005  // function and zero in all others
1006  Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1007  constant_modes(0, 0) = true;
1008  return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1009  constant_modes, std::vector<unsigned int>(1, 0));
1010 }
1011 
1012 
1013 
1014 template <int dim, int spacedim>
1015 std::string
1017 {
1018  return "FE_DGQLegendre<" + 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::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1029 }
1030 
1031 
1032 
1033 // ---------------------------------- FE_DGQHermite --------------------------
1034 
1035 template <int dim, int spacedim>
1037  : FE_DGQ<dim, spacedim>(
1038  Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1039 {}
1040 
1041 
1042 
1043 template <int dim, int spacedim>
1044 std::string
1046 {
1047  return "FE_DGQHermite<" + Utilities::dim_string(dim, spacedim) + ">(" +
1048  Utilities::int_to_string(this->degree) + ")";
1049 }
1050 
1051 
1052 
1053 template <int dim, int spacedim>
1054 std::unique_ptr<FiniteElement<dim, spacedim>>
1056 {
1057  return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1058 }
1059 
1060 
1061 
1062 // explicit instantiations
1063 #include "fe_dgq.inst"
1064 
1065 
virtual std::string get_name() const override
Definition: fe_dgq.cc:851
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
Definition: fe_dgq.cc:837
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:951
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:972
virtual std::string get_name() const override
Definition: fe_dgq.cc:1045
FE_DGQHermite(const unsigned int degree)
Definition: fe_dgq.cc:1036
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:1055
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_dgq.cc:1002
FE_DGQLegendre(const unsigned int degree)
Definition: fe_dgq.cc:993
virtual std::string get_name() const override
Definition: fe_dgq.cc:1016
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:1026
Definition: fe_dgq.h:113
friend class FE_DGQ
Definition: fe_dgq.h:377
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition: fe_dgq.cc:733
virtual bool hp_constraints_are_implemented() const override
Definition: fe_dgq.cc:562
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
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:599
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition: fe_dgq.cc:183
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition: fe_dgq.cc:824
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:486
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:355
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition: fe_dgq.cc:614
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:571
virtual std::string get_name() const override
Definition: fe_dgq.cc:133
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
Definition: fe_dgq.cc:196
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:585
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition: fe_dgq.cc:170
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:410
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:382
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition: fe_dgq.cc:273
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Definition: fe_poly.h:536
Definition: fe_q.h:551
const unsigned int degree
Definition: fe_data.h:449
unsigned int n_dofs_per_cell() const
std::vector< std::vector< FullMatrix< double > > > restriction
Definition: fe.h:2386
void reinit_restriction_and_prolongation_matrices(const bool isotropic_restriction_only=false, const bool isotropic_prolongation_only=false)
std::vector< Point< dim > > unit_support_points
Definition: fe.h:2424
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition: fe.h:2400
size_type n() const
size_type m() const
const std::vector< Point< dim > > & get_points() const
const Point< dim > & point(const unsigned int i) const
unsigned int size() const
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
Definition: vector.h:110
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
Expression fabs(const Expression &x)
void compute_embedding_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false, const double threshold=1.e-12)
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number >>> &matrices, const bool isotropic_only=false)
std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
static const char U
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:160
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< Polynomial< double > > generate_complete_Lagrange_basis(const std::vector< Point< 1 >> &points)
Definition: polynomial.cc:702
std::string dim_string(const int dim, const int spacedim)
Definition: utilities.cc:558
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
std::vector< Point< dim > > unit_support_points(const std::vector< Point< 1 >> &line_support_points, const std::vector< unsigned int > &renumbering)