Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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{
75 TensorProductPolynomials<dim> *poly_space_derived_ptr =
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
111 initialize_constraints(dofs_subcell);
112 initialize_embedding_and_restriction(dofs_cell, dofs_subcell);
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
153 const FiniteElement<dim> &source,
154 FullMatrix<double> &matrix) const
155{
156 // support interpolation between FE_Q_Hierarchical only.
157 if (const FE_Q_Hierarchical<dim> *source_fe =
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 res.reserve(std::min(this_dpl, other_dpl));
279 for (unsigned int i = 0; i < std::min(this_dpl, other_dpl); ++i)
280 res.emplace_back(i, i);
281
282 return res;
283 }
284 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
285 {
286 // the FE_Nothing has no
287 // degrees of freedom, so there
288 // are no equivalencies to be
289 // recorded
290 return std::vector<std::pair<unsigned int, unsigned int>>();
291 }
292 else
293 {
295 return std::vector<std::pair<unsigned int, unsigned int>>();
296 }
297}
298
299template <int dim>
300std::vector<std::pair<unsigned int, unsigned int>>
302 const FiniteElement<dim> &fe_other,
303 const unsigned int) const
304{
305 // we can presently only compute
306 // these identities if both FEs are
307 // FE_Q_Hierarchicals or if the other
308 // one is an FE_Nothing.
309 if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
310 {
311 // TODO: the implementation makes the assumption that all faces have the
312 // same number of dofs
313 AssertDimension(this->n_unique_faces(), 1);
314 const unsigned int face_no = 0;
315
316 const unsigned int this_dpq = this->n_dofs_per_quad(face_no);
317 const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no);
318
319 // we deal with hierarchical 1d polynomials where dofs are enumerated
320 // increasingly. Thus we return a vector of pairs for the first N-1, where
321 // N is minimum number of dofs_per_line for each FE_Q_Hierarchical.
322 std::vector<std::pair<unsigned int, unsigned int>> res;
323 res.reserve(std::min(this_dpq, other_dpq));
324 for (unsigned int i = 0; i < std::min(this_dpq, other_dpq); ++i)
325 res.emplace_back(i, i);
326
327 return res;
328 }
329 else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
330 {
331 // the FE_Nothing has no
332 // degrees of freedom, so there
333 // are no equivalencies to be
334 // recorded
335 return std::vector<std::pair<unsigned int, unsigned int>>();
336 }
337 else
338 {
340 return std::vector<std::pair<unsigned int, unsigned int>>();
341 }
342}
343
344template <int dim>
347 const FiniteElement<dim> &fe_other,
348 const unsigned int codim) const
349{
350 Assert(codim <= dim, ExcImpossibleInDim(dim));
351
352 // vertex/line/face domination
353 // (if fe_other is derived from FE_DGQ)
354 // ------------------------------------
355 if (codim > 0)
356 if (dynamic_cast<const FE_DGQ<dim> *>(&fe_other) != nullptr)
357 // there are no requirements between continuous and discontinuous elements
359
360 // vertex/line/face domination
361 // (if fe_other is not derived from FE_DGQ)
362 // & cell domination
363 // ----------------------------------------
364 if (const FE_Q_Hierarchical<dim> *fe_hierarchical_other =
365 dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other))
366 {
367 if (this->degree < fe_hierarchical_other->degree)
369 else if (this->degree == fe_hierarchical_other->degree)
371 else
373 }
374 else if (const FE_Nothing<dim> *fe_nothing =
375 dynamic_cast<const FE_Nothing<dim> *>(&fe_other))
376 {
377 if (fe_nothing->is_dominating())
379 else
380 // the FE_Nothing has no degrees of freedom and it is typically used
381 // in a context where we don't require any continuity along the
382 // interface
384 }
385
388}
389
390
391//---------------------------------------------------------------------------
392// Auxiliary functions
393//---------------------------------------------------------------------------
394
395
396template <int dim>
397void
399 std::vector<FullMatrix<double>> &dofs_cell,
400 std::vector<FullMatrix<double>> &dofs_subcell) const
401{
402 const unsigned int dofs_1d =
403 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
404
405 // The dofs_subcell matrices are transposed
406 // (4.19), (4.21) and (4.27),(4.28),(4.30) in
407 // Demkowicz, Oden, Rachowicz, Hardy, CMAMAE 77, 79-112, 1989
408 // so that
409 // DoFs_c(j) = dofs_subcell[c](j,k) dofs_cell(k)
410
411 // TODO: The dofs_subcell shall differ by a factor 2^p due to shape functions
412 // defined on [0,1] instead of [-1,1]. However, that does not seem to be
413 // the case. Perhaps this factor is added later on in auxiliary functions
414 // which use these matrices.
415
416 // dofs_cells[0](j,k):
417 // 0 1 | 2 3 4.
418 // 0 1 0 | .
419 // 1 0 0 | .
420 // -------------------
421 // 2 \ .
422 // 3 \ 2^k * k! / (k-j)!j!
423 // 4 \ .
424
425 // dofs_subcells[0](j,k):
426 // 0 1 | 2 3 4 5 6 .
427 // 0 1 0 | .
428 // 1 1/2 1/2 | -1 0 -1 0 -1.
429 // -------------------------------
430 // 2 \ .
431 // 3 \ (-1)^(k+j)/ 2^k * k!/(k-j)!j!
432 // 4 \ .
433
434 // dofs_cells[1](j,k):
435 // 0 1 | 2 3 4.
436 // 0 0 0 | .
437 // 1 0 1 | .
438 // -------------------
439 // 2 \ .
440 // 3 \ (-1)^(k+j) * 2^k * k!/(k-j)!j!
441 // 4 \.
442
443 // dofs_subcells[1](j,k):
444 // 0 1 | 2 3 4 5 6 .
445 // 0 1/2 1/2 | -1 0 -1 0 -1.
446 // 1 0 1 | .
447 // -------------------------------
448 // 2 \ .
449 // 3 \ 1/ 2^k * k!/(k-j)!j!
450 // 4 .
451
452 for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
453 for (unsigned int j = 0; j < dofs_1d; ++j)
454 for (unsigned int k = 0; k < dofs_1d; ++k)
455 {
456 // upper diagonal block
457 if ((j <= 1) && (k <= 1))
458 {
459 if (((c == 0) && (j == 0) && (k == 0)) ||
460 ((c == 1) && (j == 1) && (k == 1)))
461 dofs_cell[c](j, k) = 1.;
462 else
463 dofs_cell[c](j, k) = 0.;
464
465 if (((c == 0) && (j == 1)) || ((c == 1) && (j == 0)))
466 dofs_subcell[c](j, k) = .5;
467 else if (((c == 0) && (k == 0)) || ((c == 1) && (k == 1)))
468 dofs_subcell[c](j, k) = 1.;
469 else
470 dofs_subcell[c](j, k) = 0.;
471 }
472 // upper right block
473 else if ((j <= 1) && (k >= 2))
474 {
475 if (((c == 0) && (j == 1) && ((k % 2) == 0)) ||
476 ((c == 1) && (j == 0) && ((k % 2) == 0)))
477 dofs_subcell[c](j, k) = -1.;
478 }
479 // upper diagonal block
480 else if ((j >= 2) && (k >= 2) && (j <= k))
481 {
482 double factor = 1.;
483 for (unsigned int i = 1; i <= j; ++i)
484 factor *=
485 (static_cast<double>(k - i + 1)) / (static_cast<double>(i));
486 // factor == k * (k-1) * ... * (k-j+1) / j! = k! / (k-j)! / j!
487 if (c == 0)
488 {
489 dofs_subcell[c](j, k) = ((k + j) % 2 == 0) ?
490 Utilities::pow(.5, k) * factor :
491 -Utilities::pow(.5, k) * factor;
492 dofs_cell[c](j, k) = Utilities::pow(2, j) * factor;
493 }
494 else
495 {
496 dofs_subcell[c](j, k) = Utilities::pow(.5, k) * factor;
497 dofs_cell[c](j, k) = ((k + j) % 2 == 0) ?
498 Utilities::pow(2, j) * factor :
499 -Utilities::pow(2, j) * factor;
500 }
501 }
502 }
503}
504
505
506
507template <int dim>
508void
510 const std::vector<FullMatrix<double>> &dofs_subcell)
511{
512 const unsigned int dofs_1d =
513 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
514
515 this->interface_constraints.TableBase<2, double>::reinit(
516 this->interface_constraints_size());
517
518 switch (dim)
519 {
520 case 1:
521 {
522 // no constraints in 1d
523 break;
524 }
525
526 case 2:
527 {
528 // vertex node
529 for (unsigned int i = 0; i < dofs_1d; ++i)
530 this->interface_constraints(0, i) = dofs_subcell[0](1, i);
531 // edge nodes
532 for (unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
533 ++c)
534 for (unsigned int i = 0; i < dofs_1d; ++i)
535 for (unsigned int j = 2; j < dofs_1d; ++j)
536 this->interface_constraints(1 + c * (this->degree - 1) + j - 2,
537 i) = dofs_subcell[c](j, i);
538 break;
539 }
540
541 case 3:
542 {
543 for (unsigned int i = 0; i < dofs_1d * dofs_1d; ++i)
544 {
545 // center vertex node
546 this->interface_constraints(0, face_renumber[i]) =
547 dofs_subcell[0](1, i % dofs_1d) *
548 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
549
550 // boundary vertex nodes
551 this->interface_constraints(1, face_renumber[i]) =
552 dofs_subcell[0](0, i % dofs_1d) *
553 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
554 this->interface_constraints(2, face_renumber[i]) =
555 dofs_subcell[1](1, i % dofs_1d) *
556 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
557 this->interface_constraints(3, face_renumber[i]) =
558 dofs_subcell[0](1, i % dofs_1d) *
559 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
560 this->interface_constraints(4, face_renumber[i]) =
561 dofs_subcell[1](0, i % dofs_1d) *
562 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
563
564 // interior edges
565 for (unsigned int j = 0; j < (this->degree - 1); ++j)
566 {
567 this->interface_constraints(5 + j, face_renumber[i]) =
568 dofs_subcell[0](1, i % dofs_1d) *
569 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
570 this->interface_constraints(5 + (this->degree - 1) + j,
571 face_renumber[i]) =
572 dofs_subcell[0](1, i % dofs_1d) *
573 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
574 this->interface_constraints(5 + 2 * (this->degree - 1) + j,
575 face_renumber[i]) =
576 dofs_subcell[0](2 + j, i % dofs_1d) *
577 dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
578 this->interface_constraints(5 + 3 * (this->degree - 1) + j,
579 face_renumber[i]) =
580 dofs_subcell[1](2 + j, i % dofs_1d) *
581 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
582 }
583
584 // boundary edges
585 for (unsigned int j = 0; j < (this->degree - 1); ++j)
586 {
587 // left edge
588 this->interface_constraints(5 + 4 * (this->degree - 1) + j,
589 face_renumber[i]) =
590 dofs_subcell[0](0, i % dofs_1d) *
591 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
592 this->interface_constraints(5 + 4 * (this->degree - 1) +
593 (this->degree - 1) + j,
594 face_renumber[i]) =
595 dofs_subcell[0](0, i % dofs_1d) *
596 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
597 // right edge
598 this->interface_constraints(5 + 4 * (this->degree - 1) +
599 2 * (this->degree - 1) + j,
600 face_renumber[i]) =
601 dofs_subcell[1](1, i % dofs_1d) *
602 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
603 this->interface_constraints(5 + 4 * (this->degree - 1) +
604 3 * (this->degree - 1) + j,
605 face_renumber[i]) =
606 dofs_subcell[1](1, i % dofs_1d) *
607 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
608 // bottom edge
609 this->interface_constraints(5 + 4 * (this->degree - 1) +
610 4 * (this->degree - 1) + j,
611 face_renumber[i]) =
612 dofs_subcell[0](2 + j, i % dofs_1d) *
613 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
614 this->interface_constraints(5 + 4 * (this->degree - 1) +
615 5 * (this->degree - 1) + j,
616 face_renumber[i]) =
617 dofs_subcell[1](2 + j, i % dofs_1d) *
618 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
619 // top edge
620 this->interface_constraints(5 + 4 * (this->degree - 1) +
621 6 * (this->degree - 1) + j,
622 face_renumber[i]) =
623 dofs_subcell[0](2 + j, i % dofs_1d) *
624 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
625 this->interface_constraints(5 + 4 * (this->degree - 1) +
626 7 * (this->degree - 1) + j,
627 face_renumber[i]) =
628 dofs_subcell[1](2 + j, i % dofs_1d) *
629 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
630 }
631
632 // interior faces
633 for (unsigned int j = 0; j < (this->degree - 1); ++j)
634 for (unsigned int k = 0; k < (this->degree - 1); ++k)
635 {
636 // subcell 0
637 this->interface_constraints(5 + 12 * (this->degree - 1) +
638 j + k * (this->degree - 1),
639 face_renumber[i]) =
640 dofs_subcell[0](2 + j, i % dofs_1d) *
641 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
642 // subcell 1
643 this->interface_constraints(5 + 12 * (this->degree - 1) +
644 j + k * (this->degree - 1) +
645 (this->degree - 1) *
646 (this->degree - 1),
647 face_renumber[i]) =
648 dofs_subcell[1](2 + j, i % dofs_1d) *
649 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
650 // subcell 2
651 this->interface_constraints(5 + 12 * (this->degree - 1) +
652 j + k * (this->degree - 1) +
653 2 * (this->degree - 1) *
654 (this->degree - 1),
655 face_renumber[i]) =
656 dofs_subcell[0](2 + j, i % dofs_1d) *
657 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
658 // subcell 3
659 this->interface_constraints(5 + 12 * (this->degree - 1) +
660 j + k * (this->degree - 1) +
661 3 * (this->degree - 1) *
662 (this->degree - 1),
663 face_renumber[i]) =
664 dofs_subcell[1](2 + j, i % dofs_1d) *
665 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
666 }
667 }
668 break;
669 }
670
671 default:
673 }
674}
675
676
677
678template <int dim>
679void
681 const std::vector<FullMatrix<double>> &dofs_cell,
682 const std::vector<FullMatrix<double>> &dofs_subcell)
683{
684 unsigned int iso = RefinementCase<dim>::isotropic_refinement - 1;
685
686 const unsigned int dofs_1d =
687 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
688 TensorProductPolynomials<dim> *poly_space_derived_ptr =
689 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
690 const std::vector<unsigned int> &renumber =
691 poly_space_derived_ptr->get_numbering();
692
693 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
694 {
695 this->prolongation[iso][c].reinit(this->n_dofs_per_cell(),
696 this->n_dofs_per_cell());
697 this->restriction[iso][c].reinit(this->n_dofs_per_cell(),
698 this->n_dofs_per_cell());
699 }
700
701 // the 1d case is particularly
702 // simple, so special case it:
703 if (dim == 1)
704 {
705 for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
706 ++c)
707 {
708 this->prolongation[iso][c].fill(dofs_subcell[c]);
709 this->restriction[iso][c].fill(dofs_cell[c]);
710 }
711 return;
712 }
713
714 // for higher dimensions, things
715 // are a little more tricky:
716
717 // j loops over dofs in the
718 // subcell. These are the rows in
719 // the embedding matrix.
720 //
721 // i loops over the dofs in the
722 // parent cell. These are the
723 // columns in the embedding matrix.
724 for (unsigned int j = 0; j < this->n_dofs_per_cell(); ++j)
725 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
726 switch (dim)
727 {
728 case 2:
729 {
730 for (unsigned int c = 0;
731 c < GeometryInfo<2>::max_children_per_cell;
732 ++c)
733 {
734 // left/right line: 0/1
735 const unsigned int c0 = c % 2;
736 // bottom/top line: 0/1
737 const unsigned int c1 = c / 2;
738
739 this->prolongation[iso][c](j, i) =
740 dofs_subcell[c0](renumber[j] % dofs_1d,
741 renumber[i] % dofs_1d) *
742 dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
743 dofs_1d,
744 (renumber[i] - (renumber[i] % dofs_1d)) /
745 dofs_1d);
746
747 this->restriction[iso][c](j, i) =
748 dofs_cell[c0](renumber[j] % dofs_1d,
749 renumber[i] % dofs_1d) *
750 dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
751 dofs_1d,
752 (renumber[i] - (renumber[i] % dofs_1d)) /
753 dofs_1d);
754 }
755 break;
756 }
757
758 case 3:
759 {
760 for (unsigned int c = 0;
761 c < GeometryInfo<3>::max_children_per_cell;
762 ++c)
763 {
764 // left/right face: 0/1
765 const unsigned int c0 = c % 2;
766 // front/back face: 0/1
767 const unsigned int c1 = (c % 4) / 2;
768 // bottom/top face: 0/1
769 const unsigned int c2 = c / 4;
770
771 this->prolongation[iso][c](j, i) =
772 dofs_subcell[c0](renumber[j] % dofs_1d,
773 renumber[i] % dofs_1d) *
774 dofs_subcell[c1](
775 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
776 dofs_1d,
777 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
778 dofs_1d) *
779 dofs_subcell[c2](
780 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
781 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
782 dofs_1d)) /
783 dofs_1d,
784 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
785 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
786 dofs_1d)) /
787 dofs_1d);
788
789 this->restriction[iso][c](j, i) =
790 dofs_cell[c0](renumber[j] % dofs_1d,
791 renumber[i] % dofs_1d) *
792 dofs_cell[c1](
793 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
794 dofs_1d,
795 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
796 dofs_1d) *
797 dofs_cell[c2](
798 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
799 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
800 dofs_1d)) /
801 dofs_1d,
802 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
803 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
804 dofs_1d)) /
805 dofs_1d);
806 }
807 break;
808 }
809
810 default:
812 }
813}
814
815
816
817template <int dim>
818void
820{
821 // number of points: (degree+1)^dim
822 unsigned int n = this->degree + 1;
823 for (unsigned int i = 1; i < dim; ++i)
824 n *= this->degree + 1;
825
826 this->generalized_support_points.resize(n);
827
828 TensorProductPolynomials<dim> *poly_space_derived_ptr =
829 dynamic_cast<TensorProductPolynomials<dim> *>(this->poly_space.get());
830 const std::vector<unsigned int> &index_map_inverse =
831 poly_space_derived_ptr->get_numbering_inverse();
832
833 Point<dim> p;
834 // the method of numbering allows
835 // each dof to be associated with a
836 // support point. There is
837 // only one support point per
838 // vertex, line, quad, hex, etc.
839 //
840 // note, however, that the support
841 // points thus associated with
842 // shape functions are not unique:
843 // the linear shape functions are
844 // associated with the vertices,
845 // but all others are associated
846 // with either line, quad, or hex
847 // midpoints, and there may be
848 // multiple shape functions
849 // associated with them. there
850 // really is no other useful
851 // numbering, since the
852 // hierarchical shape functions do
853 // not vanish at all-but-one
854 // interpolation points (like the
855 // Lagrange functions used in
856 // FE_Q), so there's not much we
857 // can do here.
858
859 // TODO shouldn't we just at least make support points unique,
860 // even though the delta property is not satisfied for this FE?
861 unsigned int k = 0;
862 for (unsigned int iz = 0; iz <= ((dim > 2) ? this->degree : 0); ++iz)
863 for (unsigned int iy = 0; iy <= ((dim > 1) ? this->degree : 0); ++iy)
864 for (unsigned int ix = 0; ix <= this->degree; ++ix)
865 {
866 if (ix == 0)
867 p[0] = 0.;
868 else if (ix == 1)
869 p[0] = 1.;
870 else
871 p[0] = .5;
872 if (dim > 1)
873 {
874 if (iy == 0)
875 p[1] = 0.;
876 else if (iy == 1)
877 p[1] = 1.;
878 else
879 p[1] = .5;
880 }
881 if (dim > 2)
882 {
883 if (iz == 0)
884 p[2] = 0.;
885 else if (iz == 1)
886 p[2] = 1.;
887 else
888 p[2] = .5;
889 }
890 this->generalized_support_points[index_map_inverse[k++]] = p;
891 }
892}
893
894
895
896template <>
897void
899{
900 // no faces in 1d, so nothing to do
901}
902
903
904template <>
905void
907 const FiniteElement<1, 1> & /*x_source_fe*/,
908 FullMatrix<double> & /*interpolation_matrix*/,
909 const unsigned int) const
910{
911 Assert(false, ExcImpossibleInDim(1));
912}
913
914
915template <>
916void
918 const FiniteElement<1, 1> & /*x_source_fe*/,
919 const unsigned int /*subface*/,
920 FullMatrix<double> & /*interpolation_matrix*/,
921 const unsigned int) const
922{
923 Assert(false, ExcImpossibleInDim(1));
924}
925
926
927
928template <int dim>
929void
931 const FiniteElement<dim> &x_source_fe,
932 FullMatrix<double> &interpolation_matrix,
933 const unsigned int face_no) const
934{
935 // this is only implemented, if the
936 // source FE is also a
937 // Q_Hierarchical element
938 using FEQHierarchical = FE_Q_Hierarchical<dim>;
939 AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
940 (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
941 nullptr),
943
944 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
945 ExcDimensionMismatch(interpolation_matrix.n(),
946 this->n_dofs_per_face(face_no)));
947 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
948 ExcDimensionMismatch(interpolation_matrix.m(),
949 x_source_fe.n_dofs_per_face(face_no)));
950
951 // ok, source is a Q_Hierarchical element, so
952 // we will be able to do the work
953 const FE_Q_Hierarchical<dim> &source_fe =
954 dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
955 (void)source_fe;
956
957 // Make sure, that the element,
958 // for which the DoFs should be
959 // constrained is the one with
960 // the higher polynomial degree.
961 // Actually the procedure will work
962 // also if this assertion is not
963 // satisfied. But the matrices
964 // produced in that case might
965 // lead to problems in the
966 // hp-procedures, which use this
967 // method.
968 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
970 interpolation_matrix = 0;
971
972 switch (dim)
973 {
974 case 2:
975 {
976 // In 2-dimension the constraints are trivial.
977 // First this->dofs_per_face DoFs of the constrained
978 // element are made equal to the current (dominating)
979 // element, which corresponds to 1 on diagonal of the matrix.
980 // DoFs which correspond to higher polynomials
981 // are zeroed (zero rows in the matrix).
982 for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
983 interpolation_matrix(i, i) = 1;
984
985 break;
986 }
987
988 case 3:
989 {
990 for (unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i)
991 interpolation_matrix(i, i) = 1;
992
993 for (unsigned int i = 0; i < this->degree - 1; ++i)
994 {
995 for (unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j)
996 interpolation_matrix(i + j * (x_source_fe.degree - 1) +
998 i + j * (this->degree - 1) +
1000
1001 for (unsigned int j = 0; j < this->degree - 1; ++j)
1002 interpolation_matrix((i + GeometryInfo<3>::lines_per_face) *
1003 (x_source_fe.degree - 1) +
1006 (this->degree - 1) +
1008 1;
1009 }
1010 }
1011 }
1012}
1013
1014
1015
1016template <int dim>
1017void
1019 const FiniteElement<dim> &x_source_fe,
1020 const unsigned int subface,
1021 FullMatrix<double> &interpolation_matrix,
1022 const unsigned int face_no) const
1023{
1024 // this is only implemented, if the
1025 // source FE is also a
1026 // Q_Hierarchical element
1027 using FEQHierarchical = FE_Q_Hierarchical<dim>;
1028 AssertThrow((x_source_fe.get_name().find("FE_Q_Hierarchical<") == 0) ||
1029 (dynamic_cast<const FEQHierarchical *>(&x_source_fe) !=
1030 nullptr),
1032
1033 Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
1034 ExcDimensionMismatch(interpolation_matrix.n(),
1035 this->n_dofs_per_face(face_no)));
1036 Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
1037 ExcDimensionMismatch(interpolation_matrix.m(),
1038 x_source_fe.n_dofs_per_face(face_no)));
1039
1040 // ok, source is a Q_Hierarchical element, so
1041 // we will be able to do the work
1042 const FE_Q_Hierarchical<dim> &source_fe =
1043 dynamic_cast<const FE_Q_Hierarchical<dim> &>(x_source_fe);
1044
1045 // Make sure, that the element,
1046 // for which the DoFs should be
1047 // constrained is the one with
1048 // the higher polynomial degree.
1049 // Actually the procedure will work
1050 // also if this assertion is not
1051 // satisfied. But the matrices
1052 // produced in that case might
1053 // lead to problems in the
1054 // hp-procedures, which use this
1055 // method.
1056 Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
1058
1059 switch (dim)
1060 {
1061 case 2:
1062 {
1063 switch (subface)
1064 {
1065 case 0:
1066 {
1067 interpolation_matrix(0, 0) = 1.0;
1068 interpolation_matrix(1, 0) = 0.5;
1069 interpolation_matrix(1, 1) = 0.5;
1070
1071 for (unsigned int dof = 2;
1072 dof < this->n_dofs_per_face(face_no);)
1073 {
1074 interpolation_matrix(1, dof) = -1.0;
1075 dof = dof + 2;
1076 }
1077
1078 int factorial_i = 1;
1079
1080 for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1081 ++i)
1082 {
1083 interpolation_matrix(i, i) = Utilities::pow(0.5, i);
1084 factorial_i *= i;
1085 int factorial_j = factorial_i;
1086 int factorial_ij = 1;
1087
1088 for (unsigned int j = i + 1;
1089 j < this->n_dofs_per_face(face_no);
1090 ++j)
1091 {
1092 factorial_ij *= j - i;
1093 factorial_j *= j;
1094
1095 if (((i + j) & 1) != 0u)
1096 interpolation_matrix(i, j) =
1097 -1.0 * Utilities::pow(0.5, j) * factorial_j /
1098 (factorial_i * factorial_ij);
1099
1100 else
1101 interpolation_matrix(i, j) =
1102 Utilities::pow(0.5, j) * factorial_j /
1103 (factorial_i * factorial_ij);
1104 }
1105 }
1106
1107 break;
1108 }
1109
1110 case 1:
1111 {
1112 interpolation_matrix(0, 0) = 0.5;
1113 interpolation_matrix(0, 1) = 0.5;
1114
1115 for (unsigned int dof = 2;
1116 dof < this->n_dofs_per_face(face_no);)
1117 {
1118 interpolation_matrix(0, dof) = -1.0;
1119 dof = dof + 2;
1120 }
1121
1122 interpolation_matrix(1, 1) = 1.0;
1123
1124 int factorial_i = 1;
1125
1126 for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1127 ++i)
1128 {
1129 interpolation_matrix(i, i) = Utilities::pow(0.5, i);
1130 factorial_i *= i;
1131 int factorial_j = factorial_i;
1132 int factorial_ij = 1;
1133
1134 for (unsigned int j = i + 1;
1135 j < this->n_dofs_per_face(face_no);
1136 ++j)
1137 {
1138 factorial_ij *= j - i;
1139 factorial_j *= j;
1140 interpolation_matrix(i, j) =
1141 Utilities::pow(0.5, j) * factorial_j /
1142 (factorial_i * factorial_ij);
1143 }
1144 }
1145 }
1146 }
1147
1148 break;
1149 }
1150
1151 case 3:
1152 {
1153 switch (subface)
1154 {
1155 case 0:
1156 {
1157 interpolation_matrix(0, 0) = 1.0;
1158 interpolation_matrix(1, 0) = 0.5;
1159 interpolation_matrix(1, 1) = 0.5;
1160 interpolation_matrix(2, 0) = 0.5;
1161 interpolation_matrix(2, 2) = 0.5;
1162
1163 for (unsigned int i = 0; i < this->degree - 1;)
1164 {
1165 for (unsigned int line = 0;
1166 line < GeometryInfo<3>::lines_per_face;
1167 ++line)
1168 interpolation_matrix(3,
1169 i + line * (this->degree - 1) +
1170 4) = -0.5;
1171
1172 for (unsigned int j = 0; j < this->degree - 1;)
1173 {
1174 interpolation_matrix(3,
1175 i + (j + 4) * this->degree - j) =
1176 1.0;
1177 j = j + 2;
1178 }
1179
1180 interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1181 -1.0;
1182 interpolation_matrix(2, i + 4) = -1.0;
1183 i = i + 2;
1184 }
1185
1186 for (unsigned int vertex = 0;
1187 vertex < GeometryInfo<3>::vertices_per_face;
1188 ++vertex)
1189 interpolation_matrix(3, vertex) = 0.25;
1190
1191 int factorial_i = 1;
1192
1193 for (unsigned int i = 2; i <= this->degree; ++i)
1194 {
1195 double tmp = Utilities::pow(0.5, i);
1196 interpolation_matrix(i + 2, i + 2) = tmp;
1197 interpolation_matrix(i + 2 * source_fe.degree,
1198 i + 2 * this->degree) = tmp;
1199 tmp *= 0.5;
1200 interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1201 tmp;
1202 interpolation_matrix(i + source_fe.degree + 1,
1203 i + this->degree + 1) = tmp;
1204 interpolation_matrix(i + 3 * source_fe.degree - 1,
1205 i + 2 * this->degree) = tmp;
1206 interpolation_matrix(i + 3 * source_fe.degree - 1,
1207 i + 3 * this->degree - 1) = tmp;
1208 tmp *= -2.0;
1209
1210 for (unsigned int j = 0; j < this->degree - 1;)
1211 {
1212 interpolation_matrix(i + source_fe.degree + 1,
1213 (i + 2) * this->degree + j + 2 -
1214 i) = tmp;
1215 interpolation_matrix(i + 3 * source_fe.degree - 1,
1216 i + (j + 4) * this->degree - j -
1217 2) = tmp;
1218 j = j + 2;
1219 }
1220
1221 int factorial_k = 1;
1222
1223 for (unsigned int j = 2; j <= this->degree; ++j)
1224 {
1225 interpolation_matrix(i + (j + 2) * source_fe.degree -
1226 j,
1227 i + (j + 2) * this->degree - j) =
1228 Utilities::pow(0.5, i + j);
1229 factorial_k *= j;
1230 int factorial_kl = 1;
1231 int factorial_l = factorial_k;
1232
1233 for (unsigned int k = j + 1; k < this->degree; ++k)
1234 {
1235 factorial_kl *= k - j;
1236 factorial_l *= k;
1237
1238 if (((j + k) & 1) != 0u)
1239 interpolation_matrix(
1240 i + (j + 2) * source_fe.degree - j,
1241 i + (k + 2) * this->degree - k) =
1242 -1.0 * Utilities::pow(0.5, i + k) *
1243 factorial_l / (factorial_k * factorial_kl);
1244
1245 else
1246 interpolation_matrix(
1247 i + (j + 2) * source_fe.degree - j,
1248 i + (k + 2) * this->degree - k) =
1249 Utilities::pow(0.5, i + k) * factorial_l /
1250 (factorial_k * factorial_kl);
1251 }
1252 }
1253
1254 factorial_i *= i;
1255 int factorial_j = factorial_i;
1256 int factorial_ij = 1;
1257
1258 for (unsigned int j = i + 1; j <= this->degree; ++j)
1259 {
1260 factorial_ij *= j - i;
1261 factorial_j *= j;
1262
1263 if (((i + j) & 1) != 0u)
1264 {
1265 tmp = -1.0 * Utilities::pow(0.5, j) *
1266 factorial_j / (factorial_i * factorial_ij);
1267 interpolation_matrix(i + 2, j + 2) = tmp;
1268 interpolation_matrix(i + 2 * source_fe.degree,
1269 j + 2 * this->degree) = tmp;
1270 factorial_k = 1;
1271
1272 for (unsigned int k = 2; k <= this->degree; ++k)
1273 {
1274 interpolation_matrix(
1275 i + (k + 2) * source_fe.degree - k,
1276 j + (k + 2) * this->degree - k) =
1277 tmp * Utilities::pow(0.5, k);
1278 factorial_k *= k;
1279 int factorial_l = factorial_k;
1280 int factorial_kl = 1;
1281
1282 for (unsigned int l = k + 1;
1283 l <= this->degree;
1284 ++l)
1285 {
1286 factorial_kl *= l - k;
1287 factorial_l *= l;
1288
1289 if (((k + l) & 1) != 0u)
1290 interpolation_matrix(
1291 i + (k + 2) * source_fe.degree - k,
1292 j + (l + 2) * this->degree - l) =
1293 -1.0 * tmp * Utilities::pow(0.5, l) *
1294 factorial_l /
1295 (factorial_k * factorial_kl);
1296
1297 else
1298 interpolation_matrix(
1299 i + (k + 2) * source_fe.degree - k,
1300 j + (l + 2) * this->degree - l) =
1301 tmp * Utilities::pow(0.5, l) *
1302 factorial_l /
1303 (factorial_k * factorial_kl);
1304 }
1305 }
1306
1307 tmp *= 0.5;
1308 interpolation_matrix(i + source_fe.degree + 1,
1309 j + 2) = tmp;
1310 interpolation_matrix(i + source_fe.degree + 1,
1311 j + this->degree + 1) = tmp;
1312 interpolation_matrix(i + 3 * source_fe.degree - 1,
1313 j + 2 * this->degree) = tmp;
1314 interpolation_matrix(i + 3 * source_fe.degree - 1,
1315 j + 3 * this->degree - 1) =
1316 tmp;
1317 tmp *= -2.0;
1318
1319 for (unsigned int k = 0; k < this->degree - 1;)
1320 {
1321 interpolation_matrix(i + source_fe.degree + 1,
1322 (j + 2) * this->degree +
1323 k + 2 - j) = tmp;
1324 interpolation_matrix(
1325 i + 3 * source_fe.degree - 1,
1326 j + (k + 4) * this->degree - k - 2) = tmp;
1327 k = k + 2;
1328 }
1329 }
1330 else
1331 {
1332 tmp = Utilities::pow(0.5, j) * factorial_j /
1333 (factorial_i * factorial_ij);
1334 interpolation_matrix(i + 2, j + 2) = tmp;
1335 interpolation_matrix(i + 2 * source_fe.degree,
1336 j + 2 * this->degree) = tmp;
1337 factorial_k = 1;
1338
1339 for (unsigned int k = 2; k <= this->degree; ++k)
1340 {
1341 interpolation_matrix(
1342 i + (k + 2) * source_fe.degree - k,
1343 j + (k + 2) * this->degree - k) =
1344 tmp * Utilities::pow(0.5, k);
1345 factorial_k *= k;
1346 int factorial_l = factorial_k;
1347 int factorial_kl = 1;
1348
1349 for (unsigned int l = k + 1;
1350 l <= this->degree;
1351 ++l)
1352 {
1353 factorial_kl *= l - k;
1354 factorial_l *= l;
1355
1356 if (((k + l) & 1) != 0u)
1357 interpolation_matrix(
1358 i + (k + 2) * source_fe.degree - k,
1359 j + (l + 2) * this->degree - l) =
1360 -1.0 * tmp * Utilities::pow(0.5, l) *
1361 factorial_l /
1362 (factorial_k * factorial_kl);
1363
1364 else
1365 interpolation_matrix(
1366 i + (k + 2) * source_fe.degree - k,
1367 j + (l + 2) * this->degree - l) =
1368 tmp * Utilities::pow(0.5, l) *
1369 factorial_l /
1370 (factorial_k * factorial_kl);
1371 }
1372 }
1373
1374 tmp *= 0.5;
1375 interpolation_matrix(i + source_fe.degree + 1,
1376 j + 2) = tmp;
1377 interpolation_matrix(i + source_fe.degree + 1,
1378 j + this->degree + 1) = tmp;
1379 interpolation_matrix(i + 3 * source_fe.degree - 1,
1380 j + 2 * this->degree) = tmp;
1381 interpolation_matrix(i + 3 * source_fe.degree - 1,
1382 j + 3 * this->degree - 1) =
1383 tmp;
1384 tmp *= -2.0;
1385
1386 for (unsigned int k = 0; k < this->degree - 1;)
1387 {
1388 interpolation_matrix(i + source_fe.degree + 1,
1389 (j + 2) * this->degree +
1390 k + 2 - j) = tmp;
1391 interpolation_matrix(
1392 i + 3 * source_fe.degree - 1,
1393 j + (k + 4) * this->degree - k - 2) = tmp;
1394 k = k + 2;
1395 }
1396 }
1397 }
1398 }
1399
1400 break;
1401 }
1402
1403 case 1:
1404 {
1405 interpolation_matrix(0, 0) = 0.5;
1406 interpolation_matrix(0, 1) = 0.5;
1407 interpolation_matrix(1, 1) = 1.0;
1408 interpolation_matrix(3, 1) = 0.5;
1409 interpolation_matrix(3, 3) = 0.5;
1410
1411 for (unsigned int i = 0; i < this->degree - 1;)
1412 {
1413 for (unsigned int line = 0;
1414 line < GeometryInfo<3>::lines_per_face;
1415 ++line)
1416 interpolation_matrix(2,
1417 i + line * (this->degree - 1) +
1418 4) = -0.5;
1419
1420 for (unsigned int j = 0; j < this->degree - 1;)
1421 {
1422 interpolation_matrix(2,
1423 i + (j + 4) * this->degree - j) =
1424 1.0;
1425 j = j + 2;
1426 }
1427
1428 interpolation_matrix(0, i + 2 * (this->degree + 1)) =
1429 -1.0;
1430 interpolation_matrix(3, i + this->degree + 3) = -1.0;
1431 i = i + 2;
1432 }
1433
1434 for (unsigned int vertex = 0;
1435 vertex < GeometryInfo<3>::vertices_per_face;
1436 ++vertex)
1437 interpolation_matrix(2, vertex) = 0.25;
1438
1439 int factorial_i = 1;
1440
1441 for (unsigned int i = 2; i <= this->degree; ++i)
1442 {
1443 double tmp = Utilities::pow(0.5, i + 1);
1444 interpolation_matrix(i + 2, i + 2) = tmp;
1445 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1446 interpolation_matrix(i + 3 * source_fe.degree - 1,
1447 i + 2 * this->degree) = tmp;
1448 interpolation_matrix(i + 3 * source_fe.degree - 1,
1449 i + 3 * this->degree - 1) = tmp;
1450 tmp *= -2.0;
1451
1452 for (unsigned int j = 0; j < this->degree - 1;)
1453 {
1454 interpolation_matrix(i + 2,
1455 j + (i + 2) * this->degree + 2 -
1456 i) = tmp;
1457 interpolation_matrix(i + 3 * source_fe.degree - 1,
1458 i + (j + 4) * this->degree - j -
1459 2) = tmp;
1460 j = j + 2;
1461 }
1462
1463 tmp *= -1.0;
1464 interpolation_matrix(i + source_fe.degree + 1,
1465 i + this->degree + 1) = tmp;
1466 interpolation_matrix(i + 2 * source_fe.degree,
1467 i + 2 * this->degree) = tmp;
1468 factorial_i *= i;
1469 int factorial_j = factorial_i;
1470 int factorial_ij = 1;
1471
1472 for (unsigned int j = i + 1; j <= this->degree; ++j)
1473 {
1474 factorial_ij *= j - i;
1475 factorial_j *= j;
1476 tmp = Utilities::pow(0.5, j) * factorial_j /
1477 (factorial_i * factorial_ij);
1478 interpolation_matrix(i + 2 * source_fe.degree,
1479 j + 2 * this->degree) = tmp;
1480 int factorial_k = 1;
1481
1482 for (unsigned int k = 2; k <= this->degree; ++k)
1483 {
1484 interpolation_matrix(
1485 i + (k + 2) * source_fe.degree - k,
1486 j + (k + 2) * this->degree - k) =
1487 tmp * Utilities::pow(0.5, k);
1488 factorial_k *= k;
1489 int factorial_l = factorial_k;
1490 int factorial_kl = 1;
1491
1492 for (unsigned int l = k + 1; l <= this->degree;
1493 ++l)
1494 {
1495 factorial_kl *= l - k;
1496 factorial_l *= l;
1497
1498 if (((k + l) & 1) != 0u)
1499 interpolation_matrix(
1500 i + (k + 2) * source_fe.degree - k,
1501 j + (l + 2) * this->degree - l) =
1502 -1.0 * tmp * Utilities::pow(0.5, l) *
1503 factorial_l /
1504 (factorial_k * factorial_kl);
1505
1506 else
1507 interpolation_matrix(
1508 i + (k + 2) * source_fe.degree - k,
1509 j + (l + 2) * this->degree - l) =
1510 tmp * Utilities::pow(0.5, l) *
1511 factorial_l /
1512 (factorial_k * factorial_kl);
1513 }
1514 }
1515
1516 tmp *= -1.0;
1517
1518 for (unsigned int k = 0; k < this->degree - 1;)
1519 {
1520 interpolation_matrix(i + 3 * source_fe.degree - 1,
1521 j + (k + 4) * this->degree -
1522 k - 2) = tmp;
1523 k = k + 2;
1524 }
1525
1526 tmp *= -0.5;
1527 interpolation_matrix(i + 3 * source_fe.degree - 1,
1528 j + 2 * this->degree) = tmp;
1529 interpolation_matrix(i + 3 * source_fe.degree - 1,
1530 j + 3 * this->degree - 1) = tmp;
1531
1532 if (((i + j) & 1) != 0u)
1533 tmp *= -1.0;
1534
1535 interpolation_matrix(i + 2, j + 2) = tmp;
1536 interpolation_matrix(i + 2, j + this->degree + 1) =
1537 tmp;
1538 interpolation_matrix(i + source_fe.degree + 1,
1539 j + this->degree + 1) =
1540 2.0 * tmp;
1541 tmp *= -2.0;
1542
1543 for (unsigned int k = 0; k < this->degree - 1;)
1544 {
1545 interpolation_matrix(i + 2,
1546 k + (j + 2) * this->degree +
1547 2 - j) = tmp;
1548 k = k + 2;
1549 }
1550 }
1551
1552 int factorial_k = 1;
1553
1554 for (unsigned int j = 2; j <= this->degree; ++j)
1555 {
1556 interpolation_matrix(i + (j + 2) * source_fe.degree -
1557 j,
1558 i + (j + 2) * this->degree - j) =
1559 Utilities::pow(0.5, i + j);
1560 factorial_k *= j;
1561 int factorial_l = factorial_k;
1562 int factorial_kl = 1;
1563
1564 for (unsigned int k = j + 1; k <= this->degree; ++k)
1565 {
1566 factorial_kl *= k - j;
1567 factorial_l *= k;
1568
1569 if (((j + k) & 1) != 0u)
1570 interpolation_matrix(
1571 i + (j + 2) * source_fe.degree - j,
1572 i + (k + 2) * this->degree - k) =
1573 -1.0 * Utilities::pow(0.5, i + k) *
1574 factorial_l / (factorial_k * factorial_kl);
1575
1576 else
1577 interpolation_matrix(
1578 i + (j + 2) * source_fe.degree - j,
1579 i + (k + 2) * this->degree - k) =
1580 Utilities::pow(0.5, i + k) * factorial_l /
1581 (factorial_k * factorial_kl);
1582 }
1583 }
1584 }
1585
1586 break;
1587 }
1588
1589 case 2:
1590 {
1591 interpolation_matrix(0, 0) = 0.5;
1592 interpolation_matrix(0, 2) = 0.5;
1593 interpolation_matrix(2, 2) = 1.0;
1594 interpolation_matrix(3, 2) = 0.5;
1595 interpolation_matrix(3, 3) = 0.5;
1596
1597 for (unsigned int i = 0; i < this->degree - 1;)
1598 {
1599 for (unsigned int line = 0;
1600 line < GeometryInfo<3>::lines_per_face;
1601 ++line)
1602 interpolation_matrix(1,
1603 i + line * (this->degree - 1) +
1604 4) = -0.5;
1605
1606 for (unsigned int j = 0; j < this->degree - 1;)
1607 {
1608 interpolation_matrix(1,
1609 i + (j + 4) * this->degree - j) =
1610 1.0;
1611 j = j + 2;
1612 }
1613
1614 interpolation_matrix(0, i + 4) = -1.0;
1615 interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1616 i = i + 2;
1617 }
1618
1619 for (unsigned int vertex = 0;
1620 vertex < GeometryInfo<3>::vertices_per_face;
1621 ++vertex)
1622 interpolation_matrix(1, vertex) = 0.25;
1623
1624 int factorial_i = 1;
1625
1626 for (unsigned int i = 2; i <= this->degree; ++i)
1627 {
1628 double tmp = Utilities::pow(0.5, i);
1629 interpolation_matrix(i + 2, i + 2) = tmp;
1630 interpolation_matrix(i + 3 * source_fe.degree - 1,
1631 i + 3 * this->degree - 1) = tmp;
1632 tmp *= 0.5;
1633 interpolation_matrix(i + source_fe.degree + 1, i + 2) =
1634 tmp;
1635 interpolation_matrix(i + source_fe.degree + 1,
1636 i + this->degree + 1) = tmp;
1637 interpolation_matrix(i + 2 * source_fe.degree,
1638 i + 2 * this->degree) = tmp;
1639 interpolation_matrix(i + 2 * source_fe.degree,
1640 i + 3 * this->degree - 1) = tmp;
1641 tmp *= -2.0;
1642
1643 for (unsigned int j = 0; j < this->degree - 1;)
1644 {
1645 interpolation_matrix(i + source_fe.degree + 1,
1646 j + (i + 2) * this->degree + 2 -
1647 i) = tmp;
1648 interpolation_matrix(i + 2 * source_fe.degree,
1649 i + (j + 4) * this->degree - j -
1650 2) = tmp;
1651 j = j + 2;
1652 }
1653
1654 int factorial_k = 1;
1655
1656 for (unsigned int j = 2; j <= this->degree; ++j)
1657 {
1658 interpolation_matrix(i + (j + 2) * source_fe.degree -
1659 j,
1660 i + (j + 2) * this->degree - j) =
1661 Utilities::pow(0.5, i + j);
1662 factorial_k *= j;
1663 int factorial_kl = 1;
1664 int factorial_l = factorial_k;
1665
1666 for (unsigned int k = j + 1; k <= this->degree; ++k)
1667 {
1668 factorial_kl *= k - j;
1669 factorial_l *= k;
1670 interpolation_matrix(
1671 i + (j + 2) * source_fe.degree - j,
1672 i + (k + 2) * this->degree - k) =
1673 Utilities::pow(0.5, i + k) * factorial_l /
1674 (factorial_k * factorial_kl);
1675 }
1676 }
1677
1678 factorial_i *= i;
1679 int factorial_j = factorial_i;
1680 int factorial_ij = 1;
1681
1682 for (unsigned int j = i + 1; j <= this->degree; ++j)
1683 {
1684 factorial_ij *= j - i;
1685 factorial_j *= j;
1686 tmp = Utilities::pow(0.5, j) * factorial_j /
1687 (factorial_i * factorial_ij);
1688 interpolation_matrix(i + 2, j + 2) = tmp;
1689 tmp *= -1.0;
1690
1691 for (unsigned int k = 0; k < this->degree - 1;)
1692 {
1693 interpolation_matrix(i + source_fe.degree + 1,
1694 k + (j + 2) * this->degree +
1695 2 - j) = tmp;
1696 k = k + 2;
1697 }
1698
1699 tmp *= -0.5;
1700 interpolation_matrix(i + source_fe.degree + 1,
1701 j + 2) = tmp;
1702 interpolation_matrix(i + source_fe.degree + 1,
1703 j + this->degree + 1) = tmp;
1704
1705 if (((i + j) & 1) != 0u)
1706 tmp *= -1.0;
1707
1708 interpolation_matrix(i + 2 * source_fe.degree,
1709 j + 2 * this->degree) = tmp;
1710 interpolation_matrix(i + 2 * source_fe.degree,
1711 j + 3 * this->degree - 1) = tmp;
1712 tmp *= 2.0;
1713 interpolation_matrix(i + 3 * source_fe.degree - 1,
1714 j + 3 * this->degree - 1) = tmp;
1715 int factorial_k = 1;
1716
1717 for (unsigned int k = 2; k <= this->degree; ++k)
1718 {
1719 interpolation_matrix(
1720 i + (k + 2) * source_fe.degree - k,
1721 j + (k + 2) * this->degree - k) =
1722 tmp * Utilities::pow(0.5, k);
1723 factorial_k *= k;
1724 int factorial_l = factorial_k;
1725 int factorial_kl = 1;
1726
1727 for (unsigned int l = k + 1; l <= this->degree;
1728 ++l)
1729 {
1730 factorial_kl *= l - k;
1731 factorial_l *= l;
1732 interpolation_matrix(
1733 i + (k + 2) * source_fe.degree - k,
1734 j + (l + 2) * this->degree - l) =
1735 tmp * Utilities::pow(0.5, l) * factorial_l /
1736 (factorial_k * factorial_kl);
1737 }
1738 }
1739
1740 tmp *= -1.0;
1741
1742 for (unsigned int k = 0; k < this->degree - 1;)
1743 {
1744 interpolation_matrix(i + 2 * source_fe.degree,
1745 j + (k + 4) * this->degree -
1746 k - 2) = tmp;
1747 k = k + 2;
1748 }
1749 }
1750 }
1751
1752 break;
1753 }
1754
1755 case 3:
1756 {
1757 for (unsigned int vertex = 0;
1758 vertex < GeometryInfo<3>::vertices_per_face;
1759 ++vertex)
1760 interpolation_matrix(0, vertex) = 0.25;
1761
1762 for (unsigned int i = 0; i < this->degree - 1;)
1763 {
1764 for (unsigned int line = 0;
1765 line < GeometryInfo<3>::lines_per_face;
1766 ++line)
1767 interpolation_matrix(0,
1768 i + line * (this->degree - 1) +
1769 4) = -0.5;
1770
1771 for (unsigned int j = 0; j < this->degree - 1;)
1772 {
1773 interpolation_matrix(0,
1774 i + (j + 4) * this->degree - j) =
1775 1.0;
1776 j = j + 2;
1777 }
1778
1779 interpolation_matrix(1, i + 4) = -1.0;
1780 interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
1781 i = i + 2;
1782 }
1783
1784 interpolation_matrix(1, 0) = 0.5;
1785 interpolation_matrix(1, 1) = 0.5;
1786 interpolation_matrix(2, 2) = 0.5;
1787 interpolation_matrix(2, 3) = 0.5;
1788 interpolation_matrix(3, 3) = 1.0;
1789
1790 int factorial_i = 1;
1791
1792 for (unsigned int i = 2; i <= this->degree; ++i)
1793 {
1794 double tmp = Utilities::pow(0.5, i + 1);
1795 interpolation_matrix(i + 2, i + 2) = tmp;
1796 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1797 interpolation_matrix(i + 2 * source_fe.degree,
1798 i + 2 * this->degree) = tmp;
1799 interpolation_matrix(i + 2 * source_fe.degree,
1800 i + 3 * this->degree - 1) = tmp;
1801 tmp *= -2.0;
1802
1803 for (unsigned int j = 0; j < this->degree - 1;)
1804 {
1805 interpolation_matrix(i + 2,
1806 j + (i + 2) * this->degree + 2 -
1807 i) = tmp;
1808 interpolation_matrix(i + 2 * source_fe.degree,
1809 i + (j + 4) * this->degree - 2) =
1810 tmp;
1811 j = j + 2;
1812 }
1813
1814 tmp *= -1.0;
1815 interpolation_matrix(i + source_fe.degree + 1,
1816 i + this->degree + 1) = tmp;
1817 interpolation_matrix(i + 3 * source_fe.degree - 1,
1818 i + 3 * this->degree - 1) = tmp;
1819 int factorial_k = 1;
1820
1821 for (unsigned int j = 2; j <= this->degree; ++j)
1822 {
1823 interpolation_matrix(i + (j + 2) * source_fe.degree -
1824 j,
1825 i + (j + 2) * this->degree - j) =
1826 Utilities::pow(0.5, i + j);
1827 factorial_k *= j;
1828 int factorial_l = factorial_k;
1829 int factorial_kl = 1;
1830
1831 for (unsigned int k = j + 1; k <= this->degree; ++k)
1832 {
1833 factorial_kl *= k - j;
1834 factorial_l *= k;
1835 interpolation_matrix(
1836 i + (j + 2) * source_fe.degree - j,
1837 i + (k + 2) * this->degree - k) =
1838 Utilities::pow(0.5, i + k) * factorial_l /
1839 (factorial_k * factorial_kl);
1840 }
1841 }
1842
1843 factorial_i *= i;
1844 int factorial_j = factorial_i;
1845 int factorial_ij = 1;
1846
1847 for (unsigned int j = i + 1; j <= this->degree; ++j)
1848 {
1849 factorial_ij *= j - i;
1850 factorial_j *= j;
1851 tmp = Utilities::pow(0.5, j + 1) * factorial_j /
1852 (factorial_i * factorial_ij);
1853 interpolation_matrix(i + 2, j + 2) = tmp;
1854 interpolation_matrix(i + 2, j + this->degree + 1) =
1855 tmp;
1856 interpolation_matrix(i + 2 * source_fe.degree,
1857 j + 2 * this->degree) = tmp;
1858 interpolation_matrix(i + 2 * source_fe.degree,
1859 j + 3 * this->degree - 1) = tmp;
1860 tmp *= 2.0;
1861 interpolation_matrix(i + source_fe.degree + 1,
1862 j + this->degree + 1) = tmp;
1863 interpolation_matrix(i + 3 * source_fe.degree - 1,
1864 j + 3 * this->degree - 1) = tmp;
1865 int factorial_k = 1;
1866
1867 for (unsigned int k = 2; k <= this->degree; ++k)
1868 {
1869 interpolation_matrix(
1870 i + (k + 2) * source_fe.degree - k,
1871 j + (k + 2) * this->degree - k) =
1872 tmp * Utilities::pow(0.5, k);
1873 factorial_k *= k;
1874 int factorial_l = factorial_k;
1875 int factorial_kl = 1;
1876
1877 for (unsigned int l = k + 1; l <= this->degree;
1878 ++l)
1879 {
1880 factorial_kl *= l - k;
1881 factorial_l *= l;
1882 interpolation_matrix(
1883 i + (k + 2) * source_fe.degree - k,
1884 j + (l + 2) * this->degree - l) =
1885 tmp * Utilities::pow(0.5, l) * factorial_l /
1886 (factorial_k * factorial_kl);
1887 }
1888 }
1889
1890 tmp *= -1.0;
1891
1892 for (unsigned int k = 0; k < this->degree - 1;)
1893 {
1894 interpolation_matrix(i + 2,
1895 k + (j + 2) * this->degree +
1896 2 - j) = tmp;
1897 interpolation_matrix(i + 2 * source_fe.degree,
1898 j + (k + 4) * this->degree -
1899 2) = tmp;
1900 k = k + 2;
1901 }
1902 }
1903 }
1904 }
1905 }
1906 }
1907 }
1908}
1909
1910
1911
1912template <int dim>
1913void
1915{
1916 const unsigned int codim = dim - 1;
1917
1918 // TODO: the implementation makes the assumption that all faces have the
1919 // same number of dofs
1920 AssertDimension(this->n_unique_faces(), 1);
1921 const unsigned int face_no = 0;
1922
1923 // number of points: (degree+1)^codim
1924 unsigned int n = this->degree + 1;
1925 for (unsigned int i = 1; i < codim; ++i)
1926 n *= this->degree + 1;
1927
1928 this->generalized_face_support_points[face_no].resize(n);
1929
1930 Point<codim> p;
1931
1932 unsigned int k = 0;
1933 for (unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
1934 for (unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
1935 for (unsigned int ix = 0; ix <= this->degree; ++ix)
1936 {
1937 if (ix == 0)
1938 p[0] = 0.;
1939 else if (ix == 1)
1940 p[0] = 1.;
1941 else
1942 p[0] = .5;
1943 if (codim > 1)
1944 {
1945 if (iy == 0)
1946 p[1] = 0.;
1947 else if (iy == 1)
1948 p[1] = 1.;
1949 else
1950 p[1] = .5;
1951 }
1952 if (codim > 2)
1953 {
1954 if (iz == 0)
1955 p[2] = 0.;
1956 else if (iz == 1)
1957 p[2] = 1.;
1958 else
1959 p[2] = .5;
1960 }
1961 this->generalized_face_support_points[face_no][face_renumber[k++]] =
1962 p;
1963 }
1964}
1965
1966
1967// we use same dpo_vector as FE_Q
1968template <int dim>
1969std::vector<unsigned int>
1971{
1972 std::vector<unsigned int> dpo(dim + 1, 1U);
1973 for (unsigned int i = 1; i < dpo.size(); ++i)
1974 dpo[i] = dpo[i - 1] * (deg - 1);
1975 return dpo;
1976}
1977
1978
1979
1980template <int dim>
1981std::vector<unsigned int>
1983 const FiniteElementData<dim> &fe)
1984{
1985 Assert(fe.n_components() == 1, ExcInternalError());
1986 std::vector<unsigned int> h2l(fe.n_dofs_per_cell());
1987
1988 // polynomial degree
1989 const unsigned int degree = fe.n_dofs_per_line() + 1;
1990 // number of grid points in each
1991 // direction
1992 const unsigned int n = degree + 1;
1993
1994 // the following lines of code are
1995 // somewhat odd, due to the way the
1996 // hierarchic numbering is
1997 // organized. if someone would
1998 // really want to understand these
1999 // lines, you better draw some
2000 // pictures where you indicate the
2001 // indices and orders of vertices,
2002 // lines, etc, along with the
2003 // numbers of the degrees of
2004 // freedom in hierarchical and
2005 // lexicographical order
2006 switch (dim)
2007 {
2008 case 1:
2009 {
2010 for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i)
2011 h2l[i] = i;
2012
2013 break;
2014 }
2015
2016 case 2:
2017 {
2018 // Example: degree=3
2019 //
2020 // hierarchical numbering:
2021 // 2 10 11 3
2022 // 5 14 15 7
2023 // 4 12 13 6
2024 // 0 8 9 1
2025 //
2026 // fe_q_hierarchical numbering:
2027 // 4 6 7 5
2028 // 12 14 15 13
2029 // 8 10 11 9
2030 // 0 2 3 1
2031 unsigned int next_index = 0;
2032 // first the four vertices
2033 h2l[next_index++] = 0;
2034 h2l[next_index++] = 1;
2035 h2l[next_index++] = n;
2036 h2l[next_index++] = n + 1;
2037 // left line
2038 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2039 h2l[next_index++] = (2 + i) * n;
2040 // right line
2041 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2042 h2l[next_index++] = (2 + i) * n + 1;
2043 // bottom line
2044 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2045 h2l[next_index++] = 2 + i;
2046 // top line
2047 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2048 h2l[next_index++] = n + 2 + i;
2049 // inside quad
2050 Assert(fe.n_dofs_per_quad(0 /*only one quad in 2d*/) ==
2051 fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2053 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2054 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2055 h2l[next_index++] = (2 + i) * n + 2 + j;
2056
2057 Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2058
2059 break;
2060 }
2061
2062 case 3:
2063 {
2064 unsigned int next_index = 0;
2065 const unsigned int n2 = n * n;
2066 // first the eight vertices
2067 // bottom face, lexicographic
2068 h2l[next_index++] = 0;
2069 h2l[next_index++] = 1;
2070 h2l[next_index++] = n;
2071 h2l[next_index++] = n + 1;
2072 // top face, lexicographic
2073 h2l[next_index++] = n2;
2074 h2l[next_index++] = n2 + 1;
2075 h2l[next_index++] = n2 + n;
2076 h2l[next_index++] = n2 + n + 1;
2077
2078 // now the lines
2079 // bottom face
2080 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2081 h2l[next_index++] = (2 + i) * n;
2082 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2083 h2l[next_index++] = (2 + i) * n + 1;
2084 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2085 h2l[next_index++] = 2 + i;
2086 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2087 h2l[next_index++] = n + 2 + i;
2088 // top face
2089 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2090 h2l[next_index++] = n2 + (2 + i) * n;
2091 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2092 h2l[next_index++] = n2 + (2 + i) * n + 1;
2093 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2094 h2l[next_index++] = n2 + 2 + i;
2095 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2096 h2l[next_index++] = n2 + n + 2 + i;
2097 // lines in z-direction
2098 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2099 h2l[next_index++] = (2 + i) * n2;
2100 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2101 h2l[next_index++] = (2 + i) * n2 + 1;
2102 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2103 h2l[next_index++] = (2 + i) * n2 + n;
2104 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2105 h2l[next_index++] = (2 + i) * n2 + n + 1;
2106
2107 // TODO: the implementation makes the assumption that all faces have
2108 // the same number of dofs
2110 const unsigned int face_no = 0;
2111 (void)face_no;
2112
2113 // inside quads
2114 Assert(fe.n_dofs_per_quad(face_no) ==
2115 fe.n_dofs_per_line() * fe.n_dofs_per_line(),
2117 // left face
2118 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2119 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2120 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n;
2121 // right face
2122 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2123 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2124 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 1;
2125 // front face
2126 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2127 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2128 h2l[next_index++] = (2 + i) * n2 + 2 + j;
2129 // back face
2130 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2131 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2132 h2l[next_index++] = (2 + i) * n2 + n + 2 + j;
2133 // bottom face
2134 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2135 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2136 h2l[next_index++] = (2 + i) * n + 2 + j;
2137 // top face
2138 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2139 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2140 h2l[next_index++] = n2 + (2 + i) * n + 2 + j;
2141
2142 // inside hex
2143 Assert(fe.n_dofs_per_hex() ==
2144 fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(),
2146 for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
2147 for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
2148 for (unsigned int k = 0; k < fe.n_dofs_per_line(); ++k)
2149 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 2 + k;
2150
2151 Assert(next_index == fe.n_dofs_per_cell(), ExcInternalError());
2152
2153 break;
2154 }
2155
2156 default:
2158 }
2159 return h2l;
2160}
2161
2162
2163template <int dim>
2164std::vector<unsigned int>
2166 const unsigned int degree)
2167{
2168 FiniteElementData<dim - 1> fe_data(
2170 return internal::FE_Q_Hierarchical::invert_numbering(
2172 fe_data));
2173}
2174
2175
2176
2177template <>
2178std::vector<unsigned int>
2180 const unsigned int)
2181{
2182 return std::vector<unsigned int>();
2183}
2184
2185
2186template <>
2187bool
2188FE_Q_Hierarchical<1>::has_support_on_face(const unsigned int shape_index,
2189 const unsigned int face_index) const
2190{
2191 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2193
2194
2195 // in 1d, things are simple. since
2196 // there is only one degree of
2197 // freedom per vertex in this
2198 // class, the first is on vertex 0
2199 // (==face 0 in some sense), the
2200 // second on face 1:
2201 return (((shape_index == 0) && (face_index == 0)) ||
2202 ((shape_index == 1) && (face_index == 1)));
2203}
2204
2205
2206
2207template <int dim>
2208bool
2209FE_Q_Hierarchical<dim>::has_support_on_face(const unsigned int shape_index,
2210 const unsigned int face_index) const
2211{
2212 AssertIndexRange(shape_index, this->n_dofs_per_cell());
2214
2215 // first, special-case interior
2216 // shape functions, since they
2217 // have no support no-where on
2218 // the boundary
2219 if (((dim == 2) && (shape_index >=
2220 this->get_first_quad_index(0 /*only one quad in 2d*/))) ||
2221 ((dim == 3) && (shape_index >= this->get_first_hex_index())))
2222 return false;
2223
2224 // let's see whether this is a
2225 // vertex
2226 if (shape_index < this->get_first_line_index())
2227 {
2228 // for Q elements, there is
2229 // one dof per vertex, so
2230 // shape_index==vertex_number. check
2231 // whether this vertex is
2232 // on the given face.
2233 const unsigned int vertex_no = shape_index;
2236 for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
2237 if (GeometryInfo<dim>::face_to_cell_vertices(face_index, i) ==
2238 vertex_no)
2239 return true;
2240 return false;
2241 }
2242 else if (shape_index < this->get_first_quad_index(0))
2243 // ok, dof is on a line
2244 {
2245 const unsigned int line_index =
2246 (shape_index - this->get_first_line_index()) / this->n_dofs_per_line();
2249
2250 for (unsigned int i = 0; i < GeometryInfo<dim>::lines_per_face; ++i)
2251 if (GeometryInfo<dim>::face_to_cell_lines(face_index, i) == line_index)
2252 return true;
2253 return false;
2254 }
2255 else if (shape_index < this->get_first_hex_index())
2256 // dof is on a quad
2257 {
2258 const unsigned int quad_index =
2259 (shape_index - this->get_first_quad_index(0 /*first quad*/)) /
2260 this->n_dofs_per_quad(face_index);
2261 Assert(static_cast<signed int>(quad_index) <
2262 static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
2264
2265 // in 2d, cell bubble are
2266 // zero on all faces. but
2267 // we have treated this
2268 // case above already
2269 Assert(dim != 2, ExcInternalError());
2270
2271 // in 3d,
2272 // quad_index=face_index
2273 if (dim == 3)
2274 return (quad_index == face_index);
2275 else
2277 }
2278 else
2279 // dof on hex
2280 {
2281 // can only happen in 3d, but
2282 // this case has already been
2283 // covered above
2285 return false;
2286 }
2287
2288 // we should not have gotten here
2290 return false;
2291}
2292
2293
2294
2295template <int dim>
2296std::vector<unsigned int>
2297FE_Q_Hierarchical<dim>::get_embedding_dofs(const unsigned int sub_degree) const
2298{
2299 Assert((sub_degree > 0) && (sub_degree <= this->degree),
2300 ExcIndexRange(sub_degree, 1, this->degree));
2301
2302 if (dim == 1)
2303 {
2304 std::vector<unsigned int> embedding_dofs(sub_degree + 1);
2305 for (unsigned int i = 0; i < (sub_degree + 1); ++i)
2306 embedding_dofs[i] = i;
2307
2308 return embedding_dofs;
2309 }
2310
2311 if (sub_degree == 1)
2312 {
2313 std::vector<unsigned int> embedding_dofs(
2315 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2316 embedding_dofs[i] = i;
2317
2318 return embedding_dofs;
2319 }
2320 else if (sub_degree == this->degree)
2321 {
2322 std::vector<unsigned int> embedding_dofs(this->n_dofs_per_cell());
2323 for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
2324 embedding_dofs[i] = i;
2325
2326 return embedding_dofs;
2327 }
2328
2329 if ((dim == 2) || (dim == 3))
2330 {
2331 std::vector<unsigned int> embedding_dofs(
2332 (dim == 2) ? (sub_degree + 1) * (sub_degree + 1) :
2333 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2334
2335 for (unsigned int i = 0;
2336 i < ((dim == 2) ?
2337 (sub_degree + 1) * (sub_degree + 1) :
2338 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2339 ++i)
2340 {
2341 // vertex
2343 embedding_dofs[i] = i;
2344 // line
2346 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)))
2347 {
2348 const unsigned int j =
2349 (i - GeometryInfo<dim>::vertices_per_cell) % (sub_degree - 1);
2350 const unsigned int line =
2352 (sub_degree - 1);
2353
2354 embedding_dofs[i] = GeometryInfo<dim>::vertices_per_cell +
2355 line * (this->degree - 1) + j;
2356 }
2357 // quad
2359 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) +
2360 GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2361 (sub_degree - 1))
2362 {
2363 const unsigned int j =
2365 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) %
2366 (sub_degree - 1);
2367 const unsigned int k =
2369 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) - j) /
2370 (sub_degree - 1)) %
2371 (sub_degree - 1);
2372 const unsigned int face =
2374 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2375 k * (sub_degree - 1) - j) /
2376 ((sub_degree - 1) * (sub_degree - 1));
2377
2378 embedding_dofs[i] =
2380 GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2381 face * (this->degree - 1) * (this->degree - 1) +
2382 k * (this->degree - 1) + j;
2383 }
2384 // hex
2386 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1)) +
2387 GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2388 (sub_degree - 1) +
2389 GeometryInfo<dim>::hexes_per_cell * (sub_degree - 1) *
2390 (sub_degree - 1) * (sub_degree - 1))
2391 {
2392 const unsigned int j =
2394 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2395 GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2396 (sub_degree - 1)) %
2397 (sub_degree - 1);
2398 const unsigned int k =
2400 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2401 GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2402 (sub_degree - 1) -
2403 j) /
2404 (sub_degree - 1)) %
2405 (sub_degree - 1);
2406 const unsigned int l =
2408 GeometryInfo<dim>::lines_per_cell * (sub_degree - 1) -
2409 GeometryInfo<dim>::quads_per_cell * (sub_degree - 1) *
2410 (sub_degree - 1) -
2411 j - k * (sub_degree - 1)) /
2412 ((sub_degree - 1) * (sub_degree - 1));
2413
2414 embedding_dofs[i] =
2416 GeometryInfo<dim>::lines_per_cell * (this->degree - 1) +
2417 GeometryInfo<dim>::quads_per_cell * (this->degree - 1) *
2418 (this->degree - 1) +
2419 l * (this->degree - 1) * (this->degree - 1) +
2420 k * (this->degree - 1) + j;
2421 }
2422 }
2423
2424 return embedding_dofs;
2425 }
2426 else
2427 {
2429 return std::vector<unsigned int>();
2430 }
2431}
2432
2433
2434
2435template <int dim>
2436std::pair<Table<2, bool>, std::vector<unsigned int>>
2438{
2439 Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
2440 for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
2441 constant_modes(0, i) = true;
2442 for (unsigned int i = GeometryInfo<dim>::vertices_per_cell;
2443 i < this->n_dofs_per_cell();
2444 ++i)
2445 constant_modes(0, i) = false;
2446 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
2447 constant_modes, std::vector<unsigned int>(1, 0));
2448}
2449
2450
2451
2452template <int dim>
2453std::size_t
2459
2460
2461
2462// explicit instantiations
2463#include "fe/fe_q_hierarchical.inst"
2464
2465
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
const unsigned int degree
Definition fe_data.h:452
unsigned int n_dofs_per_cell() const
unsigned int n_dofs_per_line() const
unsigned int n_dofs_per_face(unsigned int face_no=0, unsigned int child=0) const
unsigned int n_components() const
unsigned int n_unique_faces() const
unsigned int n_dofs_per_quad(unsigned int face_no=0) const
unsigned int n_dofs_per_hex() const
virtual std::string get_name() const =0
size_type n() const
size_type m() const
Definition point.h:113
const std::vector< unsigned int > & get_numbering() const
void set_numbering(const std::vector< unsigned int > &renumber)
const std::vector< unsigned int > & get_numbering_inverse() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
#define Assert(cond, exc)
static ::ExceptionBase & ExcImpossibleInDim(int arg1)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & 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)
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()