Reference documentation for deal.II version 9.5.1
\(\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// 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
17#include <deal.II/base/quadrature.h>
18#include <deal.II/base/quadrature_lib.h>
19
20#include <deal.II/fe/fe.h>
21#include <deal.II/fe/fe_bernstein.h>
22#include <deal.II/fe/fe_dgq.h>
23#include <deal.II/fe/fe_nothing.h>
24#include <deal.II/fe/fe_pyramid_p.h>
25#include <deal.II/fe/fe_q.h>
26#include <deal.II/fe/fe_q_bubbles.h>
27#include <deal.II/fe/fe_q_dg0.h>
28#include <deal.II/fe/fe_q_hierarchical.h>
29#include <deal.II/fe/fe_q_iso_q1.h>
30#include <deal.II/fe/fe_simplex_p.h>
31#include <deal.II/fe/fe_simplex_p_bubbles.h>
32#include <deal.II/fe/fe_tools.h>
33#include <deal.II/fe/fe_wedge_p.h>
34
35#include <deal.II/lac/vector.h>
36
37#include <iostream>
38#include <memory>
39#include <sstream>
40
41
42DEAL_II_NAMESPACE_OPEN
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>(
68 TensorProductPolynomials<dim>(
69 Polynomials::generate_complete_Lagrange_basis(
70 internal::FE_DGQ::get_QGaussLobatto_points(degree))),
71 FiniteElementData<dim>(get_dpo_vector(degree),
72 1,
73 degree,
74 FiniteElementData<dim>::L2),
75 std::vector<bool>(
76 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
77 .n_dofs_per_cell(),
78 true),
79 std::vector<ComponentMask>(
80 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
81 .n_dofs_per_cell(),
82 std::vector<bool>(1, true)))
83{
84 // Compute support points, which are the tensor product of the Lagrange
85 // interpolation points in the constructor.
86 this->unit_support_points =
87 Quadrature<dim>(internal::FE_DGQ::get_QGaussLobatto_points(degree))
88 .get_points();
89
90 // do not initialize embedding and restriction here. these matrices are
91 // initialized on demand in get_restriction_matrix and
92 // get_prolongation_matrix
93
94 // note: no face support points for DG elements
95}
96
97
98
99template <int dim, int spacedim>
100FE_DGQ<dim, spacedim>::FE_DGQ(
101 const std::vector<Polynomials::Polynomial<double>> &polynomials)
102 : FE_Poly<dim, spacedim>(
103 TensorProductPolynomials<dim>(polynomials),
104 FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
105 1,
106 polynomials.size() - 1,
107 FiniteElementData<dim>::L2),
108 std::vector<bool>(FiniteElementData<dim>(get_dpo_vector(
109 polynomials.size() - 1),
110 1,
111 polynomials.size() - 1)
112 .n_dofs_per_cell(),
113 true),
114 std::vector<ComponentMask>(
115 FiniteElementData<dim>(get_dpo_vector(polynomials.size() - 1),
116 1,
117 polynomials.size() - 1)
118 .n_dofs_per_cell(),
119 std::vector<bool>(1, true)))
120{
121 // No support points can be defined in general. Derived classes might define
122 // support points like the class FE_DGQArbitraryNodes
123
124 // do not initialize embedding and restriction here. these matrices are
125 // initialized on demand in get_restriction_matrix and
126 // get_prolongation_matrix
127}
128
129
130
131template <int dim, int spacedim>
132std::string
133FE_DGQ<dim, spacedim>::get_name() const
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
149FE_DGQ<dim, spacedim>::convert_generalized_support_point_values_to_dof_values(
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>>
170FE_DGQ<dim, spacedim>::clone() const
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>
183FE_DGQ<dim, spacedim>::get_dpo_vector(const unsigned int deg)
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
196FE_DGQ<dim, spacedim>::rotate_indices(std::vector<unsigned int> &numbers,
197 const char direction) const
198{
199 const unsigned int n = this->degree + 1;
200 unsigned int s = n;
201 for (unsigned int i = 1; i < dim; ++i)
202 s *= n;
203 numbers.resize(s);
204
205 unsigned int l = 0;
206
207 if (dim == 1)
208 {
209 // Mirror around midpoint
210 for (unsigned int i = n; i > 0;)
211 numbers[l++] = --i;
212 }
213 else
214 {
215 switch (direction)
216 {
217 // Rotate xy-plane
218 // counter-clockwise
219 case 'z':
220 for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
221 for (unsigned int j = 0; j < n; ++j)
222 for (unsigned int i = 0; i < n; ++i)
223 {
224 unsigned int k = n * i - j + n - 1 + n * n * iz;
225 numbers[l++] = k;
226 }
227 break;
228 // Rotate xy-plane
229 // clockwise
230 case 'Z':
231 for (unsigned int iz = 0; iz < ((dim > 2) ? n : 1); ++iz)
232 for (unsigned int iy = 0; iy < n; ++iy)
233 for (unsigned int ix = 0; ix < n; ++ix)
234 {
235 unsigned int k = n * ix - iy + n - 1 + n * n * iz;
236 numbers[k] = l++;
237 }
238 break;
239 // Rotate yz-plane
240 // counter-clockwise
241 case 'x':
242 Assert(dim > 2, ExcDimensionMismatch(dim, 3));
243 for (unsigned int iz = 0; iz < n; ++iz)
244 for (unsigned int iy = 0; iy < n; ++iy)
245 for (unsigned int ix = 0; ix < n; ++ix)
246 {
247 unsigned int k = n * (n * iy - iz + n - 1) + ix;
248 numbers[l++] = k;
249 }
250 break;
251 // Rotate yz-plane
252 // clockwise
253 case 'X':
254 Assert(dim > 2, ExcDimensionMismatch(dim, 3));
255 for (unsigned int iz = 0; iz < n; ++iz)
256 for (unsigned int iy = 0; iy < n; ++iy)
257 for (unsigned int ix = 0; ix < n; ++ix)
258 {
259 unsigned int k = n * (n * iy - iz + n - 1) + ix;
260 numbers[k] = l++;
261 }
262 break;
263 default:
264 Assert(false, ExcNotImplemented());
265 }
266 }
267}
268
269
270
271template <int dim, int spacedim>
272void
273FE_DGQ<dim, spacedim>::get_interpolation_matrix(
274 const FiniteElement<dim, spacedim> &x_source_fe,
275 FullMatrix<double> & interpolation_matrix) const
276{
277 // this is only implemented, if the
278 // source FE is also a
279 // DGQ element
280 using FE = FiniteElement<dim, spacedim>;
281 AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
282 nullptr),
283 typename FE::ExcInterpolationNotImplemented());
284
285 // ok, source is a Q element, so
286 // we will be able to do the work
287 const FE_DGQ<dim, spacedim> &source_fe =
288 dynamic_cast<const FE_DGQ<dim, spacedim> &>(x_source_fe);
289
290 Assert(interpolation_matrix.m() == this->n_dofs_per_cell(),
291 ExcDimensionMismatch(interpolation_matrix.m(),
292 this->n_dofs_per_cell()));
293 Assert(interpolation_matrix.n() == source_fe.n_dofs_per_cell(),
294 ExcDimensionMismatch(interpolation_matrix.n(),
295 source_fe.n_dofs_per_cell()));
296
297
298 // compute the interpolation
299 // matrices in much the same way as
300 // we do for the embedding matrices
301 // from mother to child.
302 FullMatrix<double> cell_interpolation(this->n_dofs_per_cell(),
303 this->n_dofs_per_cell());
304 FullMatrix<double> source_interpolation(this->n_dofs_per_cell(),
305 source_fe.n_dofs_per_cell());
306 FullMatrix<double> tmp(this->n_dofs_per_cell(), source_fe.n_dofs_per_cell());
307 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
308 {
309 // generate a point on this
310 // cell and evaluate the
311 // shape functions there
312 const Point<dim> p = this->unit_support_points[j];
313 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
314 cell_interpolation(j, i) = this->poly_space->compute_value(i, p);
315
316 for (unsigned int i = 0; i < source_fe.n_dofs_per_cell(); ++i)
317 source_interpolation(j, i) = source_fe.poly_space->compute_value(i, p);
318 }
319
320 // then compute the
321 // interpolation matrix matrix
322 // for this coordinate
323 // direction
324 cell_interpolation.gauss_jordan();
325 cell_interpolation.mmult(interpolation_matrix, source_interpolation);
326
327 // cut off very small values
328 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
329 for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
330 if (std::fabs(interpolation_matrix(i, j)) < 1e-15)
331 interpolation_matrix(i, j) = 0.;
332
333#ifdef DEBUG
334 // make sure that the row sum of
335 // each of the matrices is 1 at
336 // this point. this must be so
337 // since the shape functions sum up
338 // to 1
339 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
340 {
341 double sum = 0.;
342 for (unsigned int j = 0; j < source_fe.n_dofs_per_cell(); ++j)
343 sum += interpolation_matrix(i, j);
344
345 Assert(std::fabs(sum - 1) < 5e-14 * std::max(this->degree, 1U) * dim,
346 ExcInternalError());
347 }
348#endif
349}
350
351
352
353template <int dim, int spacedim>
354void
355FE_DGQ<dim, spacedim>::get_face_interpolation_matrix(
356 const FiniteElement<dim, spacedim> &x_source_fe,
357 FullMatrix<double> & interpolation_matrix,
358 const unsigned int) const
359{
360 // this is only implemented, if the source
361 // FE is also a DGQ element. in that case,
362 // both elements have no dofs on their
363 // faces and the face interpolation matrix
364 // is necessarily empty -- i.e. there isn't
365 // much we need to do here.
366 (void)interpolation_matrix;
367 using FE = FiniteElement<dim, spacedim>;
368 AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
369 nullptr),
370 typename FE::ExcInterpolationNotImplemented());
371
372 Assert(interpolation_matrix.m() == 0,
373 ExcDimensionMismatch(interpolation_matrix.m(), 0));
374 Assert(interpolation_matrix.n() == 0,
375 ExcDimensionMismatch(interpolation_matrix.m(), 0));
376}
377
378
379
380template <int dim, int spacedim>
381void
382FE_DGQ<dim, spacedim>::get_subface_interpolation_matrix(
383 const FiniteElement<dim, spacedim> &x_source_fe,
384 const unsigned int,
385 FullMatrix<double> &interpolation_matrix,
386 const unsigned int) const
387{
388 // this is only implemented, if the source
389 // FE is also a DGQ element. in that case,
390 // both elements have no dofs on their
391 // faces and the face interpolation matrix
392 // is necessarily empty -- i.e. there isn't
393 // much we need to do here.
394 (void)interpolation_matrix;
395 using FE = FiniteElement<dim, spacedim>;
396 AssertThrow((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&x_source_fe) !=
397 nullptr),
398 typename FE::ExcInterpolationNotImplemented());
399
400 Assert(interpolation_matrix.m() == 0,
401 ExcDimensionMismatch(interpolation_matrix.m(), 0));
402 Assert(interpolation_matrix.n() == 0,
403 ExcDimensionMismatch(interpolation_matrix.m(), 0));
404}
405
406
407
408template <int dim, int spacedim>
409const FullMatrix<double> &
410FE_DGQ<dim, spacedim>::get_prolongation_matrix(
411 const unsigned int child,
412 const RefinementCase<dim> &refinement_case) const
413{
414 AssertIndexRange(refinement_case,
415 RefinementCase<dim>::isotropic_refinement + 1);
416 Assert(refinement_case != RefinementCase<dim>::no_refinement,
417 ExcMessage(
418 "Prolongation matrices are only available for refined cells!"));
419 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
420
421 // initialization upon first request
422 if (this->prolongation[refinement_case - 1][child].n() == 0)
423 {
424 std::lock_guard<std::mutex> lock(this->mutex);
425
426 // if matrix got updated while waiting for the lock
427 if (this->prolongation[refinement_case - 1][child].n() ==
428 this->n_dofs_per_cell())
429 return this->prolongation[refinement_case - 1][child];
430
431 // now do the work. need to get a non-const version of data in order to
432 // be able to modify them inside a const function
433 FE_DGQ<dim, spacedim> &this_nonconst =
434 const_cast<FE_DGQ<dim, spacedim> &>(*this);
435 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
436 {
437 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
438 RefinementCase<dim>::isotropic_refinement);
439 isotropic_matrices.back().resize(
440 GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
441 FullMatrix<double>(this->n_dofs_per_cell(),
442 this->n_dofs_per_cell()));
443 if (dim == spacedim)
444 FETools::compute_embedding_matrices(*this,
445 isotropic_matrices,
446 true);
447 else
448 FETools::compute_embedding_matrices(FE_DGQ<dim>(this->degree),
449 isotropic_matrices,
450 true);
451 this_nonconst.prolongation[refinement_case - 1].swap(
452 isotropic_matrices.back());
453 }
454 else
455 {
456 // must compute both restriction and prolongation matrices because
457 // we only check for their size and the reinit call initializes them
458 // all
459 this_nonconst.reinit_restriction_and_prolongation_matrices();
460 if (dim == spacedim)
461 {
462 FETools::compute_embedding_matrices(*this,
463 this_nonconst.prolongation);
464 FETools::compute_projection_matrices(*this,
465 this_nonconst.restriction);
466 }
467 else
468 {
469 FE_DGQ<dim> tmp(this->degree);
470 FETools::compute_embedding_matrices(tmp,
471 this_nonconst.prolongation);
472 FETools::compute_projection_matrices(tmp,
473 this_nonconst.restriction);
474 }
475 }
476 }
477
478 // finally return the matrix
479 return this->prolongation[refinement_case - 1][child];
480}
481
482
483
484template <int dim, int spacedim>
485const FullMatrix<double> &
486FE_DGQ<dim, spacedim>::get_restriction_matrix(
487 const unsigned int child,
488 const RefinementCase<dim> &refinement_case) const
489{
490 AssertIndexRange(refinement_case,
491 RefinementCase<dim>::isotropic_refinement + 1);
492 Assert(refinement_case != RefinementCase<dim>::no_refinement,
493 ExcMessage(
494 "Restriction matrices are only available for refined cells!"));
495 AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
496
497 // initialization upon first request
498 if (this->restriction[refinement_case - 1][child].n() == 0)
499 {
500 std::lock_guard<std::mutex> lock(this->mutex);
501
502 // if matrix got updated while waiting for the lock...
503 if (this->restriction[refinement_case - 1][child].n() ==
504 this->n_dofs_per_cell())
505 return this->restriction[refinement_case - 1][child];
506
507 // now do the work. need to get a non-const version of data in order to
508 // be able to modify them inside a const function
509 FE_DGQ<dim, spacedim> &this_nonconst =
510 const_cast<FE_DGQ<dim, spacedim> &>(*this);
511 if (refinement_case == RefinementCase<dim>::isotropic_refinement)
512 {
513 std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
514 RefinementCase<dim>::isotropic_refinement);
515 isotropic_matrices.back().resize(
516 GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
517 FullMatrix<double>(this->n_dofs_per_cell(),
518 this->n_dofs_per_cell()));
519 if (dim == spacedim)
520 FETools::compute_projection_matrices(*this,
521 isotropic_matrices,
522 true);
523 else
524 FETools::compute_projection_matrices(FE_DGQ<dim>(this->degree),
525 isotropic_matrices,
526 true);
527 this_nonconst.restriction[refinement_case - 1].swap(
528 isotropic_matrices.back());
529 }
530 else
531 {
532 // must compute both restriction and prolongation matrices because
533 // we only check for their size and the reinit call initializes them
534 // all
535 this_nonconst.reinit_restriction_and_prolongation_matrices();
536 if (dim == spacedim)
537 {
538 FETools::compute_embedding_matrices(*this,
539 this_nonconst.prolongation);
540 FETools::compute_projection_matrices(*this,
541 this_nonconst.restriction);
542 }
543 else
544 {
545 FE_DGQ<dim> tmp(this->degree);
546 FETools::compute_embedding_matrices(tmp,
547 this_nonconst.prolongation);
548 FETools::compute_projection_matrices(tmp,
549 this_nonconst.restriction);
550 }
551 }
552 }
553
554 // finally return the matrix
555 return this->restriction[refinement_case - 1][child];
556}
557
558
559
560template <int dim, int spacedim>
561bool
566
567
568
569template <int dim, int spacedim>
570std::vector<std::pair<unsigned int, unsigned int>>
571FE_DGQ<dim, spacedim>::hp_vertex_dof_identities(
572 const FiniteElement<dim, spacedim> & /*fe_other*/) const
573{
574 // this element is discontinuous, so by definition there can
575 // be no identities between its dofs and those of any neighbor
576 // (of whichever type the neighbor may be -- after all, we have
577 // no face dofs on this side to begin with)
578 return std::vector<std::pair<unsigned int, unsigned int>>();
579}
580
581
582
583template <int dim, int spacedim>
584std::vector<std::pair<unsigned int, unsigned int>>
585FE_DGQ<dim, spacedim>::hp_line_dof_identities(
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>>
599FE_DGQ<dim, spacedim>::hp_quad_dof_identities(
600 const FiniteElement<dim, spacedim> & /*fe_other*/,
601 const unsigned int) const
602{
603 // this element is discontinuous, so by definition there can
604 // be no identities between its dofs and those of any neighbor
605 // (of whichever type the neighbor may be -- after all, we have
606 // no face dofs on this side to begin with)
607 return std::vector<std::pair<unsigned int, unsigned int>>();
608}
609
610
611
612template <int dim, int spacedim>
613FiniteElementDomination::Domination
614FE_DGQ<dim, spacedim>::compare_for_domination(
615 const FiniteElement<dim, spacedim> &fe_other,
616 const unsigned int codim) const
617{
618 Assert(codim <= dim, ExcImpossibleInDim(dim));
619
620 // vertex/line/face domination
621 // ---------------------------
622 if (codim > 0)
623 // this is a discontinuous element, so by definition there will
624 // be no constraints wherever this element comes together with
625 // any other kind of element
626 return FiniteElementDomination::no_requirements;
627
628 // cell domination
629 // ---------------
630 // The following block of conditionals is rather ugly, but there is currently
631 // no other way how to deal with a robust comparison of FE_DGQ elements with
632 // relevant others in the current implementation.
633 if (const FE_DGQ<dim, spacedim> *fe_dgq_other =
634 dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe_other))
635 {
636 if (this->degree < fe_dgq_other->degree)
637 return FiniteElementDomination::this_element_dominates;
638 else if (this->degree == fe_dgq_other->degree)
639 return FiniteElementDomination::either_element_can_dominate;
640 else
641 return FiniteElementDomination::other_element_dominates;
642 }
643 else if (const FE_Q<dim, spacedim> *fe_q_other =
644 dynamic_cast<const FE_Q<dim, spacedim> *>(&fe_other))
645 {
646 if (this->degree < fe_q_other->degree)
647 return FiniteElementDomination::this_element_dominates;
648 else if (this->degree == fe_q_other->degree)
649 return FiniteElementDomination::either_element_can_dominate;
650 else
651 return FiniteElementDomination::other_element_dominates;
652 }
653 else if (const FE_Bernstein<dim, spacedim> *fe_bernstein_other =
654 dynamic_cast<const FE_Bernstein<dim, spacedim> *>(&fe_other))
655 {
656 if (this->degree < fe_bernstein_other->degree)
657 return FiniteElementDomination::this_element_dominates;
658 else if (this->degree == fe_bernstein_other->degree)
659 return FiniteElementDomination::either_element_can_dominate;
660 else
661 return FiniteElementDomination::other_element_dominates;
662 }
663 else if (const FE_Q_Bubbles<dim, spacedim> *fe_bubbles_other =
664 dynamic_cast<const FE_Q_Bubbles<dim, spacedim> *>(&fe_other))
665 {
666 if (this->degree < fe_bubbles_other->degree)
667 return FiniteElementDomination::this_element_dominates;
668 else if (this->degree == fe_bubbles_other->degree)
669 return FiniteElementDomination::either_element_can_dominate;
670 else
671 return FiniteElementDomination::other_element_dominates;
672 }
673 else if (const FE_Q_DG0<dim, spacedim> *fe_dg0_other =
674 dynamic_cast<const FE_Q_DG0<dim, spacedim> *>(&fe_other))
675 {
676 if (this->degree < fe_dg0_other->degree)
677 return FiniteElementDomination::this_element_dominates;
678 else if (this->degree == fe_dg0_other->degree)
679 return FiniteElementDomination::either_element_can_dominate;
680 else
681 return FiniteElementDomination::other_element_dominates;
682 }
683 else if (const FE_Q_iso_Q1<dim, spacedim> *fe_q_iso_q1_other =
684 dynamic_cast<const FE_Q_iso_Q1<dim, spacedim> *>(&fe_other))
685 {
686 if (this->degree < fe_q_iso_q1_other->degree)
687 return FiniteElementDomination::this_element_dominates;
688 else if (this->degree == fe_q_iso_q1_other->degree)
689 return FiniteElementDomination::either_element_can_dominate;
690 else
691 return FiniteElementDomination::other_element_dominates;
692 }
693 else if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
694 dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
695 {
696 if (this->degree < fe_hierarchical_other->degree)
697 return FiniteElementDomination::this_element_dominates;
698 else if (this->degree == fe_hierarchical_other->degree)
699 return FiniteElementDomination::either_element_can_dominate;
700 else
701 return FiniteElementDomination::other_element_dominates;
702 }
703 else if (const FE_SimplexDGP<dim, spacedim> *fe_dgp_other =
704 dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe_other))
705 {
706 if (this->degree < fe_dgp_other->degree)
707 return FiniteElementDomination::this_element_dominates;
708 else if (this->degree == fe_dgp_other->degree)
709 return FiniteElementDomination::either_element_can_dominate;
710 else
711 return FiniteElementDomination::other_element_dominates;
712 }
713 else if (const FE_Nothing<dim> *fe_nothing =
714 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
715 {
716 if (fe_nothing->is_dominating())
717 return FiniteElementDomination::other_element_dominates;
718 else
719 // the FE_Nothing has no degrees of freedom and it is typically used
720 // in a context where we don't require any continuity along the
721 // interface
722 return FiniteElementDomination::no_requirements;
723 }
724
725 Assert(false, ExcNotImplemented());
726 return FiniteElementDomination::neither_element_dominates;
727}
728
729
730
731template <int dim, int spacedim>
732bool
733FE_DGQ<dim, spacedim>::has_support_on_face(const unsigned int shape_index,
734 const unsigned int face_index) const
735{
736 AssertIndexRange(shape_index, this->n_dofs_per_cell());
737 AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
738
739 unsigned int n = this->degree + 1;
740
741 // For DGQ elements that do not define support points, we cannot define
742 // whether they have support at the boundary easily, so return true in any
743 // case
744 if (this->unit_support_points.empty())
745 return true;
746
747 // for DGQ(0) elements or arbitrary node DGQ with support points not located
748 // at the element boundary, the single shape functions is constant and
749 // therefore lives on the boundary
750 bool support_points_on_boundary = true;
751 for (unsigned int d = 0; d < dim; ++d)
752 if (std::abs(this->unit_support_points[0][d]) > 1e-13)
753 support_points_on_boundary = false;
754 for (unsigned int d = 0; d < dim; ++d)
755 if (std::abs(this->unit_support_points.back()[d] - 1.) > 1e-13)
756 support_points_on_boundary = false;
757 if (support_points_on_boundary == false)
758 return true;
759
760 unsigned int n2 = n * n;
761
762 switch (dim)
763 {
764 case 1:
765 {
766 // in 1d, things are simple. since
767 // there is only one degree of
768 // freedom per vertex in this
769 // class, the first is on vertex 0
770 // (==face 0 in some sense), the
771 // second on face 1:
772 return (((shape_index == 0) && (face_index == 0)) ||
773 ((shape_index == this->degree) && (face_index == 1)));
774 }
775
776 case 2:
777 {
778 if (face_index == 0 && (shape_index % n) == 0)
779 return true;
780 if (face_index == 1 && (shape_index % n) == this->degree)
781 return true;
782 if (face_index == 2 && shape_index < n)
783 return true;
784 if (face_index == 3 && shape_index >= this->n_dofs_per_cell() - n)
785 return true;
786 return false;
787 }
788
789 case 3:
790 {
791 const unsigned int in2 = shape_index % n2;
792
793 // x=0
794 if (face_index == 0 && (shape_index % n) == 0)
795 return true;
796 // x=1
797 if (face_index == 1 && (shape_index % n) == n - 1)
798 return true;
799 // y=0
800 if (face_index == 2 && in2 < n)
801 return true;
802 // y=1
803 if (face_index == 3 && in2 >= n2 - n)
804 return true;
805 // z=0
806 if (face_index == 4 && shape_index < n2)
807 return true;
808 // z=1
809 if (face_index == 5 && shape_index >= this->n_dofs_per_cell() - n2)
810 return true;
811 return false;
812 }
813
814 default:
815 Assert(false, ExcNotImplemented());
816 }
817 return true;
818}
819
820
821
822template <int dim, int spacedim>
823std::pair<Table<2, bool>, std::vector<unsigned int>>
824FE_DGQ<dim, spacedim>::get_constant_modes() const
825{
826 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
827 constant_modes.fill(true);
828 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
829 constant_modes, std::vector<unsigned int>(1, 0));
830}
831
832
833
834// ------------------------------ FE_DGQArbitraryNodes -----------------------
835
836template <int dim, int spacedim>
837FE_DGQArbitraryNodes<dim, spacedim>::FE_DGQArbitraryNodes(
838 const Quadrature<1> &points)
839 : FE_DGQ<dim, spacedim>(
840 Polynomials::generate_complete_Lagrange_basis(points.get_points()))
841{
842 Assert(points.size() > 0,
843 (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
844 this->unit_support_points = Quadrature<dim>(points).get_points();
845}
846
847
848
849template <int dim, int spacedim>
850std::string
851FE_DGQArbitraryNodes<dim, spacedim>::get_name() const
852{
853 // note that the FETools::get_fe_by_name function does not work for
854 // FE_DGQArbitraryNodes since there is no initialization by a degree value.
855 std::ostringstream namebuf;
856 bool equidistant = true;
857 std::vector<double> points(this->degree + 1);
858
859 auto *const polynomial_space =
860 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
861 Assert(polynomial_space != nullptr, ExcInternalError());
862 std::vector<unsigned int> lexicographic =
863 polynomial_space->get_numbering_inverse();
864 for (unsigned int j = 0; j <= this->degree; ++j)
865 points[j] = this->unit_support_points[lexicographic[j]][0];
866
867 // Check whether the support points are equidistant.
868 for (unsigned int j = 0; j <= this->degree; ++j)
869 if (std::abs(points[j] - static_cast<double>(j) / this->degree) > 1e-15)
870 {
871 equidistant = false;
872 break;
873 }
874 if (this->degree == 0 && std::abs(points[0] - 0.5) < 1e-15)
875 equidistant = true;
876
877 if (equidistant == true)
878 {
879 if (this->degree > 2)
880 namebuf << "FE_DGQArbitraryNodes<"
881 << Utilities::dim_string(dim, spacedim)
882 << ">(QIterated(QTrapezoid()," << this->degree << "))";
883 else
884 namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
885 << this->degree << ")";
886 return namebuf.str();
887 }
888
889 // Check whether the support points come from QGaussLobatto.
890 const QGaussLobatto<1> points_gl(this->degree + 1);
891 bool gauss_lobatto = true;
892 for (unsigned int j = 0; j <= this->degree; ++j)
893 if (points[j] != points_gl.point(j)(0))
894 {
895 gauss_lobatto = false;
896 break;
897 }
898
899 if (gauss_lobatto == true)
900 {
901 namebuf << "FE_DGQ<" << Utilities::dim_string(dim, spacedim) << ">("
902 << this->degree << ")";
903 return namebuf.str();
904 }
905
906 // Check whether the support points come from QGauss.
907 const QGauss<1> points_g(this->degree + 1);
908 bool gauss = true;
909 for (unsigned int j = 0; j <= this->degree; ++j)
910 if (points[j] != points_g.point(j)(0))
911 {
912 gauss = false;
913 break;
914 }
915
916 if (gauss == true)
917 {
918 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
919 << ">(QGauss(" << this->degree + 1 << "))";
920 return namebuf.str();
921 }
922
923 // Check whether the support points come from QGauss.
924 const QGaussLog<1> points_glog(this->degree + 1);
925 bool gauss_log = true;
926 for (unsigned int j = 0; j <= this->degree; ++j)
927 if (points[j] != points_glog.point(j)(0))
928 {
929 gauss_log = false;
930 break;
931 }
932
933 if (gauss_log == true)
934 {
935 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
936 << ">(QGaussLog(" << this->degree + 1 << "))";
937 return namebuf.str();
938 }
939
940 // All guesses exhausted
941 namebuf << "FE_DGQArbitraryNodes<" << Utilities::dim_string(dim, spacedim)
942 << ">(QUnknownNodes(" << this->degree + 1 << "))";
943 return namebuf.str();
944}
945
946
947
948template <int dim, int spacedim>
949void
950FE_DGQArbitraryNodes<dim, spacedim>::
951 convert_generalized_support_point_values_to_dof_values(
952 const std::vector<Vector<double>> &support_point_values,
953 std::vector<double> & nodal_values) const
954{
955 AssertDimension(support_point_values.size(),
956 this->get_unit_support_points().size());
957 AssertDimension(support_point_values.size(), nodal_values.size());
958 AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
959
960 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
961 {
962 AssertDimension(support_point_values[i].size(), 1);
963
964 nodal_values[i] = support_point_values[i](0);
965 }
966}
967
968
969
970template <int dim, int spacedim>
971std::unique_ptr<FiniteElement<dim, spacedim>>
972FE_DGQArbitraryNodes<dim, spacedim>::clone() const
973{
974 // Construct a dummy quadrature formula containing the FE's nodes:
975 std::vector<Point<1>> qpoints(this->degree + 1);
976 auto *const polynomial_space =
977 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
978 Assert(polynomial_space != nullptr, ExcInternalError());
979 std::vector<unsigned int> lexicographic =
980 polynomial_space->get_numbering_inverse();
981 for (unsigned int i = 0; i <= this->degree; ++i)
982 qpoints[i] = Point<1>(this->unit_support_points[lexicographic[i]][0]);
983 Quadrature<1> pquadrature(qpoints);
984
985 return std::make_unique<FE_DGQArbitraryNodes<dim, spacedim>>(pquadrature);
986}
987
988
989
990// ---------------------------------- FE_DGQLegendre -------------------------
991
992template <int dim, int spacedim>
993FE_DGQLegendre<dim, spacedim>::FE_DGQLegendre(const unsigned int degree)
994 : FE_DGQ<dim, spacedim>(
995 Polynomials::Legendre::generate_complete_basis(degree))
996{}
997
998
999
1000template <int dim, int spacedim>
1001std::pair<Table<2, bool>, std::vector<unsigned int>>
1002FE_DGQLegendre<dim, spacedim>::get_constant_modes() const
1003{
1004 // Legendre represents a constant function by one in the first basis
1005 // function and zero in all others
1006 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
1007 constant_modes(0, 0) = true;
1008 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
1009 constant_modes, std::vector<unsigned int>(1, 0));
1010}
1011
1012
1013
1014template <int dim, int spacedim>
1015std::string
1016FE_DGQLegendre<dim, spacedim>::get_name() const
1017{
1018 return "FE_DGQLegendre<" + Utilities::dim_string(dim, spacedim) + ">(" +
1019 Utilities::int_to_string(this->degree) + ")";
1020}
1021
1022
1023
1024template <int dim, int spacedim>
1025std::unique_ptr<FiniteElement<dim, spacedim>>
1026FE_DGQLegendre<dim, spacedim>::clone() const
1027{
1028 return std::make_unique<FE_DGQLegendre<dim, spacedim>>(this->degree);
1029}
1030
1031
1032
1033// ---------------------------------- FE_DGQHermite --------------------------
1034
1035template <int dim, int spacedim>
1036FE_DGQHermite<dim, spacedim>::FE_DGQHermite(const unsigned int degree)
1037 : FE_DGQ<dim, spacedim>(
1038 Polynomials::HermiteLikeInterpolation::generate_complete_basis(degree))
1039{}
1040
1041
1042
1043template <int dim, int spacedim>
1044std::string
1045FE_DGQHermite<dim, spacedim>::get_name() const
1046{
1047 return "FE_DGQHermite<" + Utilities::dim_string(dim, spacedim) + ">(" +
1048 Utilities::int_to_string(this->degree) + ")";
1049}
1050
1051
1052
1053template <int dim, int spacedim>
1054std::unique_ptr<FiniteElement<dim, spacedim>>
1055FE_DGQHermite<dim, spacedim>::clone() const
1056{
1057 return std::make_unique<FE_DGQHermite<dim, spacedim>>(this->degree);
1058}
1059
1060
1061
1062// explicit instantiations
1063#include "fe_dgq.inst"
1064
1065
1066DEAL_II_NAMESPACE_CLOSE
ArrayView::size
std::size_t size() const
Definition array_view.h:576
FE_DGQArbitraryNodes::get_name
virtual std::string get_name() const override
Definition fe_dgq.cc:851
FE_DGQArbitraryNodes::FE_DGQArbitraryNodes
FE_DGQArbitraryNodes(const Quadrature< 1 > &points)
Definition fe_dgq.cc:837
FE_DGQArbitraryNodes::convert_generalized_support_point_values_to_dof_values
virtual void convert_generalized_support_point_values_to_dof_values(const std::vector< Vector< double > > &support_point_values, std::vector< double > &nodal_values) const override
Definition fe_dgq.cc:951
FE_DGQArbitraryNodes::clone
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:972
FE_DGQHermite::get_name
virtual std::string get_name() const override
Definition fe_dgq.cc:1045
FE_DGQHermite::FE_DGQHermite
FE_DGQHermite(const unsigned int degree)
Definition fe_dgq.cc:1036
FE_DGQHermite::clone
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:1055
FE_DGQLegendre::get_constant_modes
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_dgq.cc:1002
FE_DGQLegendre::FE_DGQLegendre
FE_DGQLegendre(const unsigned int degree)
Definition fe_dgq.cc:993
FE_DGQLegendre::get_name
virtual std::string get_name() const override
Definition fe_dgq.cc:1016
FE_DGQLegendre::clone
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:1026
FE_DGQ
Definition fe_dgq.h:113
FE_DGQ::FE_DGQ
friend class FE_DGQ
Definition fe_dgq.h:377
FE_DGQ::has_support_on_face
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
Definition fe_dgq.cc:733
FE_DGQ::hp_constraints_are_implemented
virtual bool hp_constraints_are_implemented() const override
Definition fe_dgq.cc:562
FE_DGQ::hp_quad_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim, spacedim > &fe_other, const unsigned int face_no=0) const override
Definition fe_dgq.cc:599
FE_DGQ::get_dpo_vector
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
Definition fe_dgq.cc:183
FE_DGQ::get_constant_modes
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
Definition fe_dgq.cc:824
FE_DGQ::get_restriction_matrix
virtual const FullMatrix< double > & get_restriction_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition fe_dgq.cc:486
FE_DGQ::get_face_interpolation_matrix
virtual void get_face_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_dgq.cc:355
FE_DGQ::compare_for_domination
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim, spacedim > &fe_other, const unsigned int codim=0) const override final
Definition fe_dgq.cc:614
FE_DGQ::hp_vertex_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition fe_dgq.cc:571
FE_DGQ::get_name
virtual std::string get_name() const override
Definition fe_dgq.cc:133
FE_DGQ::rotate_indices
void rotate_indices(std::vector< unsigned int > &indices, const char direction) const
Definition fe_dgq.cc:196
FE_DGQ::hp_line_dof_identities
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim, spacedim > &fe_other) const override
Definition fe_dgq.cc:585
FE_DGQ::convert_generalized_support_point_values_to_dof_values
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
FE_DGQ::clone
virtual std::unique_ptr< FiniteElement< dim, spacedim > > clone() const override
Definition fe_dgq.cc:170
FE_DGQ::get_prolongation_matrix
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
Definition fe_dgq.cc:410
FE_DGQ::get_subface_interpolation_matrix
virtual void get_subface_interpolation_matrix(const FiniteElement< dim, spacedim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
Definition fe_dgq.cc:382
FE_DGQ::get_interpolation_matrix
virtual void get_interpolation_matrix(const FiniteElement< dim, spacedim > &source, FullMatrix< double > &matrix) const override
Definition fe_dgq.cc:273
FE_Poly
Definition fe_poly.h:77
FiniteElementData::degree
const unsigned int degree
Definition fe_data.h:453
FiniteElement::unit_support_points
std::vector< Point< dim > > unit_support_points
Definition fe.h:2440
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
fe_bernstein.h
fe_dgq.h
fe_nothing.h
fe_pyramid_p.h
fe_q.h
fe_q_bubbles.h
fe_q_dg0.h
fe_q_hierarchical.h
fe_q_iso_q1.h
fe_simplex_p.h
fe_simplex_p_bubbles.h
fe_tools.h
fe_wedge_p.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
StandardExceptions::ExcImpossibleInDim
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exceptions.h:1787
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
FETools::compute_embedding_matrices
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)
FETools::compute_projection_matrices
void compute_projection_matrices(const FiniteElement< dim, spacedim > &fe, std::vector< std::vector< FullMatrix< number > > > &matrices, const bool isotropic_only=false)
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Utilities::dim_string
std::string dim_string(const int dim, const int spacedim)
Definition utilities.cc:556
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
std
STL namespace.
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6024
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
quadrature_lib.h
Generated by doxygen 1.9.8