deal.II version GIT relicensing-2140-gf963fb7f8e 2024-11-11 12:50: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#ifdef DEBUG
326 // make sure that the row sum of
327 // each of the matrices is 1 at
328 // this point. this must be so
329 // since the shape functions sum up
330 // to 1
331 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
332 {
333 double sum = 0.;
334 for (unsigned int j = 0; j < source_fe->n_dofs_per_cell(); ++j)
335 sum += interpolation_matrix(i, j);
336
337 Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
339 }
340#endif
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}
364
365
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. 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 typename FE::ExcInterpolationNotImplemented());
385
386 Assert(interpolation_matrix.m() == 0,
387 ExcDimensionMismatch(interpolation_matrix.m(), 0));
388 Assert(interpolation_matrix.n() == 0,
389 ExcDimensionMismatch(interpolation_matrix.m(), 0));
390}
391
392
393
394template <int dim, int spacedim>
395void
397 const FiniteElement<dim, spacedim> &x_source_fe,
398 const unsigned int,
399 FullMatrix<double> &interpolation_matrix,
400 const unsigned int) const
401{
402 // this is only implemented, if the source
403 // FE is also a DGQ element. in that case,
404 // both elements have no dofs on their
405 // faces and the face interpolation matrix
406 // is necessarily empty -- i.e. there isn't
407 // much we need to do here.
408 (void)interpolation_matrix;
410 AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
411 nullptr),
412 typename FE::ExcInterpolationNotImplemented());
413
414 Assert(interpolation_matrix.m() == 0,
415 ExcDimensionMismatch(interpolation_matrix.m(), 0));
416 Assert(interpolation_matrix.n() == 0,
417 ExcDimensionMismatch(interpolation_matrix.m(), 0));
418}
419
420
421
422template <int dim, int spacedim>
423const FullMatrix<double> &
425 const unsigned int child,
426 const RefinementCase<dim> &refinement_case) const
427{
428 AssertIndexRange(refinement_case,
432 "Prolongation matrices are only available for refined cells!"));
433 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
434
435 // initialization upon first request
436 if (this->prolongation[refinement_case - 1][child].n() == 0)
437 {
438 std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
439
440 // if matrix got updated while waiting for the lock
441 if (this->prolongation[refinement_case - 1][child].n() ==
442 this->n_dofs_per_cell())
443 return this->prolongation[refinement_case - 1][child];
444
445 // now do the work. need to get a non-const version of data in order to
446 // be able to modify them inside a const function
447 FE_DGQ<dim, spacedim> &this_nonconst =
448 const_cast<FE_DGQ<dim, spacedim> &>(*this);
449 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
450 {
451 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
453 isotropic_matrices.back().resize(
455 FullMatrix<double>(this->n_dofs_per_cell(),
456 this->n_dofs_per_cell()));
457 if (dim == spacedim)
459 isotropic_matrices,
460 true);
461 else
463 isotropic_matrices,
464 true);
465 this_nonconst.prolongation[refinement_case - 1] =
466 std::move(isotropic_matrices.back());
467 }
468 else
469 {
470 // must compute both restriction and prolongation matrices because
471 // we only check for their size and the reinit call initializes them
472 // all
474 if (dim == spacedim)
475 {
477 this_nonconst.prolongation);
479 this_nonconst.restriction);
480 }
481 else
482 {
483 FE_DGQ<dim> tmp(this->degree);
485 this_nonconst.prolongation);
487 this_nonconst.restriction);
488 }
489 }
490 }
491
492 // finally return the matrix
493 return this->prolongation[refinement_case - 1][child];
494}
495
496
497
498template <int dim, int spacedim>
499const FullMatrix<double> &
501 const unsigned int child,
502 const RefinementCase<dim> &refinement_case) const
503{
504 AssertIndexRange(refinement_case,
508 "Restriction matrices are only available for refined cells!"));
509 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
510
511 // initialization upon first request
512 if (this->restriction[refinement_case - 1][child].n() == 0)
513 {
514 std::lock_guard<std::mutex> lock(restriction_matrix_mutex);
515
516 // if matrix got updated while waiting for the lock...
517 if (this->restriction[refinement_case - 1][child].n() ==
518 this->n_dofs_per_cell())
519 return this->restriction[refinement_case - 1][child];
520
521 // now do the work. need to get a non-const version of data in order to
522 // be able to modify them inside a const function
523 FE_DGQ<dim, spacedim> &this_nonconst =
524 const_cast<FE_DGQ<dim, spacedim> &>(*this);
525 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
526 {
527 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
529 isotropic_matrices.back().resize(
531 FullMatrix<double>(this->n_dofs_per_cell(),
532 this->n_dofs_per_cell()));
533 if (dim == spacedim)
535 isotropic_matrices,
536 true);
537 else
539 isotropic_matrices,
540 true);
541 this_nonconst.restriction[refinement_case - 1] =
542 std::move(isotropic_matrices.back());
543 }
544 else
545 {
546 // must compute both restriction and prolongation matrices because
547 // we only check for their size and the reinit call initializes them
548 // all
550 if (dim == spacedim)
551 {
553 this_nonconst.prolongation);
555 this_nonconst.restriction);
556 }
557 else
558 {
559 FE_DGQ<dim> tmp(this->degree);
561 this_nonconst.prolongation);
563 this_nonconst.restriction);
564 }
565 }
566 }
567
568 // finally return the matrix
569 return this->restriction[refinement_case - 1][child];
570}
571
572
573
574template <int dim, int spacedim>
575bool
580
581
582
583template <int dim, int spacedim>
584std::vector<std::pair<unsigned int, unsigned int>>
586 const FiniteElement<dim, spacedim> & /*fe_other*/) const
587{
588 // this element is discontinuous, so by definition there can
589 // be no identities between its dofs and those of any neighbor
590 // (of whichever type the neighbor may be -- after all, we have
591 // no face dofs on this side to begin with)
592 return std::vector<std::pair<unsigned int, unsigned int>>();
593}
594
595
596
597template <int dim, int spacedim>
598std::vector<std::pair<unsigned int, unsigned int>>
600 const FiniteElement<dim, spacedim> & /*fe_other*/) const
601{
602 // this element is discontinuous, so by definition there can
603 // be no identities between its dofs and those of any neighbor
604 // (of whichever type the neighbor may be -- after all, we have
605 // no face dofs on this side to begin with)
606 return std::vector<std::pair<unsigned int, unsigned int>>();
607}
608
609
610
611template <int dim, int spacedim>
612std::vector<std::pair<unsigned int, unsigned int>>
614 const FiniteElement<dim, spacedim> & /*fe_other*/,
615 const unsigned int) const
616{
617 // this element is discontinuous, so by definition there can
618 // be no identities between its dofs and those of any neighbor
619 // (of whichever type the neighbor may be -- after all, we have
620 // no face dofs on this side to begin with)
621 return std::vector<std::pair<unsigned int, unsigned int>>();
622}
623
624
625
626template <int dim, int spacedim>
629 const FiniteElement<dim, spacedim> &fe_other,
630 const unsigned int codim) const
631{
632 Assert(codim <= dim, ExcImpossibleInDim(dim));
633
634 // vertex/line/face domination
635 // ---------------------------
636 if (codim > 0)
637 // this is a discontinuous element, so by definition there will
638 // be no constraints wherever this element comes together with
639 // any other kind of element
641
642 // cell domination
643 // ---------------
644 // The following block of conditionals is rather ugly, but there is currently
645 // no other way how to deal with a robust comparison of FE_DGQ elements with
646 // relevant others in the current implementation.
647 if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
648 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
649 {
650 if (this->degree < fe_dgq_other->degree)
652 else if (this->degree == fe_dgq_other->degree)
654 else
656 }
657 else if (const FE_Q<dim, spacedim> *fe_q_other =
658 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
659 {
660 if (this->degree < fe_q_other->degree)
662 else if (this->degree == fe_q_other->degree)
664 else
666 }
667 else if (const FE_Bernstein<dim, spacedim> *fe_bernstein_other =
668 dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
669 {
670 if (this->degree < fe_bernstein_other->degree)
672 else if (this->degree == fe_bernstein_other->degree)
674 else
676 }
677 else if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
678 dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
679 {
680 if (this->degree < fe_bubbles_other->degree)
682 else if (this->degree == fe_bubbles_other->degree)
684 else
686 }
687 else if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
688 dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
689 {
690 if (this->degree < fe_dg0_other->degree)
692 else if (this->degree == fe_dg0_other->degree)
694 else
696 }
697 else if (const FE_Q_iso_Q1<dim, spacedim> *fe_q_iso_q1_other =
698 dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe_other))
699 {
700 if (this->degree < fe_q_iso_q1_other->degree)
702 else if (this->degree == fe_q_iso_q1_other->degree)
704 else
706 }
707 else if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
708 dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
709 {
710 if (this->degree < fe_hierarchical_other->degree)
712 else if (this->degree == fe_hierarchical_other->degree)
714 else
716 }
717 else if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
718 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
719 {
720 if (this->degree < fe_dgp_other->degree)
722 else if (this->degree == fe_dgp_other->degree)
724 else
726 }
727 else if (const FE_Nothing<dim> *fe_nothing =
728 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
729 {
730 if (fe_nothing->is_dominating())
732 else
733 // the FE_Nothing has no degrees of freedom and it is typically used
734 // in a context where we don't require any continuity along the
735 // interface
737 }
738
741}
742
743
744
745template <int dim, int spacedim>
746bool
747FE_DGQ<dim, spacedim>::has_support_on_face(const unsigned int shape_index,
748 const unsigned int face_index) const
749{
750 AssertIndexRange(shape_index, this->n_dofs_per_cell());
752
753 unsigned int n = this->degree + 1;
754
755 // For DGQ elements that do not define support points, we cannot define
756 // whether they have support at the boundary easily, so return true in any
757 // case
758 if (this->unit_support_points.empty())
759 return true;
760
761 // for DGQ(0) elements or arbitrary node DGQ with support points not located
762 // at the element boundary, the single shape functions is constant and
763 // therefore lives on the boundary
764 bool support_points_on_boundary = true;
765 for (unsigned int d = 0; d < dim; ++d)
766 if (std::abs(this->unit_support_points[0][d]) > 1e-13)
767 support_points_on_boundary = false;
768 for (unsigned int d = 0; d < dim; ++d)
769 if (std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
770 support_points_on_boundary = false;
771 if (support_points_on_boundary == false)
772 return true;
773
774 unsigned int n2 = n * n;
775
776 switch (dim)
777 {
778 case 1:
779 {
780 // in 1d, things are simple. since
781 // there is only one degree of
782 // freedom per vertex in this
783 // class, the first is on vertex 0
784 // (==face 0 in some sense), the
785 // second on face 1:
786 return (((shape_index == 0) && (face_index == 0)) ||
787 ((shape_index == this->degree) && (face_index == 1)));
788 }
789
790 case 2:
791 {
792 if (face_index == 0 && (shape_index % n) == 0)
793 return true;
794 if (face_index == 1 && (shape_index % n) == this->degree)
795 return true;
796 if (face_index == 2 && shape_index < n)
797 return true;
798 if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
799 return true;
800 return false;
801 }
802
803 case 3:
804 {
805 const unsigned int in2 = shape_index % n2;
806
807 // x=0
808 if (face_index == 0 && (shape_index % n) == 0)
809 return true;
810 // x=1
811 if (face_index == 1 && (shape_index % n) == n - 1)
812 return true;
813 // y=0
814 if (face_index == 2 && in2 < n)
815 return true;
816 // y=1
817 if (face_index == 3 && in2 >= n2 - n)
818 return true;
819 // z=0
820 if (face_index == 4 && shape_index < n2)
821 return true;
822 // z=1
823 if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
824 return true;
825 return false;
826 }
827
828 default:
830 }
831 return true;
832}
833
834
835
836template <int dim, int spacedim>
837std::pair<Table<2, bool>, std::vector<unsigned int>>
839{
840 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
841 constant_modes.fill(true);
842 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
843 constant_modes, std::vector<unsigned int>(1, 0));
844}
845
846
847
848// ------------------------------ FE_DGQArbitraryNodes -----------------------
849
850template <int dim, int spacedim>
852 const Quadrature<1> &points)
853 : FE_DGQ<dim, spacedim>(
854 Polynomials::generate_complete_Lagrange_basis(points.get_points()))
855{
856 Assert(points.size() > 0,
859}
860
861
862
863template <int dim, int spacedim>
864std::string
866{
867 // note that the FETools::get_fe_by_name function does not work for
868 // FE_DGQArbitraryNodes since there is no initialization by a degree value.
869 std::ostringstream namebuf;
870 bool equidistant = true;
871 std::vector<double> points(this->degree + 1);
872
873 const auto *const polynomial_space =
874 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
875 Assert(polynomial_space != nullptr, ExcInternalError());
876 std::vector<unsigned int> lexicographic =
877 polynomial_space->get_numbering_inverse();
878 for (unsigned int j = 0; j <= this->degree; ++j)
879 points[j] = this->unit_support_points[lexicographic[j]][0];
880
881 // Check whether the support points are equidistant.
882 for (unsigned int j = 0; j <= this->degree; ++j)
883 if (std::abs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
884 {
885 equidistant = false;
886 break;
887 }
888 if (this->degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
889 equidistant = true;
890
891 if (equidistant == true)
892 {
893 if (this->degree > 2)
894 namebuf << "FE_DGQArbitraryNodes<"
895 << Utilities::dim_string(dim, spacedim)
896 << ">(QIterated(QTrapezoid()," << this->degree << "))";
897 else
898 namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
899 << this->degree << ")";
900 return namebuf.str();
901 }
902
903 // Check whether the support points come from QGaussLobatto.
904 const QGaussLobatto<1> points_gl(this->degree + 1);
905 bool gauss_lobatto = true;
906 for (unsigned int j = 0; j <= this->degree; ++j)
907 if (points[j] != points_gl.point(j)[0])
908 {
909 gauss_lobatto = false;
910 break;
911 }
912
913 if (gauss_lobatto == true)
914 {
915 namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
916 << this->degree << ")";
917 return namebuf.str();
918 }
919
920 // Check whether the support points come from QGauss.
921 const QGauss<1> points_g(this->degree + 1);
922 bool gauss = true;
923 for (unsigned int j = 0; j <= this->degree; ++j)
924 if (points[j] != points_g.point(j)[0])
925 {
926 gauss = false;
927 break;
928 }
929
930 if (gauss == true)
931 {
932 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
933 << ">(QGauss(" << this->degree + 1 << "))";
934 return namebuf.str();
935 }
936
937 // Check whether the support points come from QGauss.
938 const QGaussLog<1> points_glog(this->degree + 1);
939 bool gauss_log = true;
940 for (unsigned int j = 0; j <= this->degree; ++j)
941 if (points[j] != points_glog.point(j)[0])
942 {
943 gauss_log = false;
944 break;
945 }
946
947 if (gauss_log == true)
948 {
949 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
950 << ">(QGaussLog(" << this->degree + 1 << "))";
951 return namebuf.str();
952 }
953
954 // All guesses exhausted
955 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
956 << ">(QUnknownNodes(" << this->degree + 1 << "))";
957 return namebuf.str();
958}
959
960
961
962template <int dim, int spacedim>
963void
966 const std::vector<Vector<double>> &support_point_values,
967 std::vector<double> &nodal_values) const
968{
969 AssertDimension(support_point_values.size(),
970 this->get_unit_support_points().size());
971 AssertDimension(support_point_values.size(), nodal_values.size());
972 AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
973
974 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
975 {
976 AssertDimension(support_point_values[i].size(), 1);
977
978 nodal_values[i] = support_point_values[i](0);
979 }
980}
981
982
983
984template <int dim, int spacedim>
985std::unique_ptr<FiniteElement<dim, spacedim>>
987{
988 // Construct a dummy quadrature formula containing the FE's nodes:
989 std::vector<Point<1>> qpoints(this->degree + 1);
990 const auto *const polynomial_space =
991 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
992 Assert(polynomial_space != nullptr, ExcInternalError());
993 std::vector<unsigned int> lexicographic =
994 polynomial_space->get_numbering_inverse();
995 for (unsigned int i = 0; i <= this->degree; ++i)
996 qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
997 const Quadrature<1> pquadrature(qpoints);
998
999 return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
1000}
1001
1002
1003
1004// ---------------------------------- FE_DGQLegendre -------------------------
1005
1006template <int dim, int spacedim>
1008 : FE_DGQ<dim, spacedim>(
1009 Polynomials::Legendre::generate_complete_basis(degree))
1010{}
1011
1012
1013
1014template <int dim, int spacedim>
1015std::pair<Table<2, bool>, std::vector<unsigned int>>
1017{
1018 // Legendre represents a constant function by one in the first basis
1019 // function and zero in all others
1020 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1021 constant_modes(0, 0) = true;
1022 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1023 constant_modes, std::vector<unsigned int>(1, 0));
1024}
1025
1026
1027
1028template <int dim, int spacedim>
1029std::string
1031{
1032 return "FE_DGQLegendre<" + Utilities::dim_string(dim, spacedim) + ">(" +
1033 Utilities::int_to_string(this->degree) + ")";
1034}
1035
1036
1037
1038template <int dim, int spacedim>
1039std::unique_ptr<FiniteElement<dim, spacedim>>
1041{
1042 return std::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1043}
1044
1045
1046
1047// ---------------------------------- FE_DGQHermite --------------------------
1048
1049template <int dim, int spacedim>
1051 : FE_DGQ<dim, spacedim>(
1052 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1053{}
1054
1055
1056
1057template <int dim, int spacedim>
1058std::string
1060{
1061 return "FE_DGQHermite<" + Utilities::dim_string(dim, spacedim) + ">(" +
1062 Utilities::int_to_string(this->degree) + ")";
1063}
1064
1065
1066
1067template <int dim, int spacedim>
1068std::unique_ptr<FiniteElement<dim, spacedim>>
1070{
1071 return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1072}
1073
1074
1075
1076// explicit instantiations
1077#include "fe_dgq.inst"
1078
1079
virtual std::string get_name() const override
Definition fe_dgq.cc:865
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
Definition fe_dgq.cc:851
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:965
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:986
virtual std::string get_name() const override
Definition fe_dgq.cc:1059
FE_DGQHermite(const unsigned int degree)
Definition fe_dgq.cc:1050
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:1069
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_dgq.cc:1016
FE_DGQLegendre(const unsigned int degree)
Definition fe_dgq.cc:1007
virtual std::string get_name() const override
Definition fe_dgq.cc:1030
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:1040
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:747
virtual bool hp_constraints_are_implemented() const override
Definition fe_dgq.cc:576
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:613
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:838
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:500
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:628
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:585
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:599
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:424
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:396
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:2431
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:2469
std::vector< std::vector< FullMatrix< double > > > prolongation
Definition fe.h:2445
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:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
#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)
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
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:555
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)