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