deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_q_hierarchical.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2002 - 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
16#include <deal.II/fe/fe_dgq.h>
19
20#include <cmath>
21#include <memory>
22#include <sstream>
23
24// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
25// adjust_line_dof_index_for_line_orientation_table fields, and write tests
26// similar to bits/face_orientation_and_fe_q_*
27
28
30
31namespace internal
32{
34 {
35 namespace
36 {
40 inline std::vector<unsigned int>
41 invert_numbering(const std::vector<unsigned int> &in)
42 {
43 std::vector<unsigned int> out(in.size());
44 for (unsigned int i = 0; i < in.size(); ++i)
45 {
46 AssertIndexRange(in[i], out.size());
47 out[in[i]] = i;
48 }
49 return out;
50 }
51 } // namespace
52 } // namespace FE_Q_Hierarchical
53} // namespace internal
54
55
56
57template <int dim>
60 Polynomials::Hierarchical::generate_complete_basis(degree)),
61 FiniteElementData<dim>(get_dpo_vector(degree),
62 1,
63 degree,
64 FiniteElementData<dim>::H1),
65 std::vector<bool>(
66 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
67 .n_dofs_per_cell(),
68 false),
69 std::vector<ComponentMask>(
70 FiniteElementData<dim>(get_dpo_vector(degree), 1, degree)
71 .n_dofs_per_cell(),
72 ComponentMask(std::vector<bool>(1, true))))
73 , face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering(degree))
74{
76 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
77 poly_space_derived_ptr->set_numbering(
79
80 // The matrix @p{dofs_cell} contains the
81 // values of the linear functionals of
82 // the parent 1d cell applied to the
83 // shape functions of the two 1d subcells.
84 // The matrix @p{dofs_subcell} contains
85 // the values of the linear functionals
86 // on each 1d subcell applied to the
87 // shape functions on the parent 1d
88 // subcell.
89 // We use @p{dofs_cell} and
90 // @p{dofs_subcell} to compute the
91 // @p{prolongation}, @p{restriction} and
92 // @p{interface_constraints} matrices
93 // for all dimensions.
94 std::vector<FullMatrix<double>> dofs_cell(
97 2 * this->n_dofs_per_vertex() +
98 this->n_dofs_per_line()));
99 std::vector<FullMatrix<double>> dofs_subcell(
102 2 * this->n_dofs_per_vertex() +
103 this->n_dofs_per_line()));
104 // build these fields, as they are
105 // needed as auxiliary fields later
106 // on
107 build_dofs_cell(dofs_cell, dofs_subcell);
108
109 // then use them to initialize
110 // other fields
113
114 // finally fill in support points
115 // on cell and face
118}
119
120
121
122template <int dim>
123std::string
125{
126 // note that the
127 // FETools::get_fe_by_name
128 // function depends on the
129 // particular format of the string
130 // this function returns, so they
131 // have to be kept in synch
132
133 std::ostringstream namebuf;
134 namebuf << "FE_Q_Hierarchical<" << dim << ">(" << this->degree << ")";
135
136 return namebuf.str();
137}
138
139
140
141template <int dim>
142std::unique_ptr<FiniteElement<dim, dim>>
144{
145 return std::make_unique<FE_Q_Hierarchical<dim>>(*this);
146}
147
148
149
150template <int dim>
151void
154 FullMatrix<double> &matrix) const
155{
156 // support interpolation between FE_Q_Hierarchical only.
158 dynamic_cast<const FE_Q_Hierarchical<dim> *>(&source))
159 {
160 // ok, source is a Q_Hierarchical element, so we will be able to do the
161 // work
162 Assert(matrix.m() == this->n_dofs_per_cell(),
163 ExcDimensionMismatch(matrix.m(), this->n_dofs_per_cell()));
164 Assert(matrix.n() == source.n_dofs_per_cell(),
165 ExcDimensionMismatch(matrix.m(), source_fe->n_dofs_per_cell()));
166
167 // Recall that DoFs are renumbered in the following order:
168 // vertices, lines, quads, hexes.
169 // As we deal with hierarchical FE, interpolation matrix is rather easy:
170 // it has 1 on pairs of dofs which are the same.
171 // To get those use get_embedding_dofs();
172
173 matrix = 0.;
174
175 // distinguish between the case when we interpolate to a richer element
176 if (this->n_dofs_per_cell() >= source_fe->n_dofs_per_cell())
177 {
178 const std::vector<unsigned int> dof_map =
179 this->get_embedding_dofs(source_fe->degree);
180 for (unsigned int j = 0; j < dof_map.size(); ++j)
181 matrix[dof_map[j]][j] = 1.;
182 }
183 // and when just truncate higher modes.
184 else
185 {
186 const std::vector<unsigned int> dof_map =
187 source_fe->get_embedding_dofs(this->degree);
188 for (unsigned int j = 0; j < dof_map.size(); ++j)
189 matrix[j][dof_map[j]] = 1.;
190 }
191 }
192 else
193 {
195 false,
197 "Interpolation is supported only between FE_Q_Hierarchical"));
198 }
199}
200
201template <int dim>
202const FullMatrix<double> &
204 const unsigned int child,
205 const RefinementCase<dim> &refinement_case) const
206{
207 Assert(
210 "Prolongation matrices are only available for isotropic refinement!"));
211
213 child, this->reference_cell().template n_children<dim>(refinement_case));
214
215 return this->prolongation[refinement_case - 1][child];
216}
217
218
219template <int dim>
220bool
225
226
227template <int dim>
228std::vector<std::pair<unsigned int, unsigned int>>
230 const FiniteElement<dim> &fe_other) const
231{
232 // we can presently only compute
233 // these identities if both FEs are
234 // FE_Q_Hierarchicals or if the other
235 // one is an FE_Nothing. in the first
236 // case, there should be exactly one
237 // single DoF of each FE at a
238 // vertex, and they should have
239 // identical value
240 if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
241 {
242 return std::vector<std::pair<unsigned int, unsigned int>>(
243 1, std::make_pair(0U, 0U));
244 }
245 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
246 {
247 // the FE_Nothing has no
248 // degrees of freedom, so there
249 // are no equivalencies to be
250 // recorded
251 return std::vector<std::pair<unsigned int, unsigned int>>();
252 }
253 else
254 {
256 return std::vector<std::pair<unsigned int, unsigned int>>();
257 }
258}
259
260template <int dim>
261std::vector<std::pair<unsigned int, unsigned int>>
263 const FiniteElement<dim> &fe_other) const
264{
265 // we can presently only compute
266 // these identities if both FEs are
267 // FE_Q_Hierarchicals or if the other
268 // one is an FE_Nothing.
269 if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
270 {
271 const unsigned int this_dpl = this->n_dofs_per_line();
272 const unsigned int other_dpl = fe_other.n_dofs_per_line();
273
274 // we deal with hierarchical 1d polynomials where dofs are enumerated
275 // increasingly. Thus we return a vector of pairs for the first N-1, where
276 // N is minimum number of dofs_per_line for each FE_Q_Hierarchical.
277 std::vector<std::pair<unsigned int, unsigned int>> res;
278 for (unsigned int i = 0; i < std::min(this_dpl, other_dpl); ++i)
279 res.emplace_back(i, i);
280
281 return res;
282 }
283 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
284 {
285 // the FE_Nothing has no
286 // degrees of freedom, so there
287 // are no equivalencies to be
288 // recorded
289 return std::vector<std::pair<unsigned int, unsigned int>>();
290 }
291 else
292 {
294 return std::vector<std::pair<unsigned int, unsigned int>>();
295 }
296}
297
298template <int dim>
299std::vector<std::pair<unsigned int, unsigned int>>
302 const unsigned int) const
303{
304 // we can presently only compute
305 // these identities if both FEs are
306 // FE_Q_Hierarchicals or if the other
307 // one is an FE_Nothing.
308 if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
309 {
310 // TODO: the implementation makes the assumption that all faces have the
311 // same number of dofs
312 AssertDimension(this->n_unique_faces(), 1);
313 const unsigned int face_no = 0;
314
315 const unsigned int this_dpq = this->n_dofs_per_quad(face_no);
316 const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no);
317
318 // we deal with hierarchical 1d polynomials where dofs are enumerated
319 // increasingly. Thus we return a vector of pairs for the first N-1, where
320 // N is minimum number of dofs_per_line for each FE_Q_Hierarchical.
321 std::vector<std::pair<unsigned int, unsigned int>> res;
322 for (unsigned int i = 0; i < std::min(this_dpq, other_dpq); ++i)
323 res.emplace_back(i, i);
324
325 return res;
326 }
327 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
328 {
329 // the FE_Nothing has no
330 // degrees of freedom, so there
331 // are no equivalencies to be
332 // recorded
333 return std::vector<std::pair<unsigned int, unsigned int>>();
334 }
335 else
336 {
338 return std::vector<std::pair<unsigned int, unsigned int>>();
339 }
340}
341
342template <int dim>
346 const unsigned int codim) const
347{
348 Assert(codim <= dim, ExcImpossibleInDim(dim));
349
350 // vertex/line/face domination
351 // (if fe_other is derived from FE_DGQ)
352 // ------------------------------------
353 if (codim > 0)
354 if (dynamic_cast<const FE_DGQ<dim> *>(&fe_other) != nullptr)
355 // there are no requirements between continuous and discontinuous elements
357
358 // vertex/line/face domination
359 // (if fe_other is not derived from FE_DGQ)
360 // & cell domination
361 // ----------------------------------------
363 dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
364 {
365 if (this->degree < fe_hierarchical_other->degree)
367 else if (this->degree == fe_hierarchical_other->degree)
369 else
371 }
372 else if (const FE_Nothing<dim> *fe_nothing =
373 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
374 {
375 if (fe_nothing->is_dominating())
377 else
378 // the FE_Nothing has no degrees of freedom and it is typically used
379 // in a context where we don't require any continuity along the
380 // interface
382 }
383
386}
387
388
389//---------------------------------------------------------------------------
390// Auxiliary functions
391//---------------------------------------------------------------------------
392
393
394template <int dim>
395void
397 std::vector<FullMatrix<double>> &dofs_cell,
398 std::vector<FullMatrix<double>> &dofs_subcell) const
399{
400 const unsigned int dofs_1d =
401 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
402
403 // The dofs_subcell matrices are transposed
404 // (4.19), (4.21) and (4.27),(4.28),(4.30) in
405 // Demkowicz, Oden, Rachowicz, Hardy, CMAMAE 77, 79-112, 1989
406 // so that
407 // DoFs_c(j) = dofs_subcell[c](j,k) dofs_cell(k)
408
409 // TODO: The dofs_subcell shall differ by a factor 2^p due to shape functions
410 // defined on [0,1] instead of [-1,1]. However, that does not seem to be
411 // the case. Perhaps this factor is added later on in auxiliary functions
412 // which use these matrices.
413
414 // dofs_cells[0](j,k):
415 // 0 1 | 2 3 4.
416 // 0 1 0 | .
417 // 1 0 0 | .
418 // -------------------
419 // 2 \ .
420 // 3 \ 2^k * k! / (k-j)!j!
421 // 4 \ .
422
423 // dofs_subcells[0](j,k):
424 // 0 1 | 2 3 4 5 6 .
425 // 0 1 0 | .
426 // 1 1/2 1/2 | -1 0 -1 0 -1.
427 // -------------------------------
428 // 2 \ .
429 // 3 \ (-1)^(k+j)/ 2^k * k!/(k-j)!j!
430 // 4 \ .
431
432 // dofs_cells[1](j,k):
433 // 0 1 | 2 3 4.
434 // 0 0 0 | .
435 // 1 0 1 | .
436 // -------------------
437 // 2 \ .
438 // 3 \ (-1)^(k+j) * 2^k * k!/(k-j)!j!
439 // 4 \.
440
441 // dofs_subcells[1](j,k):
442 // 0 1 | 2 3 4 5 6 .
443 // 0 1/2 1/2 | -1 0 -1 0 -1.
444 // 1 0 1 | .
445 // -------------------------------
446 // 2 \ .
447 // 3 \ 1/ 2^k * k!/(k-j)!j!
448 // 4 .
449
450 for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
451 for (unsigned int j = 0; j < dofs_1d; ++j)
452 for (unsigned int k = 0; k < dofs_1d; ++k)
453 {
454 // upper diagonal block
455 if ((j <= 1) && (k <= 1))
456 {
457 if (((c == 0) && (j == 0) && (k == 0)) ||
458 ((c == 1) && (j == 1) && (k == 1)))
459 dofs_cell[c](j, k) = 1.;
460 else
461 dofs_cell[c](j, k) = 0.;
462
463 if (((c == 0) && (j == 1)) || ((c == 1) && (j == 0)))
464 dofs_subcell[c](j, k) = .5;
465 else if (((c == 0) && (k == 0)) || ((c == 1) && (k == 1)))
466 dofs_subcell[c](j, k) = 1.;
467 else
468 dofs_subcell[c](j, k) = 0.;
469 }
470 // upper right block
471 else if ((j <= 1) && (k >= 2))
472 {
473 if (((c == 0) && (j == 1) && ((k % 2) == 0)) ||
474 ((c == 1) && (j == 0) && ((k % 2) == 0)))
475 dofs_subcell[c](j, k) = -1.;
476 }
477 // upper diagonal block
478 else if ((j >= 2) && (k >= 2) && (j <= k))
479 {
480 double factor = 1.;
481 for (unsigned int i = 1; i <= j; ++i)
482 factor *=
483 (static_cast<double>(k - i + 1)) / (static_cast<double>(i));
484 // factor == k * (k-1) * ... * (k-j+1) / j! = k! / (k-j)! / j!
485 if (c == 0)
486 {
487 dofs_subcell[c](j, k) = ((k + j) % 2 == 0) ?
488 Utilities::pow(.5, k) * factor :
489 -Utilities::pow(.5, k) * factor;
490 dofs_cell[c](j, k) = Utilities::pow(2, j) * factor;
491 }
492 else
493 {
494 dofs_subcell[c](j, k) = Utilities::pow(.5, k) * factor;
495 dofs_cell[c](j, k) = ((k + j) % 2 == 0) ?
496 Utilities::pow(2, j) * factor :
497 -Utilities::pow(2, j) * factor;
498 }
499 }
500 }
501}
502
503
504
505template <int dim>
506void
508 const std::vector<FullMatrix<double>> &dofs_subcell)
509{
510 const unsigned int dofs_1d =
511 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
512
513 this->interface_constraints.TableBase<2, double>::reinit(
514 this->interface_constraints_size());
515
516 switch (dim)
517 {
518 case 1:
519 {
520 // no constraints in 1d
521 break;
522 }
523
524 case 2:
525 {
526 // vertex node
527 for (unsigned int i = 0; i < dofs_1d; ++i)
528 this->interface_constraints(0, i) = dofs_subcell[0](1, i);
529 // edge nodes
530 for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
531 ++c)
532 for (unsigned int i = 0; i < dofs_1d; ++i)
533 for (unsigned int j = 2; j < dofs_1d; ++j)
534 this->interface_constraints(1 + c * (this->degree - 1) + j - 2,
535 i) = dofs_subcell[c](j, i);
536 break;
537 }
538
539 case 3:
540 {
541 for (unsigned int i = 0; i < dofs_1d * dofs_1d; ++i)
542 {
543 // center vertex node
544 this->interface_constraints(0, face_renumber[i]) =
545 dofs_subcell[0](1, i % dofs_1d) *
546 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
547
548 // boundary vertex nodes
549 this->interface_constraints(1, face_renumber[i]) =
550 dofs_subcell[0](0, i % dofs_1d) *
551 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
552 this->interface_constraints(2, face_renumber[i]) =
553 dofs_subcell[1](1, i % dofs_1d) *
554 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
555 this->interface_constraints(3, face_renumber[i]) =
556 dofs_subcell[0](1, i % dofs_1d) *
557 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
558 this->interface_constraints(4, face_renumber[i]) =
559 dofs_subcell[1](0, i % dofs_1d) *
560 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
561
562 // interior edges
563 for (unsigned int j = 0; j < (this->degree - 1); ++j)
564 {
565 this->interface_constraints(5 + j, face_renumber[i]) =
566 dofs_subcell[0](1, i % dofs_1d) *
567 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
568 this->interface_constraints(5 + (this->degree - 1) + j,
569 face_renumber[i]) =
570 dofs_subcell[0](1, i % dofs_1d) *
571 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
572 this->interface_constraints(5 + 2 * (this->degree - 1) + j,
573 face_renumber[i]) =
574 dofs_subcell[0](2 + j, i % dofs_1d) *
575 dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
576 this->interface_constraints(5 + 3 * (this->degree - 1) + j,
577 face_renumber[i]) =
578 dofs_subcell[1](2 + j, i % dofs_1d) *
579 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
580 }
581
582 // boundary edges
583 for (unsigned int j = 0; j < (this->degree - 1); ++j)
584 {
585 // left edge
586 this->interface_constraints(5 + 4 * (this->degree - 1) + j,
587 face_renumber[i]) =
588 dofs_subcell[0](0, i % dofs_1d) *
589 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
590 this->interface_constraints(5 + 4 * (this->degree - 1) +
591 (this->degree - 1) + j,
592 face_renumber[i]) =
593 dofs_subcell[0](0, i % dofs_1d) *
594 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
595 // right edge
596 this->interface_constraints(5 + 4 * (this->degree - 1) +
597 2 * (this->degree - 1) + j,
598 face_renumber[i]) =
599 dofs_subcell[1](1, i % dofs_1d) *
600 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
601 this->interface_constraints(5 + 4 * (this->degree - 1) +
602 3 * (this->degree - 1) + j,
603 face_renumber[i]) =
604 dofs_subcell[1](1, i % dofs_1d) *
605 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
606 // bottom edge
607 this->interface_constraints(5 + 4 * (this->degree - 1) +
608 4 * (this->degree - 1) + j,
609 face_renumber[i]) =
610 dofs_subcell[0](2 + j, i % dofs_1d) *
611 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
612 this->interface_constraints(5 + 4 * (this->degree - 1) +
613 5 * (this->degree - 1) + j,
614 face_renumber[i]) =
615 dofs_subcell[1](2 + j, i % dofs_1d) *
616 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
617 // top edge
618 this->interface_constraints(5 + 4 * (this->degree - 1) +
619 6 * (this->degree - 1) + j,
620 face_renumber[i]) =
621 dofs_subcell[0](2 + j, i % dofs_1d) *
622 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
623 this->interface_constraints(5 + 4 * (this->degree - 1) +
624 7 * (this->degree - 1) + j,
625 face_renumber[i]) =
626 dofs_subcell[1](2 + j, i % dofs_1d) *
627 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
628 }
629
630 // interior faces
631 for (unsigned int j = 0; j < (this->degree - 1); ++j)
632 for (unsigned int k = 0; k < (this->degree - 1); ++k)
633 {
634 // subcell 0
635 this->interface_constraints(5 + 12 * (this->degree - 1) +
636 j + k * (this->degree - 1),
637 face_renumber[i]) =
638 dofs_subcell[0](2 + j, i % dofs_1d) *
639 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
640 // subcell 1
641 this->interface_constraints(5 + 12 * (this->degree - 1) +
642 j + k * (this->degree - 1) +
643 (this->degree - 1) *
644 (this->degree - 1),
645 face_renumber[i]) =
646 dofs_subcell[1](2 + j, i % dofs_1d) *
647 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
648 // subcell 2
649 this->interface_constraints(5 + 12 * (this->degree - 1) +
650 j + k * (this->degree - 1) +
651 2 * (this->degree - 1) *
652 (this->degree - 1),
653 face_renumber[i]) =
654 dofs_subcell[0](2 + j, i % dofs_1d) *
655 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
656 // subcell 3
657 this->interface_constraints(5 + 12 * (this->degree - 1) +
658 j + k * (this->degree - 1) +
659 3 * (this->degree - 1) *
660 (this->degree - 1),
661 face_renumber[i]) =
662 dofs_subcell[1](2 + j, i % dofs_1d) *
663 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
664 }
665 }
666 break;
667 }
668
669 default:
671 }
672}
673
674
675
676template <int dim>
677void
679 const std::vector<FullMatrix<double>> &dofs_cell,
680 const std::vector<FullMatrix<double>> &dofs_subcell)
681{
683
684 const unsigned int dofs_1d =
685 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
687 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
688 const std::vector<unsigned int> &renumber =
689 poly_space_derived_ptr->get_numbering();
690
691 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
692 {
693 this->prolongation[iso][c].reinit(this->n_dofs_per_cell(),
694 this->n_dofs_per_cell());
695 this->restriction[iso][c].reinit(this->n_dofs_per_cell(),
696 this->n_dofs_per_cell());
697 }
698
699 // the 1d case is particularly
700 // simple, so special case it:
701 if (dim == 1)
702 {
703 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
704 ++c)
705 {
706 this->prolongation[iso][c].fill(dofs_subcell[c]);
707 this->restriction[iso][c].fill(dofs_cell[c]);
708 }
709 return;
710 }
711
712 // for higher dimensions, things
713 // are a little more tricky:
714
715 // j loops over dofs in the
716 // subcell. These are the rows in
717 // the embedding matrix.
718 //
719 // i loops over the dofs in the
720 // parent cell. These are the
721 // columns in the embedding matrix.
722 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
723 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
724 switch (dim)
725 {
726 case 2:
727 {
728 for (unsigned int c = 0;
730 ++c)
731 {
732 // left/right line: 0/1
733 const unsigned int c0 = c % 2;
734 // bottom/top line: 0/1
735 const unsigned int c1 = c / 2;
736
737 this->prolongation[iso][c](j, i) =
738 dofs_subcell[c0](renumber[j] % dofs_1d,
739 renumber[i] % dofs_1d) *
740 dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
741 dofs_1d,
742 (renumber[i] - (renumber[i] % dofs_1d)) /
743 dofs_1d);
744
745 this->restriction[iso][c](j, i) =
746 dofs_cell[c0](renumber[j] % dofs_1d,
747 renumber[i] % dofs_1d) *
748 dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
749 dofs_1d,
750 (renumber[i] - (renumber[i] % dofs_1d)) /
751 dofs_1d);
752 }
753 break;
754 }
755
756 case 3:
757 {
758 for (unsigned int c = 0;
760 ++c)
761 {
762 // left/right face: 0/1
763 const unsigned int c0 = c % 2;
764 // front/back face: 0/1
765 const unsigned int c1 = (c % 4) / 2;
766 // bottom/top face: 0/1
767 const unsigned int c2 = c / 4;
768
769 this->prolongation[iso][c](j, i) =
770 dofs_subcell[c0](renumber[j] % dofs_1d,
771 renumber[i] % dofs_1d) *
773 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
774 dofs_1d,
775 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
776 dofs_1d) *
778 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
779 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
780 dofs_1d)) /
781 dofs_1d,
782 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
783 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
784 dofs_1d)) /
785 dofs_1d);
786
787 this->restriction[iso][c](j, i) =
788 dofs_cell[c0](renumber[j] % dofs_1d,
789 renumber[i] % dofs_1d) *
790 dofs_cell[c1](
791 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
792 dofs_1d,
793 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
794 dofs_1d) *
795 dofs_cell[c2](
796 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
797 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
798 dofs_1d)) /
799 dofs_1d,
800 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
801 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
802 dofs_1d)) /
803 dofs_1d);
804 }
805 break;
806 }
807
808 default:
810 }
811}
812
813
814
815template <int dim>
816void
818{
819 // number of points: (degree+1)^dim
820 unsigned int n = this->degree + 1;
821 for (unsigned int i = 1; i < dim; ++i)
822 n *= this->degree + 1;
823
824 this->generalized_support_points.resize(n);
825
827 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
828 const std::vector<unsigned int> &index_map_inverse =
829 poly_space_derived_ptr->get_numbering_inverse();
830
831 Point<dim> p;
832 // the method of numbering allows
833 // each dof to be associated with a
834 // support point. There is
835 // only one support point per
836 // vertex, line, quad, hex, etc.
837 //
838 // note, however, that the support
839 // points thus associated with
840 // shape functions are not unique:
841 // the linear shape functions are
842 // associated with the vertices,
843 // but all others are associated
844 // with either line, quad, or hex
845 // midpoints, and there may be
846 // multiple shape functions
847 // associated with them. there
848 // really is no other useful
849 // numbering, since the
850 // hierarchical shape functions do
851 // not vanish at all-but-one
852 // interpolation points (like the
853 // Lagrange functions used in
854 // FE_Q), so there's not much we
855 // can do here.
856
857 // TODO shouldn't we just at least make support points unique,
858 // even though the delta property is not satisfied for this FE?
859 unsigned int k = 0;
860 for (unsigned int iz = 0; iz <= ((dim > 2) ? this->degree : 0); ++iz)
861 for (unsigned int iy = 0; iy <= ((dim > 1) ? this->degree : 0); ++iy)
862 for (unsigned int ix = 0; ix <= this->degree; ++ix)
863 {
864 if (ix == 0)
865 p[0] = 0.;
866 else if (ix == 1)
867 p[0] = 1.;
868 else
869 p[0] = .5;
870 if (dim > 1)
871 {
872 if (iy == 0)
873 p[1] = 0.;
874 else if (iy == 1)
875 p[1] = 1.;
876 else
877 p[1] = .5;
878 }
879 if (dim > 2)
880 {
881 if (iz == 0)
882 p[2] = 0.;
883 else if (iz == 1)
884 p[2] = 1.;
885 else
886 p[2] = .5;
887 }
888 this->generalized_support_points[index_map_inverse[k++]] = p;
889 }
890}
891
892
893
894template <>
895void
897{
898 // no faces in 1d, so nothing to do
899}
900
901
902template <>
903void
905 const FiniteElement<1, 1> & /*x_source_fe*/,
906 FullMatrix<double> & /*interpolation_matrix*/,
907 const unsigned int) const
908{
909 Assert(false, ExcImpossibleInDim(1));
910}
911
912
913template <>
914void
916 const FiniteElement<1, 1> & /*x_source_fe*/,
917 const unsigned int /*subface*/,
918 FullMatrix<double> & /*interpolation_matrix*/,
919 const unsigned int) const
920{
921 Assert(false, ExcImpossibleInDim(1));
922}
923
924
925
926template <int dim>
927void
931 const unsigned int face_no) const
932{
933 // this is only implemented, if the
934 // source FE is also a
935 // Q_Hierarchical element
937 AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
938 (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
939 nullptr),
941
942 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
944 this->n_dofs_per_face(face_no)));
945 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
947 x_source_fe.n_dofs_per_face(face_no)));
948
949 // ok, source is a Q_Hierarchical element, so
950 // we will be able to do the work
952 dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
954
955 // Make sure, that the element,
956 // for which the DoFs should be
957 // constrained is the one with
958 // the higher polynomial degree.
959 // Actually the procedure will work
960 // also if this assertion is not
961 // satisfied. But the matrices
962 // produced in that case might
963 // lead to problems in the
964 // hp-procedures, which use this
965 // method.
966 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
969
970 switch (dim)
971 {
972 case 2:
973 {
974 // In 2-dimension the constraints are trivial.
975 // First this->dofs_per_face DoFs of the constrained
976 // element are made equal to the current (dominating)
977 // element, which corresponds to 1 on diagonal of the matrix.
978 // DoFs which correspond to higher polynomials
979 // are zeroed (zero rows in the matrix).
980 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
981 interpolation_matrix(i, i) = 1;
982
983 break;
984 }
985
986 case 3:
987 {
988 for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i)
989 interpolation_matrix(i, i) = 1;
990
991 for (unsigned int i = 0; i < this->degree - 1; ++i)
992 {
993 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j)
994 interpolation_matrix(i + j * (x_source_fe.degree - 1) +
996 i + j * (this->degree - 1) +
998
999 for (unsigned int j = 0; j < this->degree - 1; ++j)
1001 (x_source_fe.degree - 1) +
1004 (this->degree - 1) +
1006 1;
1007 }
1008 }
1009 }
1010}
1011
1012
1013
1014template <int dim>
1015void
1018 const unsigned int subface,
1020 const unsigned int face_no) const
1021{
1022 // this is only implemented, if the
1023 // source FE is also a
1024 // Q_Hierarchical element
1026 AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
1027 (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
1028 nullptr),
1030
1031 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
1033 this->n_dofs_per_face(face_no)));
1034 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
1036 x_source_fe.n_dofs_per_face(face_no)));
1037
1038 // ok, source is a Q_Hierarchical element, so
1039 // we will be able to do the work
1041 dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
1042
1043 // Make sure, that the element,
1044 // for which the DoFs should be
1045 // constrained is the one with
1046 // the higher polynomial degree.
1047 // Actually the procedure will work
1048 // also if this assertion is not
1049 // satisfied. But the matrices
1050 // produced in that case might
1051 // lead to problems in the
1052 // hp-procedures, which use this
1053 // method.
1054 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
1056
1057 switch (dim)
1058 {
1059 case 2:
1060 {
1061 switch (subface)
1062 {
1063 case 0:
1064 {
1065 interpolation_matrix(0, 0) = 1.0;
1066 interpolation_matrix(1, 0) = 0.5;
1067 interpolation_matrix(1, 1) = 0.5;
1068
1069 for (unsigned int dof = 2;
1070 dof < this->n_dofs_per_face(face_no);)
1071 {
1072 interpolation_matrix(1, dof) = -1.0;
1073 dof = dof + 2;
1074 }
1075
1076 int factorial_i = 1;
1077
1078 for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1079 ++i)
1080 {
1081 interpolation_matrix(i, i) = Utilities::pow(0.5, i);
1082 factorial_i *= i;
1084 int factorial_ij = 1;
1085
1086 for (unsigned int j = i + 1;
1087 j < this->n_dofs_per_face(face_no);
1088 ++j)
1089 {
1090 factorial_ij *= j - i;
1091 factorial_j *= j;
1092
1093 if (((i + j) & 1) != 0u)
1095 -1.0 * Utilities::pow(0.5, j) * factorial_j /
1097
1098 else
1100 Utilities::pow(0.5, j) * factorial_j /
1102 }
1103 }
1104
1105 break;
1106 }
1107
1108 case 1:
1109 {
1110 interpolation_matrix(0, 0) = 0.5;
1111 interpolation_matrix(0, 1) = 0.5;
1112
1113 for (unsigned int dof = 2;
1114 dof < this->n_dofs_per_face(face_no);)
1115 {
1116 interpolation_matrix(0, dof) = -1.0;
1117 dof = dof + 2;
1118 }
1119
1120 interpolation_matrix(1, 1) = 1.0;
1121
1122 int factorial_i = 1;
1123
1124 for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1125 ++i)
1126 {
1127 interpolation_matrix(i, i) = Utilities::pow(0.5, i);
1128 factorial_i *= i;
1130 int factorial_ij = 1;
1131
1132 for (unsigned int j = i + 1;
1133 j < this->n_dofs_per_face(face_no);
1134 ++j)
1135 {
1136 factorial_ij *= j - i;
1137 factorial_j *= j;
1139 Utilities::pow(0.5, j) * factorial_j /
1141 }
1142 }
1143 }
1144 }
1145
1146 break;
1147 }
1148
1149 case 3:
1150 {
1151 switch (subface)
1152 {
1153 case 0:
1154 {
1155 interpolation_matrix(0, 0) = 1.0;
1156 interpolation_matrix(1, 0) = 0.5;
1157 interpolation_matrix(1, 1) = 0.5;
1158 interpolation_matrix(2, 0) = 0.5;
1159 interpolation_matrix(2, 2) = 0.5;
1160
1161 for (unsigned int i = 0; i < this->degree - 1;)
1162 {
1163 for (unsigned int line = 0;
1165 ++line)
1167 i + line * (this->degree - 1) +
1168 4) = -0.5;
1169
1170 for (unsigned int j = 0; j < this->degree - 1;)
1171 {
1173 i + (j + 4) * this->degree - j) =
1174 1.0;
1175 j = j + 2;
1176 }
1177
1178 interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1179 -1.0;
1180 interpolation_matrix(2, i + 4) = -1.0;
1181 i = i + 2;
1182 }
1183
1184 for (unsigned int vertex = 0;
1186 ++vertex)
1187 interpolation_matrix(3, vertex) = 0.25;
1188
1189 int factorial_i = 1;
1190
1191 for (unsigned int i = 2; i <= this->degree; ++i)
1192 {
1193 double tmp = Utilities::pow(0.5, i);
1194 interpolation_matrix(i + 2, i + 2) = tmp;
1195 interpolation_matrix(i + 2 * source_fe.degree,
1196 i + 2 * this->degree) = tmp;
1197 tmp *= 0.5;
1198 interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1199 tmp;
1200 interpolation_matrix(i + source_fe.degree + 1,
1201 i + this->degree + 1) = tmp;
1202 interpolation_matrix(i + 3 * source_fe.degree - 1,
1203 i + 2 * this->degree) = tmp;
1204 interpolation_matrix(i + 3 * source_fe.degree - 1,
1205 i + 3 * this->degree - 1) = tmp;
1206 tmp *= -2.0;
1207
1208 for (unsigned int j = 0; j < this->degree - 1;)
1209 {
1210 interpolation_matrix(i + source_fe.degree + 1,
1211 (i + 2) * this->degree + j + 2 -
1212 i) = tmp;
1213 interpolation_matrix(i + 3 * source_fe.degree - 1,
1214 i + (j + 4) * this->degree - j -
1215 2) = tmp;
1216 j = j + 2;
1217 }
1218
1219 int factorial_k = 1;
1220
1221 for (unsigned int j = 2; j <= this->degree; ++j)
1222 {
1223 interpolation_matrix(i + (j + 2) * source_fe.degree -
1224 j,
1225 i + (j + 2) * this->degree - j) =
1226 Utilities::pow(0.5, i + j);
1227 factorial_k *= j;
1228 int factorial_kl = 1;
1230
1231 for (unsigned int k = j + 1; k < this->degree; ++k)
1232 {
1233 factorial_kl *= k - j;
1234 factorial_l *= k;
1235
1236 if (((j + k) & 1) != 0u)
1238 i + (j + 2) * source_fe.degree - j,
1239 i + (k + 2) * this->degree - k) =
1240 -1.0 * Utilities::pow(0.5, i + k) *
1242
1243 else
1245 i + (j + 2) * source_fe.degree - j,
1246 i + (k + 2) * this->degree - k) =
1247 Utilities::pow(0.5, i + k) * factorial_l /
1249 }
1250 }
1251
1252 factorial_i *= i;
1254 int factorial_ij = 1;
1255
1256 for (unsigned int j = i + 1; j <= this->degree; ++j)
1257 {
1258 factorial_ij *= j - i;
1259 factorial_j *= j;
1260
1261 if (((i + j) & 1) != 0u)
1262 {
1263 tmp = -1.0 * Utilities::pow(0.5, j) *
1265 interpolation_matrix(i + 2, j + 2) = tmp;
1266 interpolation_matrix(i + 2 * source_fe.degree,
1267 j + 2 * this->degree) = tmp;
1268 factorial_k = 1;
1269
1270 for (unsigned int k = 2; k <= this->degree; ++k)
1271 {
1273 i + (k + 2) * source_fe.degree - k,
1274 j + (k + 2) * this->degree - k) =
1275 tmp * Utilities::pow(0.5, k);
1276 factorial_k *= k;
1278 int factorial_kl = 1;
1279
1280 for (unsigned int l = k + 1;
1281 l <= this->degree;
1282 ++l)
1283 {
1284 factorial_kl *= l - k;
1285 factorial_l *= l;
1286
1287 if (((k + l) & 1) != 0u)
1289 i + (k + 2) * source_fe.degree - k,
1290 j + (l + 2) * this->degree - l) =
1291 -1.0 * tmp * Utilities::pow(0.5, l) *
1292 factorial_l /
1294
1295 else
1297 i + (k + 2) * source_fe.degree - k,
1298 j + (l + 2) * this->degree - l) =
1299 tmp * Utilities::pow(0.5, l) *
1300 factorial_l /
1302 }
1303 }
1304
1305 tmp *= 0.5;
1306 interpolation_matrix(i + source_fe.degree + 1,
1307 j + 2) = tmp;
1308 interpolation_matrix(i + source_fe.degree + 1,
1309 j + this->degree + 1) = tmp;
1310 interpolation_matrix(i + 3 * source_fe.degree - 1,
1311 j + 2 * this->degree) = tmp;
1312 interpolation_matrix(i + 3 * source_fe.degree - 1,
1313 j + 3 * this->degree - 1) =
1314 tmp;
1315 tmp *= -2.0;
1316
1317 for (unsigned int k = 0; k < this->degree - 1;)
1318 {
1319 interpolation_matrix(i + source_fe.degree + 1,
1320 (j + 2) * this->degree +
1321 k + 2 - j) = tmp;
1323 i + 3 * source_fe.degree - 1,
1324 j + (k + 4) * this->degree - k - 2) = tmp;
1325 k = k + 2;
1326 }
1327 }
1328 else
1329 {
1330 tmp = Utilities::pow(0.5, j) * factorial_j /
1332 interpolation_matrix(i + 2, j + 2) = tmp;
1333 interpolation_matrix(i + 2 * source_fe.degree,
1334 j + 2 * this->degree) = tmp;
1335 factorial_k = 1;
1336
1337 for (unsigned int k = 2; k <= this->degree; ++k)
1338 {
1340 i + (k + 2) * source_fe.degree - k,
1341 j + (k + 2) * this->degree - k) =
1342 tmp * Utilities::pow(0.5, k);
1343 factorial_k *= k;
1345 int factorial_kl = 1;
1346
1347 for (unsigned int l = k + 1;
1348 l <= this->degree;
1349 ++l)
1350 {
1351 factorial_kl *= l - k;
1352 factorial_l *= l;
1353
1354 if (((k + l) & 1) != 0u)
1356 i + (k + 2) * source_fe.degree - k,
1357 j + (l + 2) * this->degree - l) =
1358 -1.0 * tmp * Utilities::pow(0.5, l) *
1359 factorial_l /
1361
1362 else
1364 i + (k + 2) * source_fe.degree - k,
1365 j + (l + 2) * this->degree - l) =
1366 tmp * Utilities::pow(0.5, l) *
1367 factorial_l /
1369 }
1370 }
1371
1372 tmp *= 0.5;
1373 interpolation_matrix(i + source_fe.degree + 1,
1374 j + 2) = tmp;
1375 interpolation_matrix(i + source_fe.degree + 1,
1376 j + this->degree + 1) = tmp;
1377 interpolation_matrix(i + 3 * source_fe.degree - 1,
1378 j + 2 * this->degree) = tmp;
1379 interpolation_matrix(i + 3 * source_fe.degree - 1,
1380 j + 3 * this->degree - 1) =
1381 tmp;
1382 tmp *= -2.0;
1383
1384 for (unsigned int k = 0; k < this->degree - 1;)
1385 {
1386 interpolation_matrix(i + source_fe.degree + 1,
1387 (j + 2) * this->degree +
1388 k + 2 - j) = tmp;
1390 i + 3 * source_fe.degree - 1,
1391 j + (k + 4) * this->degree - k - 2) = tmp;
1392 k = k + 2;
1393 }
1394 }
1395 }
1396 }
1397
1398 break;
1399 }
1400
1401 case 1:
1402 {
1403 interpolation_matrix(0, 0) = 0.5;
1404 interpolation_matrix(0, 1) = 0.5;
1405 interpolation_matrix(1, 1) = 1.0;
1406 interpolation_matrix(3, 1) = 0.5;
1407 interpolation_matrix(3, 3) = 0.5;
1408
1409 for (unsigned int i = 0; i < this->degree - 1;)
1410 {
1411 for (unsigned int line = 0;
1413 ++line)
1415 i + line * (this->degree - 1) +
1416 4) = -0.5;
1417
1418 for (unsigned int j = 0; j < this->degree - 1;)
1419 {
1421 i + (j + 4) * this->degree - j) =
1422 1.0;
1423 j = j + 2;
1424 }
1425
1426 interpolation_matrix(0, i + 2 * (this->degree + 1)) =
1427 -1.0;
1428 interpolation_matrix(3, i + this->degree + 3) = -1.0;
1429 i = i + 2;
1430 }
1431
1432 for (unsigned int vertex = 0;
1434 ++vertex)
1435 interpolation_matrix(2, vertex) = 0.25;
1436
1437 int factorial_i = 1;
1438
1439 for (unsigned int i = 2; i <= this->degree; ++i)
1440 {
1441 double tmp = Utilities::pow(0.5, i + 1);
1442 interpolation_matrix(i + 2, i + 2) = tmp;
1443 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1444 interpolation_matrix(i + 3 * source_fe.degree - 1,
1445 i + 2 * this->degree) = tmp;
1446 interpolation_matrix(i + 3 * source_fe.degree - 1,
1447 i + 3 * this->degree - 1) = tmp;
1448 tmp *= -2.0;
1449
1450 for (unsigned int j = 0; j < this->degree - 1;)
1451 {
1453 j + (i + 2) * this->degree + 2 -
1454 i) = tmp;
1455 interpolation_matrix(i + 3 * source_fe.degree - 1,
1456 i + (j + 4) * this->degree - j -
1457 2) = tmp;
1458 j = j + 2;
1459 }
1460
1461 tmp *= -1.0;
1462 interpolation_matrix(i + source_fe.degree + 1,
1463 i + this->degree + 1) = tmp;
1464 interpolation_matrix(i + 2 * source_fe.degree,
1465 i + 2 * this->degree) = tmp;
1466 factorial_i *= i;
1468 int factorial_ij = 1;
1469
1470 for (unsigned int j = i + 1; j <= this->degree; ++j)
1471 {
1472 factorial_ij *= j - i;
1473 factorial_j *= j;
1474 tmp = Utilities::pow(0.5, j) * factorial_j /
1476 interpolation_matrix(i + 2 * source_fe.degree,
1477 j + 2 * this->degree) = tmp;
1478 int factorial_k = 1;
1479
1480 for (unsigned int k = 2; k <= this->degree; ++k)
1481 {
1483 i + (k + 2) * source_fe.degree - k,
1484 j + (k + 2) * this->degree - k) =
1485 tmp * Utilities::pow(0.5, k);
1486 factorial_k *= k;
1488 int factorial_kl = 1;
1489
1490 for (unsigned int l = k + 1; l <= this->degree;
1491 ++l)
1492 {
1493 factorial_kl *= l - k;
1494 factorial_l *= l;
1495
1496 if (((k + l) & 1) != 0u)
1498 i + (k + 2) * source_fe.degree - k,
1499 j + (l + 2) * this->degree - l) =
1500 -1.0 * tmp * Utilities::pow(0.5, l) *
1501 factorial_l /
1503
1504 else
1506 i + (k + 2) * source_fe.degree - k,
1507 j + (l + 2) * this->degree - l) =
1508 tmp * Utilities::pow(0.5, l) *
1509 factorial_l /
1511 }
1512 }
1513
1514 tmp *= -1.0;
1515
1516 for (unsigned int k = 0; k < this->degree - 1;)
1517 {
1518 interpolation_matrix(i + 3 * source_fe.degree - 1,
1519 j + (k + 4) * this->degree -
1520 k - 2) = tmp;
1521 k = k + 2;
1522 }
1523
1524 tmp *= -0.5;
1525 interpolation_matrix(i + 3 * source_fe.degree - 1,
1526 j + 2 * this->degree) = tmp;
1527 interpolation_matrix(i + 3 * source_fe.degree - 1,
1528 j + 3 * this->degree - 1) = tmp;
1529
1530 if (((i + j) & 1) != 0u)
1531 tmp *= -1.0;
1532
1533 interpolation_matrix(i + 2, j + 2) = tmp;
1534 interpolation_matrix(i + 2, j + this->degree + 1) =
1535 tmp;
1536 interpolation_matrix(i + source_fe.degree + 1,
1537 j + this->degree + 1) =
1538 2.0 * tmp;
1539 tmp *= -2.0;
1540
1541 for (unsigned int k = 0; k < this->degree - 1;)
1542 {
1544 k + (j + 2) * this->degree +
1545 2 - j) = tmp;
1546 k = k + 2;
1547 }
1548 }
1549
1550 int factorial_k = 1;
1551
1552 for (unsigned int j = 2; j <= this->degree; ++j)
1553 {
1554 interpolation_matrix(i + (j + 2) * source_fe.degree -
1555 j,
1556 i + (j + 2) * this->degree - j) =
1557 Utilities::pow(0.5, i + j);
1558 factorial_k *= j;
1560 int factorial_kl = 1;
1561
1562 for (unsigned int k = j + 1; k <= this->degree; ++k)
1563 {
1564 factorial_kl *= k - j;
1565 factorial_l *= k;
1566
1567 if (((j + k) & 1) != 0u)
1569 i + (j + 2) * source_fe.degree - j,
1570 i + (k + 2) * this->degree - k) =
1571 -1.0 * Utilities::pow(0.5, i + k) *
1573
1574 else
1576 i + (j + 2) * source_fe.degree - j,
1577 i + (k + 2) * this->degree - k) =
1578 Utilities::pow(0.5, i + k) * factorial_l /
1580 }
1581 }
1582 }
1583
1584 break;
1585 }
1586
1587 case 2:
1588 {
1589 interpolation_matrix(0, 0) = 0.5;
1590 interpolation_matrix(0, 2) = 0.5;
1591 interpolation_matrix(2, 2) = 1.0;
1592 interpolation_matrix(3, 2) = 0.5;
1593 interpolation_matrix(3, 3) = 0.5;
1594
1595 for (unsigned int i = 0; i < this->degree - 1;)
1596 {
1597 for (unsigned int line = 0;
1599 ++line)
1601 i + line * (this->degree - 1) +
1602 4) = -0.5;
1603
1604 for (unsigned int j = 0; j < this->degree - 1;)
1605 {
1607 i + (j + 4) * this->degree - j) =
1608 1.0;
1609 j = j + 2;
1610 }
1611
1612 interpolation_matrix(0, i + 4) = -1.0;
1613 interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1614 i = i + 2;
1615 }
1616
1617 for (unsigned int vertex = 0;
1619 ++vertex)
1620 interpolation_matrix(1, vertex) = 0.25;
1621
1622 int factorial_i = 1;
1623
1624 for (unsigned int i = 2; i <= this->degree; ++i)
1625 {
1626 double tmp = Utilities::pow(0.5, i);
1627 interpolation_matrix(i + 2, i + 2) = tmp;
1628 interpolation_matrix(i + 3 * source_fe.degree - 1,
1629 i + 3 * this->degree - 1) = tmp;
1630 tmp *= 0.5;
1631 interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1632 tmp;
1633 interpolation_matrix(i + source_fe.degree + 1,
1634 i + this->degree + 1) = tmp;
1635 interpolation_matrix(i + 2 * source_fe.degree,
1636 i + 2 * this->degree) = tmp;
1637 interpolation_matrix(i + 2 * source_fe.degree,
1638 i + 3 * this->degree - 1) = tmp;
1639 tmp *= -2.0;
1640
1641 for (unsigned int j = 0; j < this->degree - 1;)
1642 {
1643 interpolation_matrix(i + source_fe.degree + 1,
1644 j + (i + 2) * this->degree + 2 -
1645 i) = tmp;
1646 interpolation_matrix(i + 2 * source_fe.degree,
1647 i + (j + 4) * this->degree - j -
1648 2) = tmp;
1649 j = j + 2;
1650 }
1651
1652 int factorial_k = 1;
1653
1654 for (unsigned int j = 2; j <= this->degree; ++j)
1655 {
1656 interpolation_matrix(i + (j + 2) * source_fe.degree -
1657 j,
1658 i + (j + 2) * this->degree - j) =
1659 Utilities::pow(0.5, i + j);
1660 factorial_k *= j;
1661 int factorial_kl = 1;
1663
1664 for (unsigned int k = j + 1; k <= this->degree; ++k)
1665 {
1666 factorial_kl *= k - j;
1667 factorial_l *= k;
1669 i + (j + 2) * source_fe.degree - j,
1670 i + (k + 2) * this->degree - k) =
1671 Utilities::pow(0.5, i + k) * factorial_l /
1673 }
1674 }
1675
1676 factorial_i *= i;
1678 int factorial_ij = 1;
1679
1680 for (unsigned int j = i + 1; j <= this->degree; ++j)
1681 {
1682 factorial_ij *= j - i;
1683 factorial_j *= j;
1684 tmp = Utilities::pow(0.5, j) * factorial_j /
1686 interpolation_matrix(i + 2, j + 2) = tmp;
1687 tmp *= -1.0;
1688
1689 for (unsigned int k = 0; k < this->degree - 1;)
1690 {
1691 interpolation_matrix(i + source_fe.degree + 1,
1692 k + (j + 2) * this->degree +
1693 2 - j) = tmp;
1694 k = k + 2;
1695 }
1696
1697 tmp *= -0.5;
1698 interpolation_matrix(i + source_fe.degree + 1,
1699 j + 2) = tmp;
1700 interpolation_matrix(i + source_fe.degree + 1,
1701 j + this->degree + 1) = tmp;
1702
1703 if (((i + j) & 1) != 0u)
1704 tmp *= -1.0;
1705
1706 interpolation_matrix(i + 2 * source_fe.degree,
1707 j + 2 * this->degree) = tmp;
1708 interpolation_matrix(i + 2 * source_fe.degree,
1709 j + 3 * this->degree - 1) = tmp;
1710 tmp *= 2.0;
1711 interpolation_matrix(i + 3 * source_fe.degree - 1,
1712 j + 3 * this->degree - 1) = tmp;
1713 int factorial_k = 1;
1714
1715 for (unsigned int k = 2; k <= this->degree; ++k)
1716 {
1718 i + (k + 2) * source_fe.degree - k,
1719 j + (k + 2) * this->degree - k) =
1720 tmp * Utilities::pow(0.5, k);
1721 factorial_k *= k;
1723 int factorial_kl = 1;
1724
1725 for (unsigned int l = k + 1; l <= this->degree;
1726 ++l)
1727 {
1728 factorial_kl *= l - k;
1729 factorial_l *= l;
1731 i + (k + 2) * source_fe.degree - k,
1732 j + (l + 2) * this->degree - l) =
1733 tmp * Utilities::pow(0.5, l) * factorial_l /
1735 }
1736 }
1737
1738 tmp *= -1.0;
1739
1740 for (unsigned int k = 0; k < this->degree - 1;)
1741 {
1742 interpolation_matrix(i + 2 * source_fe.degree,
1743 j + (k + 4) * this->degree -
1744 k - 2) = tmp;
1745 k = k + 2;
1746 }
1747 }
1748 }
1749
1750 break;
1751 }
1752
1753 case 3:
1754 {
1755 for (unsigned int vertex = 0;
1757 ++vertex)
1758 interpolation_matrix(0, vertex) = 0.25;
1759
1760 for (unsigned int i = 0; i < this->degree - 1;)
1761 {
1762 for (unsigned int line = 0;
1764 ++line)
1766 i + line * (this->degree - 1) +
1767 4) = -0.5;
1768
1769 for (unsigned int j = 0; j < this->degree - 1;)
1770 {
1772 i + (j + 4) * this->degree - j) =
1773 1.0;
1774 j = j + 2;
1775 }
1776
1777 interpolation_matrix(1, i + 4) = -1.0;
1778 interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
1779 i = i + 2;
1780 }
1781
1782 interpolation_matrix(1, 0) = 0.5;
1783 interpolation_matrix(1, 1) = 0.5;
1784 interpolation_matrix(2, 2) = 0.5;
1785 interpolation_matrix(2, 3) = 0.5;
1786 interpolation_matrix(3, 3) = 1.0;
1787
1788 int factorial_i = 1;
1789
1790 for (unsigned int i = 2; i <= this->degree; ++i)
1791 {
1792 double tmp = Utilities::pow(0.5, i + 1);
1793 interpolation_matrix(i + 2, i + 2) = tmp;
1794 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1795 interpolation_matrix(i + 2 * source_fe.degree,
1796 i + 2 * this->degree) = tmp;
1797 interpolation_matrix(i + 2 * source_fe.degree,
1798 i + 3 * this->degree - 1) = tmp;
1799 tmp *= -2.0;
1800
1801 for (unsigned int j = 0; j < this->degree - 1;)
1802 {
1804 j + (i + 2) * this->degree + 2 -
1805 i) = tmp;
1806 interpolation_matrix(i + 2 * source_fe.degree,
1807 i + (j + 4) * this->degree - 2) =
1808 tmp;
1809 j = j + 2;
1810 }
1811
1812 tmp *= -1.0;
1813 interpolation_matrix(i + source_fe.degree + 1,
1814 i + this->degree + 1) = tmp;
1815 interpolation_matrix(i + 3 * source_fe.degree - 1,
1816 i + 3 * this->degree - 1) = tmp;
1817 int factorial_k = 1;
1818
1819 for (unsigned int j = 2; j <= this->degree; ++j)
1820 {
1821 interpolation_matrix(i + (j + 2) * source_fe.degree -
1822 j,
1823 i + (j + 2) * this->degree - j) =
1824 Utilities::pow(0.5, i + j);
1825 factorial_k *= j;
1827 int factorial_kl = 1;
1828
1829 for (unsigned int k = j + 1; k <= this->degree; ++k)
1830 {
1831 factorial_kl *= k - j;
1832 factorial_l *= k;
1834 i + (j + 2) * source_fe.degree - j,
1835 i + (k + 2) * this->degree - k) =
1836 Utilities::pow(0.5, i + k) * factorial_l /
1838 }
1839 }
1840
1841 factorial_i *= i;
1843 int factorial_ij = 1;
1844
1845 for (unsigned int j = i + 1; j <= this->degree; ++j)
1846 {
1847 factorial_ij *= j - i;
1848 factorial_j *= j;
1849 tmp = Utilities::pow(0.5, j + 1) * factorial_j /
1851 interpolation_matrix(i + 2, j + 2) = tmp;
1852 interpolation_matrix(i + 2, j + this->degree + 1) =
1853 tmp;
1854 interpolation_matrix(i + 2 * source_fe.degree,
1855 j + 2 * this->degree) = tmp;
1856 interpolation_matrix(i + 2 * source_fe.degree,
1857 j + 3 * this->degree - 1) = tmp;
1858 tmp *= 2.0;
1859 interpolation_matrix(i + source_fe.degree + 1,
1860 j + this->degree + 1) = tmp;
1861 interpolation_matrix(i + 3 * source_fe.degree - 1,
1862 j + 3 * this->degree - 1) = tmp;
1863 int factorial_k = 1;
1864
1865 for (unsigned int k = 2; k <= this->degree; ++k)
1866 {
1868 i + (k + 2) * source_fe.degree - k,
1869 j + (k + 2) * this->degree - k) =
1870 tmp * Utilities::pow(0.5, k);
1871 factorial_k *= k;
1873 int factorial_kl = 1;
1874
1875 for (unsigned int l = k + 1; l <= this->degree;
1876 ++l)
1877 {
1878 factorial_kl *= l - k;
1879 factorial_l *= l;
1881 i + (k + 2) * source_fe.degree - k,
1882 j + (l + 2) * this->degree - l) =
1883 tmp * Utilities::pow(0.5, l) * factorial_l /
1885 }
1886 }
1887
1888 tmp *= -1.0;
1889
1890 for (unsigned int k = 0; k < this->degree - 1;)
1891 {
1893 k + (j + 2) * this->degree +
1894 2 - j) = tmp;
1895 interpolation_matrix(i + 2 * source_fe.degree,
1896 j + (k + 4) * this->degree -
1897 2) = tmp;
1898 k = k + 2;
1899 }
1900 }
1901 }
1902 }
1903 }
1904 }
1905 }
1906}
1907
1908
1909
1910template <int dim>
1911void
1913{
1914 const unsigned int codim = dim - 1;
1915
1916 // TODO: the implementation makes the assumption that all faces have the
1917 // same number of dofs
1918 AssertDimension(this->n_unique_faces(), 1);
1919 const unsigned int face_no = 0;
1920
1921 // number of points: (degree+1)^codim
1922 unsigned int n = this->degree + 1;
1923 for (unsigned int i = 1; i < codim; ++i)
1924 n *= this->degree + 1;
1925
1926 this->generalized_face_support_points[face_no].resize(n);
1927
1928 Point<codim> p;
1929
1930 unsigned int k = 0;
1931 for (unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
1932 for (unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
1933 for (unsigned int ix = 0; ix <= this->degree; ++ix)
1934 {
1935 if (ix == 0)
1936 p[0] = 0.;
1937 else if (ix == 1)
1938 p[0] = 1.;
1939 else
1940 p[0] = .5;
1941 if (codim > 1)
1942 {
1943 if (iy == 0)
1944 p[1] = 0.;
1945 else if (iy == 1)
1946 p[1] = 1.;
1947 else
1948 p[1] = .5;
1949 }
1950 if (codim > 2)
1951 {
1952 if (iz == 0)
1953 p[2] = 0.;
1954 else if (iz == 1)
1955 p[2] = 1.;
1956 else
1957 p[2] = .5;
1958 }
1959 this->generalized_face_support_points[face_no][face_renumber[k++]] =
1960 p;
1961 }
1962}
1963
1964
1965// we use same dpo_vector as FE_Q
1966template <int dim>
1967std::vector<unsigned int>
1969{
1970 std::vector<unsigned int> dpo(dim + 1, 1U);
1971 for (unsigned int i = 1; i < dpo.size(); ++i)
1972 dpo[i] = dpo[i - 1] * (deg - 1);
1973 return dpo;
1974}
1975
1976
1977
1978template <int dim>
1979std::vector<unsigned int>
1981 const FiniteElementData<dim> &fe)
1982{
1983 Assert(fe.n_components() == 1, ExcInternalError());
1984 std::vector<unsigned int> h2l(fe.n_dofs_per_cell());
1985
1986 // polynomial degree
1987 const unsigned int degree = fe.n_dofs_per_line() + 1;
1988 // number of grid points in each
1989 // direction
1990 const unsigned int n = degree + 1;
1991
1992 // the following lines of code are
1993 // somewhat odd, due to the way the
1994 // hierarchic numbering is
1995 // organized. if someone would
1996 // really want to understand these
1997 // lines, you better draw some
1998 // pictures where you indicate the
1999 // indices and orders of vertices,
2000 // lines, etc, along with the
2001 // numbers of the degrees of
2002 // freedom in hierarchical and
2003 // lexicographical order
2004 switch (dim)
2005 {
2006 case 1:
2007 {
2008 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2009 h2l[i] = i;
2010
2011 break;
2012 }
2013
2014 case 2:
2015 {
2016 // Example: degree=3
2017 //
2018 // hierarchical numbering:
2019 // 2 10 11 3
2020 // 5 14 15 7
2021 // 4 12 13 6
2022 // 0 8 9 1
2023 //
2024 // fe_q_hierarchical numbering:
2025 // 4 6 7 5
2026 // 12 14 15 13
2027 // 8 10 11 9
2028 // 0 2 3 1
2029 unsigned int next_index = 0;
2030 // first the four vertices
2031 h2l[next_index++] = 0;
2032 h2l[next_index++] = 1;
2033 h2l[next_index++] = n;
2034 h2l[next_index++] = n + 1;
2035 // left line
2036 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2037 h2l[next_index++] = (2 + i) * n;
2038 // right line
2039 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2040 h2l[next_index++] = (2 + i) * n + 1;
2041 // bottom line
2042 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2043 h2l[next_index++] = 2 + i;
2044 // top line
2045 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2046 h2l[next_index++] = n + 2 + i;
2047 // inside quad
2048 Assert(fe.n_dofs_per_quad(0 /*only one quad in 2d*/) ==
2049 fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2051 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2052 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2053 h2l[next_index++] = (2 + i) * n + 2 + j;
2054
2055 Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2056
2057 break;
2058 }
2059
2060 case 3:
2061 {
2062 unsigned int next_index = 0;
2063 const unsigned int n2 = n * n;
2064 // first the eight vertices
2065 // bottom face, lexicographic
2066 h2l[next_index++] = 0;
2067 h2l[next_index++] = 1;
2068 h2l[next_index++] = n;
2069 h2l[next_index++] = n + 1;
2070 // top face, lexicographic
2071 h2l[next_index++] = n2;
2072 h2l[next_index++] = n2 + 1;
2073 h2l[next_index++] = n2 + n;
2074 h2l[next_index++] = n2 + n + 1;
2075
2076 // now the lines
2077 // bottom face
2078 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2079 h2l[next_index++] = (2 + i) * n;
2080 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2081 h2l[next_index++] = (2 + i) * n + 1;
2082 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2083 h2l[next_index++] = 2 + i;
2084 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2085 h2l[next_index++] = n + 2 + i;
2086 // top face
2087 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2088 h2l[next_index++] = n2 + (2 + i) * n;
2089 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2090 h2l[next_index++] = n2 + (2 + i) * n + 1;
2091 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2092 h2l[next_index++] = n2 + 2 + i;
2093 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2094 h2l[next_index++] = n2 + n + 2 + i;
2095 // lines in z-direction
2096 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2097 h2l[next_index++] = (2 + i) * n2;
2098 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2099 h2l[next_index++] = (2 + i) * n2 + 1;
2100 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2101 h2l[next_index++] = (2 + i) * n2 + n;
2102 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2103 h2l[next_index++] = (2 + i) * n2 + n + 1;
2104
2105 // TODO: the implementation makes the assumption that all faces have
2106 // the same number of dofs
2107 AssertDimension(fe.n_unique_faces(), 1);
2108 const unsigned int face_no = 0;
2109 (void)face_no;
2110
2111 // inside quads
2112 Assert(fe.n_dofs_per_quad(face_no) ==
2113 fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2115 // left face
2116 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2117 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2118 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n;
2119 // right face
2120 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2121 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2122 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 1;
2123 // front face
2124 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2125 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2126 h2l[next_index++] = (2 + i) * n2 + 2 + j;
2127 // back face
2128 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2129 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2130 h2l[next_index++] = (2 + i) * n2 + n + 2 + j;
2131 // bottom face
2132 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2133 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2134 h2l[next_index++] = (2 + i) * n + 2 + j;
2135 // top face
2136 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2137 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2138 h2l[next_index++] = n2 + (2 + i) * n + 2 + j;
2139
2140 // inside hex
2141 Assert(fe.n_dofs_per_hex() ==
2142 fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(),
2144 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2145 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2146 for (unsigned int k = 0; k < fe.n_dofs_per_line(); ++k)
2147 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 2 + k;
2148
2149 Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2150
2151 break;
2152 }
2153
2154 default:
2156 }
2157 return h2l;
2158}
2159
2160
2161template <int dim>
2162std::vector<unsigned int>
2164 const unsigned int degree)
2165{
2166 FiniteElementData<dim - 1> fe_data(
2168 return internal::FE_Q_Hierarchical::invert_numbering(
2170 fe_data));
2171}
2172
2173
2174
2175template <>
2176std::vector<unsigned int>
2178 const unsigned int)
2179{
2180 return std::vector<unsigned int>();
2181}
2182
2183
2184template <>
2185bool
2187 const unsigned int face_index) const
2188{
2189 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2191
2192
2193 // in 1d, things are simple. since
2194 // there is only one degree of
2195 // freedom per vertex in this
2196 // class, the first is on vertex 0
2197 // (==face 0 in some sense), the
2198 // second on face 1:
2199 return (((shape_index == 0) && (face_index == 0)) ||
2200 ((shape_index == 1) && (face_index == 1)));
2201}
2202
2203
2204
2205template <int dim>
2206bool
2208 const unsigned int face_index) const
2209{
2210 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2212
2213 // first, special-case interior
2214 // shape functions, since they
2215 // have no support no-where on
2216 // the boundary
2217 if (((dim == 2) && (shape_index >=
2218 this->get_first_quad_index(0 /*only one quad in 2d*/))) ||
2219 ((dim == 3) && (shape_index >= this->get_first_hex_index())))
2220 return false;
2221
2222 // let's see whether this is a
2223 // vertex
2224 if (shape_index < this->get_first_line_index())
2225 {
2226 // for Q elements, there is
2227 // one dof per vertex, so
2228 // shape_index==vertex_number. check
2229 // whether this vertex is
2230 // on the given face.
2231 const unsigned int vertex_no = shape_index;
2234 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
2235 if (GeometryInfo<dim>::face_to_cell_vertices(face_index, i) ==
2236 vertex_no)
2237 return true;
2238 return false;
2239 }
2240 else if (shape_index < this->get_first_quad_index(0))
2241 // ok, dof is on a line
2242 {
2243 const unsigned int line_index =
2244 (shape_index - this->get_first_line_index()) / this->n_dofs_per_line();
2247
2248 for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_face; ++i)
2249 if (GeometryInfo<dim>::face_to_cell_lines(face_index, i) == line_index)
2250 return true;
2251 return false;
2252 }
2253 else if (shape_index < this->get_first_hex_index())
2254 // dof is on a quad
2255 {
2256 const unsigned int quad_index =
2257 (shape_index - this->get_first_quad_index(0 /*first quad*/)) /
2258 this->n_dofs_per_quad(face_index);
2259 Assert(static_cast<signed int>(quad_index) <
2260 static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
2262
2263 // in 2d, cell bubble are
2264 // zero on all faces. but
2265 // we have treated this
2266 // case above already
2267 Assert(dim != 2, ExcInternalError());
2268
2269 // in 3d,
2270 // quad_index=face_index
2271 if (dim == 3)
2272 return (quad_index == face_index);
2273 else
2275 }
2276 else
2277 // dof on hex
2278 {
2279 // can only happen in 3d, but
2280 // this case has already been
2281 // covered above
2283 return false;
2284 }
2285
2286 // we should not have gotten here
2288 return false;
2289}
2290
2291
2292
2293template <int dim>
2294std::vector<unsigned int>
2296{
2297 Assert((sub_degree > 0) && (sub_degree <= this->degree),
2298 ExcIndexRange(sub_degree, 1, this->degree));
2299
2300 if (dim == 1)
2301 {
2302 std::vector<unsigned int> embedding_dofs(sub_degree + 1);
2303 for (unsigned int i = 0; i < (sub_degree + 1); ++i)
2304 embedding_dofs[i] = i;
2305
2306 return embedding_dofs;
2307 }
2308
2309 if (sub_degree == 1)
2310 {
2311 std::vector<unsigned int> embedding_dofs(
2313 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2314 embedding_dofs[i] = i;
2315
2316 return embedding_dofs;
2317 }
2318 else if (sub_degree == this->degree)
2319 {
2320 std::vector<unsigned int> embedding_dofs(this->n_dofs_per_cell());
2321 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2322 embedding_dofs[i] = i;
2323
2324 return embedding_dofs;
2325 }
2326
2327 if ((dim == 2) || (dim == 3))
2328 {
2329 std::vector<unsigned int> embedding_dofs(
2330 (dim == 2) ? (sub_degree + 1) * (sub_degree + 1) :
2331 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2332
2333 for (unsigned int i = 0;
2334 i < ((dim == 2) ?
2335 (sub_degree + 1) * (sub_degree + 1) :
2336 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2337 ++i)
2338 {
2339 // vertex
2341 embedding_dofs[i] = i;
2342 // line
2345 {
2346 const unsigned int j =
2348 const unsigned int line =
2350 (sub_degree - 1);
2351
2353 line * (this->degree - 1) + j;
2354 }
2355 // quad
2359 (sub_degree - 1))
2360 {
2361 const unsigned int j =
2364 (sub_degree - 1);
2365 const unsigned int k =
2368 (sub_degree - 1)) %
2369 (sub_degree - 1);
2370 const unsigned int face =
2373 k * (sub_degree - 1) - j) /
2374 ((sub_degree - 1) * (sub_degree - 1));
2375
2376 embedding_dofs[i] =
2378 GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2379 face * (this->degree - 1) * (this->degree - 1) +
2380 k * (this->degree - 1) + j;
2381 }
2382 // hex
2386 (sub_degree - 1) +
2388 (sub_degree - 1) * (sub_degree - 1))
2389 {
2390 const unsigned int j =
2394 (sub_degree - 1)) %
2395 (sub_degree - 1);
2396 const unsigned int k =
2400 (sub_degree - 1) -
2401 j) /
2402 (sub_degree - 1)) %
2403 (sub_degree - 1);
2404 const unsigned int l =
2408 (sub_degree - 1) -
2409 j - k * (sub_degree - 1)) /
2410 ((sub_degree - 1) * (sub_degree - 1));
2411
2412 embedding_dofs[i] =
2414 GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2415 GeometryInfo<dim>::quads_per_cell * (this->degree - 1) *
2416 (this->degree - 1) +
2417 l * (this->degree - 1) * (this->degree - 1) +
2418 k * (this->degree - 1) + j;
2419 }
2420 }
2421
2422 return embedding_dofs;
2423 }
2424 else
2425 {
2427 return std::vector<unsigned int>();
2428 }
2429}
2430
2431
2432
2433template <int dim>
2434std::pair<Table<2, bool>, std::vector<unsigned int>>
2436{
2437 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
2438 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2439 constant_modes(0, i) = true;
2440 for (unsigned int i = GeometryInfo<dim>::vertices_per_cell;
2441 i < this->n_dofs_per_cell();
2442 ++i)
2443 constant_modes(0, i) = false;
2444 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
2445 constant_modes, std::vector<unsigned int>(1, 0));
2446}
2447
2448
2449
2450template <int dim>
2451std::size_t
2457
2458
2459
2460// explicit instantiations
2461#include "fe/fe_q_hierarchical.inst"
2462
2463
const std::unique_ptr< ScalarPolynomialsBase< dim > > poly_space
Definition fe_poly.h:532
void build_dofs_cell(std::vector< FullMatrix< double > > &dofs_cell, std::vector< FullMatrix< double > > &dofs_subcell) const
virtual std::unique_ptr< FiniteElement< dim, dim > > clone() const override
virtual void get_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix) const override
virtual std::pair< Table< 2, bool >, std::vector< unsigned int > > get_constant_modes() const override
void initialize_embedding_and_restriction(const std::vector< FullMatrix< double > > &dofs_cell, const std::vector< FullMatrix< double > > &dofs_subcell)
static std::vector< unsigned int > get_dpo_vector(const unsigned int degree)
friend class FE_Q_Hierarchical
std::vector< unsigned int > get_embedding_dofs(const unsigned int sub_degree) const
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_quad_dof_identities(const FiniteElement< dim > &fe_other, const unsigned int face_no=0) const override
void initialize_generalized_support_points()
virtual FiniteElementDomination::Domination compare_for_domination(const FiniteElement< dim > &fe_other, const unsigned int codim=0) const override final
static std::vector< unsigned int > hierarchic_to_fe_q_hierarchical_numbering(const FiniteElementData< dim > &fe)
virtual const FullMatrix< double > & get_prolongation_matrix(const unsigned int child, const RefinementCase< dim > &refinement_case=RefinementCase< dim >::isotropic_refinement) const override
virtual void get_subface_interpolation_matrix(const FiniteElement< dim > &source, const unsigned int subface, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_vertex_dof_identities(const FiniteElement< dim > &fe_other) const override
static std::vector< unsigned int > face_fe_q_hierarchical_to_hierarchic_numbering(const unsigned int degree)
void initialize_generalized_face_support_points()
void initialize_constraints(const std::vector< FullMatrix< double > > &dofs_subcell)
virtual std::string get_name() const override
virtual void get_face_interpolation_matrix(const FiniteElement< dim > &source, FullMatrix< double > &matrix, const unsigned int face_no=0) const override
virtual std::size_t memory_consumption() const override
virtual bool has_support_on_face(const unsigned int shape_index, const unsigned int face_index) const override
virtual bool hp_constraints_are_implemented() const override
virtual std::vector< std::pair< unsigned int, unsigned int > > hp_line_dof_identities(const FiniteElement< dim > &fe_other) const override
unsigned int n_dofs_per_vertex() const
unsigned int n_dofs_per_line() const
friend class Tensor
Definition tensor.h:865
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > vertex_indices()