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