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