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