deal.II version GIT relicensing-2846-g6fb608615f 2025-03-15 04:10:00+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: LGPL-2.1-or-later
4// Copyright (C) 2001 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18
19#include <deal.II/fe/fe.h>
21#include <deal.II/fe/fe_dgq.h>
24#include <deal.II/fe/fe_q.h>
26#include <deal.II/fe/fe_q_dg0.h>
31#include <deal.II/fe/fe_tools.h>
33
34#include <deal.II/lac/vector.h>
35
36#include <iostream>
37#include <memory>
38#include <sstream>
39
40
42
43
44namespace internal
45{
46 namespace FE_DGQ
47 {
48 namespace
49 {
50 std::vector<Point<1>>
51 get_QGaussLobatto_points(const unsigned int degree)
52 {
53 if (degree > 0)
54 return QGaussLobatto<1>(degree + 1).get_points();
55 else
56 return std::vector<Point<1>>(1, Point<1>(0.5));
57 }
58 } // namespace
59 } // namespace FE_DGQ
60} // namespace internal
61
62
63
64template <int dim, int spacedim>
65FE_DGQ<dim, spacedim>::FE_DGQ(const unsigned int degree)
66 : FE_Poly<dim, spacedim>(
68 Polynomials::generate_complete_Lagrange_basis(
69 internal::FE_DGQ::get_QGaussLobatto_points(degree))),
70 FiniteElementData<dim>(get_dpo_vector(degree),
71 1,
72 degree,
73 FiniteElementData<dim>::L2),
74 std::vector<bool>(
75 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
76 .n_dofs_per_cell(),
77 true),
78 std::vector<ComponentMask>(
79 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
80 .n_dofs_per_cell(),
81 ComponentMask(std::vector<bool>(1, true))))
82{
83 // Compute support points, which are the tensor product of the Lagrange
84 // interpolation points in the constructor.
86 Quadrature<dim>(internal::FE_DGQ::get_QGaussLobatto_points(degree))
87 .get_points();
88
89 // do not initialize embedding and restriction here. these matrices are
90 // initialized on demand in get_restriction_matrix and
91 // get_prolongation_matrix
92
93 // note: no face support points for DG elements
94}
95
96
97
98template <int dim, int spacedim>
100 const std::vector<Polynomials::Polynomial<double>> &polynomials)
101 : FE_Poly<dim, spacedim>(
102 TensorProductPolynomials<dim>(polynomials),
103 FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
104 1,
105 polynomials.size() - 1,
106 FiniteElementData<dim>::L2),
107 std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(
108 polynomials.size() - 1),
109 1,
110 polynomials.size() - 1)
111 .n_dofs_per_cell(),
112 true),
113 std::vector<ComponentMask>(
114 FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
115 1,
116 polynomials.size() - 1)
117 .n_dofs_per_cell(),
118 ComponentMask(std::vector<bool>(1, true))))
119{
120 // No support points can be defined in general. Derived classes might define
121 // support points like the class FE_DGQArbitraryNodes
122
123 // do not initialize embedding and restriction here. these matrices are
124 // initialized on demand in get_restriction_matrix and
125 // get_prolongation_matrix
126}
127
129
130template <int dim, int spacedim>
131std::string
133{
134 // note that the FETools::get_fe_by_name function depends on the
135 // particular format of the string this function returns, so they have to be
136 // kept in sync
137
138 std::ostringstream namebuf;
139 namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
140 << this->degree << ")";
141 return namebuf.str();
142}
143
144
145
146template <int dim, int spacedim>
147void
149 const std::vector<Vector<double>> &support_point_values,
150 std::vector<double> &nodal_values) const
151{
152 AssertDimension(support_point_values.size(),
153 this->get_unit_support_points().size());
154 AssertDimension(support_point_values.size(), nodal_values.size());
155 AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
156
157 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
158 {
159 AssertDimension(support_point_values[i].size(), 1);
160
161 nodal_values[i] = support_point_values[i](0);
162 }
163}
164
165
166
167template <int dim, int spacedim>
168std::unique_ptr<FiniteElement<dim, spacedim>>
170{
171 return std::make_unique<FE_DGQ<dim, spacedim>>(*this);
172}
173
174
175//---------------------------------------------------------------------------
176// Auxiliary functions
177//---------------------------------------------------------------------------
178
179
180template <int dim, int spacedim>
181std::vector<unsigned int>
183{
184 std::vector<unsigned int> dpo(dim + 1, 0U);
185 dpo[dim] = deg + 1;
186 for (unsigned int i = 1; i < dim; ++i)
187 dpo[dim] *= deg + 1;
188 return dpo;
189}
190
191
192
193template <int dim, int spacedim>
194void
196 const char direction) const
197{
198 const unsigned int n = this->degree + 1;
199 unsigned int s = n;
200 for (unsigned int i = 1; i < dim; ++i)
201 s *= n;
202 numbers.resize(s);
203
204 unsigned int l = 0;
205
206 if (dim == 1)
207 {
208 // Mirror around midpoint
209 for (unsigned int i = n; i > 0;)
210 numbers[l++] = --i;
211 }
212 else
213 {
214 switch (direction)
215 {
216 // Rotate xy-plane
217 // counter-clockwise
218 case 'z':
219 for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
220 for (unsigned int j = 0; j < n; ++j)
221 for (unsigned int i = 0; i < n; ++i)
223 unsigned int k = n * i - j + n - 1 + n * n * iz;
224 numbers[l++] = k;
225 }
226 break;
227 // Rotate xy-plane
228 // clockwise
229 case 'Z':
230 for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
231 for (unsigned int iy = 0; iy < n; ++iy)
232 for (unsigned int ix = 0; ix < n; ++ix)
233 {
234 unsigned int k = n * ix - iy + n - 1 + n * n * iz;
235 numbers[k] = l++;
236 }
237 break;
238 // Rotate yz-plane
239 // counter-clockwise
240 case 'x':
241 Assert(dim > 2, ExcDimensionMismatch(dim, 3));
242 for (unsigned int iz = 0; iz < n; ++iz)
243 for (unsigned int iy = 0; iy < n; ++iy)
244 for (unsigned int ix = 0; ix < n; ++ix)
245 {
246 unsigned int k = n * (n * iy - iz + n - 1) + ix;
247 numbers[l++] = k;
248 }
249 break;
250 // Rotate yz-plane
251 // clockwise
252 case 'X':
253 Assert(dim > 2, ExcDimensionMismatch(dim, 3));
254 for (unsigned int iz = 0; iz < n; ++iz)
255 for (unsigned int iy = 0; iy < n; ++iy)
256 for (unsigned int ix = 0; ix < n; ++ix)
257 {
258 unsigned int k = n * (n * iy - iz + n - 1) + ix;
259 numbers[k] = l++;
260 }
261 break;
262 default:
264 }
265 }
266}
267
268
269
270template <int dim, int spacedim>
271void
274 FullMatrix<double> &interpolation_matrix) const
275{
276 // go through the list of elements we can interpolate from
277 if (const FE_DGQ<dim, spacedim> *source_fe =
278 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe))
279 {
280 Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
281 ExcDimensionMismatch(interpolation_matrix.m(),
282 this->n_dofs_per_cell()));
283 Assert(interpolation_matrix.n() == source_fe->n_dofs_per_cell(),
284 ExcDimensionMismatch(interpolation_matrix.n(),
285 source_fe->n_dofs_per_cell()));
286
287
288 // compute the interpolation
289 // matrices in much the same way as
290 // we do for the embedding matrices
291 // from mother to child.
292 FullMatrix<double> cell_interpolation(this->n_dofs_per_cell(),
293 this->n_dofs_per_cell());
294 FullMatrix<double> source_interpolation(this->n_dofs_per_cell(),
295 source_fe->n_dofs_per_cell());
296 FullMatrix<double> tmp(this->n_dofs_per_cell(),
297 source_fe->n_dofs_per_cell());
298 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
299 {
300 // generate a point on this
301 // cell and evaluate the
302 // shape functions there
303 const Point<dim> p = this->unit_support_points[j];
304 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
305 cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
306
307 for (unsigned int i = 0; i < source_fe->n_dofs_per_cell(); ++i)
308 source_interpolation(j, i) =
309 source_fe->poly_space->compute_value(i, p);
310 }
312 // then compute the
313 // interpolation matrix matrix
314 // for this coordinate
315 // direction
316 cell_interpolation.gauss_jordan();
317 cell_interpolation.mmult(interpolation_matrix, source_interpolation);
318
319 // cut off very small values
320 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
321 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
322 if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
323 interpolation_matrix(i, j) = 0.;
324
325 if constexpr (running_in_debug_mode())
327 // make sure that the row sum of
328 // each of the matrices is 1 at
329 // this point. this must be so
330 // since the shape functions sum up
331 // to 1
332 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
333 {
334 double sum = 0.;
335 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
336 sum += interpolation_matrix(i, j);
338 Assert(std::fabs(sum - 1) <
339 5e-14 * std::max(this->degree, 1U) * dim,
341 }
342 }
343 }
344 else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe))
345 {
346 // the element we want to interpolate from is an FE_Nothing. this
347 // element represents a function that is constant zero and has no
348 // degrees of freedom, so the interpolation is simply a multiplication
349 // with a n_dofs x 0 matrix. there is nothing to do here
350
351 // we would like to verify that the number of rows and columns of
352 // the matrix equals this->n_dofs_per_cell() and zero. unfortunately,
353 // whenever we do FullMatrix::reinit(m,0), it sets both rows and
354 // columns to zero, instead of m and zero. thus, only test the
355 // number of columns
356 Assert(interpolation_matrix.n() == x_source_fe.n_dofs_per_cell(),
357 ExcDimensionMismatch(interpolation_matrix.m(),
358 x_source_fe.n_dofs_per_cell()));
359 }
360 else
362 false,
363 (typename FiniteElement<dim,
364 spacedim>::ExcInterpolationNotImplemented()));
365}
367
368
369template <int dim, int spacedim>
370void
372 const FiniteElement<dim, spacedim> &x_source_fe,
373 FullMatrix<double> &interpolation_matrix,
374 const unsigned int) const
375{
376 // this is only implemented, if the source
377 // FE is also a DGQ element. in that case,
378 // both elements have no dofs on their
379 // faces and the face interpolation matrix
380 // is necessarily empty -- i.e. there isn't
381 // much we need to do here.
382 (void)interpolation_matrix;
384 AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
385 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. 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 typename FE::ExcInterpolationNotImplemented());
415
416 Assert(interpolation_matrix.m() == 0,
417 ExcDimensionMismatch(interpolation_matrix.m(), 0));
418 Assert(interpolation_matrix.n() == 0,
419 ExcDimensionMismatch(interpolation_matrix.m(), 0));
420}
421
422
423
424template <int dim, int spacedim>
425const FullMatrix<double> &
427 const unsigned int child,
428 const RefinementCase<dim> &refinement_case) const
429{
430 AssertIndexRange(refinement_case,
434 "Prolongation matrices are only available for refined cells!"));
436 child, this->reference_cell().template n_children<dim>(refinement_case));
437
438 // initialization upon first request
439 if (this->prolongation[refinement_case - 1][child].n() == 0)
440 {
441 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
442
443 // if matrix got updated while waiting for the lock
444 if (this->prolongation[refinement_case - 1][child].n() ==
445 this->n_dofs_per_cell())
446 return this->prolongation[refinement_case - 1][child];
447
448 // now do the work. need to get a non-const version of data in order to
449 // be able to modify them inside a const function
450 FE_DGQ<dim, spacedim> &this_nonconst =
451 const_cast<FE_DGQ<dim, spacedim> &>(*this);
452 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
453 {
454 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
456 isotropic_matrices.back().resize(
457 this->reference_cell().template n_children<dim>(
458 RefinementCase<dim>(refinement_case)),
459 FullMatrix<double>(this->n_dofs_per_cell(),
460 this->n_dofs_per_cell()));
461 if (dim == spacedim)
463 isotropic_matrices,
464 true);
465 else
467 isotropic_matrices,
468 true);
469 this_nonconst.prolongation[refinement_case - 1] =
470 std::move(isotropic_matrices.back());
471 }
472 else
473 {
474 // must compute both restriction and prolongation matrices because
475 // we only check for their size and the reinit call initializes them
476 // all
478 if (dim == spacedim)
479 {
481 this_nonconst.prolongation);
483 this_nonconst.restriction);
484 }
485 else
486 {
487 FE_DGQ<dim> tmp(this->degree);
489 this_nonconst.prolongation);
491 this_nonconst.restriction);
492 }
493 }
494 }
495
496 // finally return the matrix
497 return this->prolongation[refinement_case - 1][child];
498}
499
500
501
502template <int dim, int spacedim>
503const FullMatrix<double> &
505 const unsigned int child,
506 const RefinementCase<dim> &refinement_case) const
507{
508 AssertIndexRange(refinement_case,
512 "Restriction matrices are only available for refined cells!"));
514 child, this->reference_cell().template n_children<dim>(refinement_case));
515
516 // initialization upon first request
517 if (this->restriction[refinement_case - 1][child].n() == 0)
518 {
519 std::lock_guard<std::mutex> 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().template n_children<dim>(
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:378
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:182
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:504
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:371
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:132
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
Definition fe_dgq.cc:195
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:148
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:169
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:426
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:272
Definition fe_q.h:554
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_cell() const
std::vector< std::vector< FullMatrix< double > > > restriction
Definition fe.h:2523
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:2561
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2537
size_type n() const
size_type m() const
Definition point.h:113
const std::vector< Point< dim > > & get_points() const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
constexpr bool running_in_debug_mode()
Definition config.h:78
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#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:745
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:551
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:466
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)