42 inline std::vector<unsigned int>
43 invert_numbering(
const std::vector<unsigned int> &in)
45 std::vector<unsigned int> out(in.size());
46 for (
unsigned int i = 0; i < in.size(); ++i)
62 Polynomials::Hierarchical::generate_complete_basis(degree),
72 std::vector<
bool>(1, true)))
73 , face_renumber(face_fe_q_hierarchical_to_hierarchic_numbering(degree))
92 std::vector<FullMatrix<double>> dofs_cell(
96 std::vector<FullMatrix<double>> dofs_subcell(
129 std::ostringstream namebuf;
130 namebuf <<
"FE_Q_Hierarchical<" << dim <<
">(" << this->degree <<
")";
132 return namebuf.str();
138 std::unique_ptr<FiniteElement<dim, dim>>
141 return std_cxx14::make_unique<FE_Q_Hierarchical<dim>>(*this);
172 if (this->dofs_per_cell >= source_fe->dofs_per_cell)
174 const std::vector<unsigned int> dof_map =
175 this->get_embedding_dofs(source_fe->degree);
176 for (
unsigned int j = 0; j < dof_map.size(); j++)
177 matrix[dof_map[j]][j] = 1.;
182 const std::vector<unsigned int> dof_map =
183 source_fe->get_embedding_dofs(this->degree);
184 for (
unsigned int j = 0; j < dof_map.size(); j++)
185 matrix[j][dof_map[j]] = 1.;
193 "Interpolation is supported only between FE_Q_Hierarchical"));
200 const unsigned int child,
206 "Prolongation matrices are only available for isotropic refinement!"));
210 return this->prolongation[refinement_case - 1][child];
223 std::vector<std::pair<unsigned int, unsigned int>>
237 return std::vector<std::pair<unsigned int, unsigned int>>(
238 1, std::make_pair(0
U, 0
U));
246 return std::vector<std::pair<unsigned int, unsigned int>>();
251 return std::vector<std::pair<unsigned int, unsigned int>>();
256 std::vector<std::pair<unsigned int, unsigned int>>
266 const unsigned int this_dpl = this->dofs_per_line;
272 std::vector<std::pair<unsigned int, unsigned int>> res;
273 for (
unsigned int i = 0; i <
std::min(this_dpl, other_dpl); i++)
274 res.emplace_back(i, i);
284 return std::vector<std::pair<unsigned int, unsigned int>>();
289 return std::vector<std::pair<unsigned int, unsigned int>>();
294 std::vector<std::pair<unsigned int, unsigned int>>
304 const unsigned int this_dpq = this->dofs_per_quad;
310 std::vector<std::pair<unsigned int, unsigned int>> res;
311 for (
unsigned int i = 0; i <
std::min(this_dpq, other_dpq); i++)
312 res.emplace_back(i, i);
322 return std::vector<std::pair<unsigned int, unsigned int>>();
327 return std::vector<std::pair<unsigned int, unsigned int>>();
335 const unsigned int codim)
const
343 if (
dynamic_cast<const FE_DGQ<dim> *
>(&fe_other) !=
nullptr)
354 if (this->degree < fe_hierarchical_other->degree)
356 else if (this->degree == fe_hierarchical_other->degree)
364 if (fe_nothing->is_dominating())
389 const unsigned int dofs_1d = 2 * this->dofs_per_vertex + this->dofs_per_line;
438 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell; ++c)
439 for (
unsigned int j = 0; j < dofs_1d; ++j)
440 for (
unsigned int k = 0; k < dofs_1d; ++k)
443 if ((j <= 1) && (k <= 1))
445 if (((c == 0) && (j == 0) && (k == 0)) ||
446 ((c == 1) && (j == 1) && (k == 1)))
447 dofs_cell[c](j, k) = 1.;
449 dofs_cell[c](j, k) = 0.;
451 if (((c == 0) && (j == 1)) || ((c == 1) && (j == 0)))
452 dofs_subcell[c](j, k) = .5;
453 else if (((c == 0) && (k == 0)) || ((c == 1) && (k == 1)))
454 dofs_subcell[c](j, k) = 1.;
456 dofs_subcell[c](j, k) = 0.;
459 else if ((j <= 1) && (k >= 2))
461 if (((c == 0) && (j == 1) && ((k % 2) == 0)) ||
462 ((c == 1) && (j == 0) && ((k % 2) == 0)))
463 dofs_subcell[c](j, k) = -1.;
466 else if ((j >= 2) && (k >= 2) && (j <= k))
469 for (
unsigned int i = 1; i <= j; ++i)
471 (
static_cast<double>(k - i + 1)) / (
static_cast<double>(i));
475 dofs_subcell[c](j, k) =
477 std::pow(.5,
static_cast<double>(k)) * factor :
478 -std::pow(.5,
static_cast<double>(k)) * factor;
480 std::pow(2.,
static_cast<double>(j)) * factor;
484 dofs_subcell[c](j, k) =
485 std::pow(.5,
static_cast<double>(k)) * factor;
488 std::pow(2.,
static_cast<double>(j)) * factor :
489 -std::pow(2.,
static_cast<double>(j)) * factor;
502 const unsigned int dofs_1d = 2 * this->dofs_per_vertex + this->dofs_per_line;
504 this->interface_constraints.TableBase<2, double>::reinit(
505 this->interface_constraints_size());
518 for (
unsigned int i = 0; i < dofs_1d; ++i)
519 this->interface_constraints(0, i) = dofs_subcell[0](1, i);
521 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
523 for (
unsigned int i = 0; i < dofs_1d; ++i)
524 for (
unsigned int j = 2; j < dofs_1d; ++j)
525 this->interface_constraints(1 + c * (this->degree - 1) + j - 2,
526 i) = dofs_subcell[c](j, i);
532 for (
unsigned int i = 0; i < dofs_1d * dofs_1d; i++)
535 this->interface_constraints(0, face_renumber[i]) =
536 dofs_subcell[0](1, i % dofs_1d) *
537 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
540 this->interface_constraints(1, face_renumber[i]) =
541 dofs_subcell[0](0, i % dofs_1d) *
542 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
543 this->interface_constraints(2, face_renumber[i]) =
544 dofs_subcell[1](1, i % dofs_1d) *
545 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
546 this->interface_constraints(3, face_renumber[i]) =
547 dofs_subcell[0](1, i % dofs_1d) *
548 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
549 this->interface_constraints(4, face_renumber[i]) =
550 dofs_subcell[1](0, i % dofs_1d) *
551 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
554 for (
unsigned int j = 0; j < (this->degree - 1); j++)
556 this->interface_constraints(5 + j, face_renumber[i]) =
557 dofs_subcell[0](1, i % dofs_1d) *
558 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
559 this->interface_constraints(5 + (this->degree - 1) + j,
561 dofs_subcell[0](1, i % dofs_1d) *
562 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
563 this->interface_constraints(5 + 2 * (this->degree - 1) + j,
565 dofs_subcell[0](2 + j, i % dofs_1d) *
566 dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
567 this->interface_constraints(5 + 3 * (this->degree - 1) + j,
569 dofs_subcell[1](2 + j, i % dofs_1d) *
570 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
574 for (
unsigned int j = 0; j < (this->degree - 1); j++)
577 this->interface_constraints(5 + 4 * (this->degree - 1) + j,
579 dofs_subcell[0](0, i % dofs_1d) *
580 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
581 this->interface_constraints(5 + 4 * (this->degree - 1) +
582 (this->degree - 1) + j,
584 dofs_subcell[0](0, i % dofs_1d) *
585 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
587 this->interface_constraints(5 + 4 * (this->degree - 1) +
588 2 * (this->degree - 1) + j,
590 dofs_subcell[1](1, i % dofs_1d) *
591 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
592 this->interface_constraints(5 + 4 * (this->degree - 1) +
593 3 * (this->degree - 1) + j,
595 dofs_subcell[1](1, i % dofs_1d) *
596 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
598 this->interface_constraints(5 + 4 * (this->degree - 1) +
599 4 * (this->degree - 1) + j,
601 dofs_subcell[0](2 + j, i % dofs_1d) *
602 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
603 this->interface_constraints(5 + 4 * (this->degree - 1) +
604 5 * (this->degree - 1) + j,
606 dofs_subcell[1](2 + j, i % dofs_1d) *
607 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
609 this->interface_constraints(5 + 4 * (this->degree - 1) +
610 6 * (this->degree - 1) + j,
612 dofs_subcell[0](2 + j, i % dofs_1d) *
613 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
614 this->interface_constraints(5 + 4 * (this->degree - 1) +
615 7 * (this->degree - 1) + j,
617 dofs_subcell[1](2 + j, i % dofs_1d) *
618 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
622 for (
unsigned int j = 0; j < (this->degree - 1); j++)
623 for (
unsigned int k = 0; k < (this->degree - 1); k++)
626 this->interface_constraints(5 + 12 * (this->degree - 1) +
627 j + k * (this->degree - 1),
629 dofs_subcell[0](2 + j, i % dofs_1d) *
630 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
632 this->interface_constraints(5 + 12 * (this->degree - 1) +
633 j + k * (this->degree - 1) +
637 dofs_subcell[1](2 + j, i % dofs_1d) *
638 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
640 this->interface_constraints(5 + 12 * (this->degree - 1) +
641 j + k * (this->degree - 1) +
642 2 * (this->degree - 1) *
645 dofs_subcell[0](2 + j, i % dofs_1d) *
646 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
648 this->interface_constraints(5 + 12 * (this->degree - 1) +
649 j + k * (this->degree - 1) +
650 3 * (this->degree - 1) *
653 dofs_subcell[1](2 + j, i % dofs_1d) *
654 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
675 const unsigned int dofs_1d = 2 * this->dofs_per_vertex + this->dofs_per_line;
676 const std::vector<unsigned int> &renumber = this->poly_space.get_numbering();
678 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell; ++c)
680 this->prolongation[iso][c].reinit(this->dofs_per_cell,
681 this->dofs_per_cell);
682 this->restriction[iso][c].reinit(this->dofs_per_cell,
683 this->dofs_per_cell);
690 for (
unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
693 this->prolongation[iso][c].fill(dofs_subcell[c]);
694 this->restriction[iso][c].fill(dofs_cell[c]);
709 for (
unsigned int j = 0; j < this->dofs_per_cell; ++j)
710 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
715 for (
unsigned int c = 0;
716 c < GeometryInfo<2>::max_children_per_cell;
720 const unsigned int c0 = c % 2;
722 const unsigned int c1 = c / 2;
724 this->prolongation[iso][c](j, i) =
725 dofs_subcell[c0](renumber[j] % dofs_1d,
726 renumber[i] % dofs_1d) *
727 dofs_subcell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
729 (renumber[i] - (renumber[i] % dofs_1d)) /
732 this->restriction[iso][c](j, i) =
733 dofs_cell[c0](renumber[j] % dofs_1d,
734 renumber[i] % dofs_1d) *
735 dofs_cell[c1]((renumber[j] - (renumber[j] % dofs_1d)) /
737 (renumber[i] - (renumber[i] % dofs_1d)) /
745 for (
unsigned int c = 0;
746 c < GeometryInfo<3>::max_children_per_cell;
750 const unsigned int c0 = c % 2;
752 const unsigned int c1 = (c % 4) / 2;
754 const unsigned int c2 = c / 4;
756 this->prolongation[iso][c](j, i) =
757 dofs_subcell[c0](renumber[j] % dofs_1d,
758 renumber[i] % dofs_1d) *
760 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
762 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
765 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
766 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
769 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
770 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
774 this->restriction[iso][c](j, i) =
775 dofs_cell[c0](renumber[j] % dofs_1d,
776 renumber[i] % dofs_1d) *
778 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
780 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
783 ((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d -
784 (((renumber[j] - (renumber[j] % dofs_1d)) / dofs_1d) %
787 ((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d -
788 (((renumber[i] - (renumber[i] % dofs_1d)) / dofs_1d) %
807 unsigned int n = this->degree + 1;
808 for (
unsigned int i = 1; i < dim; ++i)
809 n *= this->degree + 1;
811 this->generalized_support_points.resize(n);
813 const std::vector<unsigned int> &index_map_inverse =
814 this->poly_space.get_numbering_inverse();
845 for (
unsigned int iz = 0; iz <= ((dim > 2) ? this->degree : 0); ++iz)
846 for (
unsigned int iy = 0; iy <= ((dim > 1) ? this->degree : 0); ++iy)
847 for (
unsigned int ix = 0; ix <= this->degree; ++ix)
873 this->generalized_support_points[index_map_inverse[k++]] = p;
920 (
dynamic_cast<const FEQHierarchical *
>(&x_source_fe) !=
924 Assert(interpolation_matrix.
n() == this->dofs_per_face,
949 interpolation_matrix = 0;
961 for (
unsigned int i = 0; i < this->dofs_per_face; ++i)
962 interpolation_matrix(i, i) = 1;
969 for (
unsigned int i = 0; i < GeometryInfo<3>::vertices_per_face; ++i)
970 interpolation_matrix(i, i) = 1;
972 for (
unsigned int i = 0; i < this->degree - 1; ++i)
974 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face; ++j)
975 interpolation_matrix(i + j * (x_source_fe.
degree - 1) +
977 i + j * (this->degree - 1) +
980 for (
unsigned int j = 0; j < this->degree - 1; ++j)
982 (x_source_fe.
degree - 1) +
999 const unsigned int subface,
1007 (
dynamic_cast<const FEQHierarchical *
>(&x_source_fe) !=
1011 Assert(interpolation_matrix.
n() == this->dofs_per_face,
1044 interpolation_matrix(0, 0) = 1.0;
1045 interpolation_matrix(1, 0) = 0.5;
1046 interpolation_matrix(1, 1) = 0.5;
1048 for (
unsigned int dof = 2; dof < this->dofs_per_face;)
1050 interpolation_matrix(1, dof) = -1.0;
1054 int factorial_i = 1;
1056 for (
unsigned int i = 2; i < this->dofs_per_face; ++i)
1058 interpolation_matrix(i, i) = std::pow(0.5, i);
1060 int factorial_j = factorial_i;
1061 int factorial_ij = 1;
1063 for (
unsigned int j = i + 1; j < this->dofs_per_face; ++j)
1065 factorial_ij *= j - i;
1069 interpolation_matrix(i, j) =
1070 -1.0 * std::pow(0.5, j) * factorial_j /
1071 (factorial_i * factorial_ij);
1074 interpolation_matrix(i, j) =
1075 std::pow(0.5, j) * factorial_j /
1076 (factorial_i * factorial_ij);
1085 interpolation_matrix(0, 0) = 0.5;
1086 interpolation_matrix(0, 1) = 0.5;
1088 for (
unsigned int dof = 2; dof < this->dofs_per_face;)
1090 interpolation_matrix(0, dof) = -1.0;
1094 interpolation_matrix(1, 1) = 1.0;
1096 int factorial_i = 1;
1098 for (
unsigned int i = 2; i < this->dofs_per_face; ++i)
1100 interpolation_matrix(i, i) = std::pow(0.5, i);
1102 int factorial_j = factorial_i;
1103 int factorial_ij = 1;
1105 for (
unsigned int j = i + 1; j < this->dofs_per_face; ++j)
1107 factorial_ij *= j - i;
1109 interpolation_matrix(i, j) =
1110 std::pow(0.5, j) * factorial_j /
1111 (factorial_i * factorial_ij);
1126 interpolation_matrix(0, 0) = 1.0;
1127 interpolation_matrix(1, 0) = 0.5;
1128 interpolation_matrix(1, 1) = 0.5;
1129 interpolation_matrix(2, 0) = 0.5;
1130 interpolation_matrix(2, 2) = 0.5;
1132 for (
unsigned int i = 0; i < this->degree - 1;)
1134 for (
unsigned int line = 0;
1135 line < GeometryInfo<3>::lines_per_face;
1137 interpolation_matrix(3,
1138 i + line * (this->degree - 1) +
1141 for (
unsigned int j = 0; j < this->degree - 1;)
1143 interpolation_matrix(3,
1144 i + (j + 4) * this->degree - j) =
1149 interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1151 interpolation_matrix(2, i + 4) = -1.0;
1155 for (
unsigned int vertex = 0;
1156 vertex < GeometryInfo<3>::vertices_per_face;
1158 interpolation_matrix(3, vertex) = 0.25;
1160 int factorial_i = 1;
1162 for (
unsigned int i = 2; i <= this->degree; ++i)
1164 double tmp = std::pow(0.5, i);
1165 interpolation_matrix(i + 2, i + 2) = tmp;
1166 interpolation_matrix(i + 2 * source_fe.
degree,
1167 i + 2 * this->degree) = tmp;
1169 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
1171 interpolation_matrix(i + source_fe.
degree + 1,
1172 i + this->degree + 1) = tmp;
1173 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1174 i + 2 * this->degree) = tmp;
1175 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1176 i + 3 * this->degree - 1) = tmp;
1179 for (
unsigned int j = 0; j < this->degree - 1;)
1181 interpolation_matrix(i + source_fe.
degree + 1,
1182 (i + 2) * this->degree + j + 2 -
1184 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1185 i + (j + 4) * this->degree - j -
1190 int factorial_k = 1;
1192 for (
unsigned int j = 2; j <= this->degree; ++j)
1194 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1196 i + (j + 2) * this->degree - j) =
1197 std::pow(0.5, i + j);
1199 int factorial_kl = 1;
1200 int factorial_l = factorial_k;
1202 for (
unsigned int k = j + 1; k < this->degree; ++k)
1204 factorial_kl *= k - j;
1208 interpolation_matrix(
1209 i + (j + 2) * source_fe.
degree - j,
1210 i + (k + 2) * this->degree - k) =
1211 -1.0 * std::pow(0.5, i + k) * factorial_l /
1212 (factorial_k * factorial_kl);
1215 interpolation_matrix(
1216 i + (j + 2) * source_fe.
degree - j,
1217 i + (k + 2) * this->degree - k) =
1218 std::pow(0.5, i + k) * factorial_l /
1219 (factorial_k * factorial_kl);
1224 int factorial_j = factorial_i;
1225 int factorial_ij = 1;
1227 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1229 factorial_ij *= j - i;
1234 tmp = -1.0 * std::pow(0.5, j) * factorial_j /
1235 (factorial_i * factorial_ij);
1236 interpolation_matrix(i + 2, j + 2) = tmp;
1237 interpolation_matrix(i + 2 * source_fe.
degree,
1238 j + 2 * this->degree) = tmp;
1241 for (
unsigned int k = 2; k <= this->degree; ++k)
1243 interpolation_matrix(
1244 i + (k + 2) * source_fe.
degree - k,
1245 j + (k + 2) * this->degree - k) =
1246 tmp * std::pow(0.5, k);
1248 int factorial_l = factorial_k;
1249 int factorial_kl = 1;
1251 for (
unsigned int l = k + 1;
1255 factorial_kl *=
l - k;
1259 interpolation_matrix(
1260 i + (k + 2) * source_fe.
degree - k,
1261 j + (
l + 2) * this->degree -
l) =
1262 -1.0 * tmp * std::pow(0.5,
l) *
1264 (factorial_k * factorial_kl);
1267 interpolation_matrix(
1268 i + (k + 2) * source_fe.
degree - k,
1269 j + (
l + 2) * this->degree -
l) =
1270 tmp * std::pow(0.5,
l) * factorial_l /
1271 (factorial_k * factorial_kl);
1276 interpolation_matrix(i + source_fe.
degree + 1,
1278 interpolation_matrix(i + source_fe.
degree + 1,
1279 j + this->degree + 1) = tmp;
1280 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1281 j + 2 * this->degree) = tmp;
1282 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1283 j + 3 * this->degree - 1) =
1287 for (
unsigned int k = 0; k < this->degree - 1;)
1289 interpolation_matrix(i + source_fe.
degree + 1,
1290 (j + 2) * this->degree +
1292 interpolation_matrix(
1293 i + 3 * source_fe.
degree - 1,
1294 j + (k + 4) * this->degree - k - 2) = tmp;
1300 tmp = std::pow(0.5, j) * factorial_j /
1301 (factorial_i * factorial_ij);
1302 interpolation_matrix(i + 2, j + 2) = tmp;
1303 interpolation_matrix(i + 2 * source_fe.
degree,
1304 j + 2 * this->degree) = tmp;
1307 for (
unsigned int k = 2; k <= this->degree; ++k)
1309 interpolation_matrix(
1310 i + (k + 2) * source_fe.
degree - k,
1311 j + (k + 2) * this->degree - k) =
1312 tmp * std::pow(0.5, k);
1314 int factorial_l = factorial_k;
1315 int factorial_kl = 1;
1317 for (
unsigned int l = k + 1;
1321 factorial_kl *=
l - k;
1325 interpolation_matrix(
1326 i + (k + 2) * source_fe.
degree - k,
1327 j + (
l + 2) * this->degree -
l) =
1328 -1.0 * tmp * std::pow(0.5,
l) *
1330 (factorial_k * factorial_kl);
1333 interpolation_matrix(
1334 i + (k + 2) * source_fe.
degree - k,
1335 j + (
l + 2) * this->degree -
l) =
1336 tmp * std::pow(0.5,
l) * factorial_l /
1337 (factorial_k * factorial_kl);
1342 interpolation_matrix(i + source_fe.
degree + 1,
1344 interpolation_matrix(i + source_fe.
degree + 1,
1345 j + this->degree + 1) = tmp;
1346 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1347 j + 2 * this->degree) = tmp;
1348 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1349 j + 3 * this->degree - 1) =
1353 for (
unsigned int k = 0; k < this->degree - 1;)
1355 interpolation_matrix(i + source_fe.
degree + 1,
1356 (j + 2) * this->degree +
1358 interpolation_matrix(
1359 i + 3 * source_fe.
degree - 1,
1360 j + (k + 4) * this->degree - k - 2) = tmp;
1372 interpolation_matrix(0, 0) = 0.5;
1373 interpolation_matrix(0, 1) = 0.5;
1374 interpolation_matrix(1, 1) = 1.0;
1375 interpolation_matrix(3, 1) = 0.5;
1376 interpolation_matrix(3, 3) = 0.5;
1378 for (
unsigned int i = 0; i < this->degree - 1;)
1380 for (
unsigned int line = 0;
1381 line < GeometryInfo<3>::lines_per_face;
1383 interpolation_matrix(2,
1384 i + line * (this->degree - 1) +
1387 for (
unsigned int j = 0; j < this->degree - 1;)
1389 interpolation_matrix(2,
1390 i + (j + 4) * this->degree - j) =
1395 interpolation_matrix(0, i + 2 * (this->degree + 1)) =
1397 interpolation_matrix(3, i + this->degree + 3) = -1.0;
1401 for (
unsigned int vertex = 0;
1402 vertex < GeometryInfo<3>::vertices_per_face;
1404 interpolation_matrix(2, vertex) = 0.25;
1406 int factorial_i = 1;
1408 for (
unsigned int i = 2; i <= this->degree; ++i)
1410 double tmp = std::pow(0.5, i + 1);
1411 interpolation_matrix(i + 2, i + 2) = tmp;
1412 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1413 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1414 i + 2 * this->degree) = tmp;
1415 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1416 i + 3 * this->degree - 1) = tmp;
1419 for (
unsigned int j = 0; j < this->degree - 1;)
1421 interpolation_matrix(i + 2,
1422 j + (i + 2) * this->degree + 2 -
1424 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1425 i + (j + 4) * this->degree - j -
1431 interpolation_matrix(i + source_fe.
degree + 1,
1432 i + this->degree + 1) = tmp;
1433 interpolation_matrix(i + 2 * source_fe.
degree,
1434 i + 2 * this->degree) = tmp;
1436 int factorial_j = factorial_i;
1437 int factorial_ij = 1;
1439 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1441 factorial_ij *= j - i;
1443 tmp = std::pow(0.5, j) * factorial_j /
1444 (factorial_i * factorial_ij);
1445 interpolation_matrix(i + 2 * source_fe.
degree,
1446 j + 2 * this->degree) = tmp;
1447 int factorial_k = 1;
1449 for (
unsigned int k = 2; k <= this->degree; ++k)
1451 interpolation_matrix(
1452 i + (k + 2) * source_fe.
degree - k,
1453 j + (k + 2) * this->degree - k) =
1454 tmp * std::pow(0.5, k);
1456 int factorial_l = factorial_k;
1457 int factorial_kl = 1;
1459 for (
unsigned int l = k + 1;
l <= this->degree;
1462 factorial_kl *=
l - k;
1466 interpolation_matrix(
1467 i + (k + 2) * source_fe.
degree - k,
1468 j + (
l + 2) * this->degree -
l) =
1469 -1.0 * tmp * std::pow(0.5,
l) *
1471 (factorial_k * factorial_kl);
1474 interpolation_matrix(
1475 i + (k + 2) * source_fe.
degree - k,
1476 j + (
l + 2) * this->degree -
l) =
1477 tmp * std::pow(0.5,
l) * factorial_l /
1478 (factorial_k * factorial_kl);
1484 for (
unsigned int k = 0; k < this->degree - 1;)
1486 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1487 j + (k + 4) * this->degree -
1493 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1494 j + 2 * this->degree) = tmp;
1495 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1496 j + 3 * this->degree - 1) = tmp;
1501 interpolation_matrix(i + 2, j + 2) = tmp;
1502 interpolation_matrix(i + 2, j + this->degree + 1) =
1504 interpolation_matrix(i + source_fe.
degree + 1,
1505 j + this->degree + 1) =
1509 for (
unsigned int k = 0; k < this->degree - 1;)
1511 interpolation_matrix(i + 2,
1512 k + (j + 2) * this->degree +
1518 int factorial_k = 1;
1520 for (
unsigned int j = 2; j <= this->degree; ++j)
1522 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1524 i + (j + 2) * this->degree - j) =
1525 std::pow(0.5, i + j);
1527 int factorial_l = factorial_k;
1528 int factorial_kl = 1;
1530 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1532 factorial_kl *= k - j;
1536 interpolation_matrix(
1537 i + (j + 2) * source_fe.
degree - j,
1538 i + (k + 2) * this->degree - k) =
1539 -1.0 * std::pow(0.5, i + k) * factorial_l /
1540 (factorial_k * factorial_kl);
1543 interpolation_matrix(
1544 i + (j + 2) * source_fe.
degree - j,
1545 i + (k + 2) * this->degree - k) =
1546 std::pow(0.5, i + k) * factorial_l /
1547 (factorial_k * factorial_kl);
1557 interpolation_matrix(0, 0) = 0.5;
1558 interpolation_matrix(0, 2) = 0.5;
1559 interpolation_matrix(2, 2) = 1.0;
1560 interpolation_matrix(3, 2) = 0.5;
1561 interpolation_matrix(3, 3) = 0.5;
1563 for (
unsigned int i = 0; i < this->degree - 1;)
1565 for (
unsigned int line = 0;
1566 line < GeometryInfo<3>::lines_per_face;
1568 interpolation_matrix(1,
1569 i + line * (this->degree - 1) +
1572 for (
unsigned int j = 0; j < this->degree - 1;)
1574 interpolation_matrix(1,
1575 i + (j + 4) * this->degree - j) =
1580 interpolation_matrix(0, i + 4) = -1.0;
1581 interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1585 for (
unsigned int vertex = 0;
1586 vertex < GeometryInfo<3>::vertices_per_face;
1588 interpolation_matrix(1, vertex) = 0.25;
1590 int factorial_i = 1;
1592 for (
unsigned int i = 2; i <= this->degree; ++i)
1594 double tmp = std::pow(0.5, i);
1595 interpolation_matrix(i + 2, i + 2) = tmp;
1596 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1597 i + 3 * this->degree - 1) = tmp;
1599 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
1601 interpolation_matrix(i + source_fe.
degree + 1,
1602 i + this->degree + 1) = tmp;
1603 interpolation_matrix(i + 2 * source_fe.
degree,
1604 i + 2 * this->degree) = tmp;
1605 interpolation_matrix(i + 2 * source_fe.
degree,
1606 i + 3 * this->degree - 1) = tmp;
1609 for (
unsigned int j = 0; j < this->degree - 1;)
1611 interpolation_matrix(i + source_fe.
degree + 1,
1612 j + (i + 2) * this->degree + 2 -
1614 interpolation_matrix(i + 2 * source_fe.
degree,
1615 i + (j + 4) * this->degree - j -
1620 int factorial_k = 1;
1622 for (
unsigned int j = 2; j <= this->degree; ++j)
1624 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1626 i + (j + 2) * this->degree - j) =
1627 std::pow(0.5, i + j);
1629 int factorial_kl = 1;
1630 int factorial_l = factorial_k;
1632 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1634 factorial_kl *= k - j;
1636 interpolation_matrix(
1637 i + (j + 2) * source_fe.
degree - j,
1638 i + (k + 2) * this->degree - k) =
1639 std::pow(0.5, i + k) * factorial_l /
1640 (factorial_k * factorial_kl);
1645 int factorial_j = factorial_i;
1646 int factorial_ij = 1;
1648 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1650 factorial_ij *= j - i;
1652 tmp = std::pow(0.5, j) * factorial_j /
1653 (factorial_i * factorial_ij);
1654 interpolation_matrix(i + 2, j + 2) = tmp;
1657 for (
unsigned int k = 0; k < this->degree - 1;)
1659 interpolation_matrix(i + source_fe.
degree + 1,
1660 k + (j + 2) * this->degree +
1666 interpolation_matrix(i + source_fe.
degree + 1,
1668 interpolation_matrix(i + source_fe.
degree + 1,
1669 j + this->degree + 1) = tmp;
1674 interpolation_matrix(i + 2 * source_fe.
degree,
1675 j + 2 * this->degree) = tmp;
1676 interpolation_matrix(i + 2 * source_fe.
degree,
1677 j + 3 * this->degree - 1) = tmp;
1679 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1680 j + 3 * this->degree - 1) = tmp;
1681 int factorial_k = 1;
1683 for (
unsigned int k = 2; k <= this->degree; ++k)
1685 interpolation_matrix(
1686 i + (k + 2) * source_fe.
degree - k,
1687 j + (k + 2) * this->degree - k) =
1688 tmp * std::pow(0.5, k);
1690 int factorial_l = factorial_k;
1691 int factorial_kl = 1;
1693 for (
unsigned int l = k + 1;
l <= this->degree;
1696 factorial_kl *=
l - k;
1698 interpolation_matrix(
1699 i + (k + 2) * source_fe.
degree - k,
1700 j + (
l + 2) * this->degree -
l) =
1701 tmp * std::pow(0.5,
l) * factorial_l /
1702 (factorial_k * factorial_kl);
1708 for (
unsigned int k = 0; k < this->degree - 1;)
1710 interpolation_matrix(i + 2 * source_fe.
degree,
1711 j + (k + 4) * this->degree -
1723 for (
unsigned int vertex = 0;
1724 vertex < GeometryInfo<3>::vertices_per_face;
1726 interpolation_matrix(0, vertex) = 0.25;
1728 for (
unsigned int i = 0; i < this->degree - 1;)
1730 for (
unsigned int line = 0;
1731 line < GeometryInfo<3>::lines_per_face;
1733 interpolation_matrix(0,
1734 i + line * (this->degree - 1) +
1737 for (
unsigned int j = 0; j < this->degree - 1;)
1739 interpolation_matrix(0,
1740 i + (j + 4) * this->degree - j) =
1745 interpolation_matrix(1, i + 4) = -1.0;
1746 interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
1750 interpolation_matrix(1, 0) = 0.5;
1751 interpolation_matrix(1, 1) = 0.5;
1752 interpolation_matrix(2, 2) = 0.5;
1753 interpolation_matrix(2, 3) = 0.5;
1754 interpolation_matrix(3, 3) = 1.0;
1756 int factorial_i = 1;
1758 for (
unsigned int i = 2; i <= this->degree; ++i)
1760 double tmp = std::pow(0.5, i + 1);
1761 interpolation_matrix(i + 2, i + 2) = tmp;
1762 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1763 interpolation_matrix(i + 2 * source_fe.
degree,
1764 i + 2 * this->degree) = tmp;
1765 interpolation_matrix(i + 2 * source_fe.
degree,
1766 i + 3 * this->degree - 1) = tmp;
1769 for (
unsigned int j = 0; j < this->degree - 1;)
1771 interpolation_matrix(i + 2,
1772 j + (i + 2) * this->degree + 2 -
1774 interpolation_matrix(i + 2 * source_fe.
degree,
1775 i + (j + 4) * this->degree - 2) =
1781 interpolation_matrix(i + source_fe.
degree + 1,
1782 i + this->degree + 1) = tmp;
1783 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1784 i + 3 * this->degree - 1) = tmp;
1785 int factorial_k = 1;
1787 for (
unsigned int j = 2; j <= this->degree; ++j)
1789 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1791 i + (j + 2) * this->degree - j) =
1792 std::pow(0.5, i + j);
1794 int factorial_l = factorial_k;
1795 int factorial_kl = 1;
1797 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1799 factorial_kl *= k - j;
1801 interpolation_matrix(
1802 i + (j + 2) * source_fe.
degree - j,
1803 i + (k + 2) * this->degree - k) =
1804 std::pow(0.5, i + k) * factorial_l /
1805 (factorial_k * factorial_kl);
1810 int factorial_j = factorial_i;
1811 int factorial_ij = 1;
1813 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1815 factorial_ij *= j - i;
1817 tmp = std::pow(0.5, j + 1) * factorial_j /
1818 (factorial_i * factorial_ij);
1819 interpolation_matrix(i + 2, j + 2) = tmp;
1820 interpolation_matrix(i + 2, j + this->degree + 1) =
1822 interpolation_matrix(i + 2 * source_fe.
degree,
1823 j + 2 * this->degree) = tmp;
1824 interpolation_matrix(i + 2 * source_fe.
degree,
1825 j + 3 * this->degree - 1) = tmp;
1827 interpolation_matrix(i + source_fe.
degree + 1,
1828 j + this->degree + 1) = tmp;
1829 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1830 j + 3 * this->degree - 1) = tmp;
1831 int factorial_k = 1;
1833 for (
unsigned int k = 2; k <= this->degree; ++k)
1835 interpolation_matrix(
1836 i + (k + 2) * source_fe.
degree - k,
1837 j + (k + 2) * this->degree - k) =
1838 tmp * std::pow(0.5, k);
1840 int factorial_l = factorial_k;
1841 int factorial_kl = 1;
1843 for (
unsigned int l = k + 1;
l <= this->degree;
1846 factorial_kl *=
l - k;
1848 interpolation_matrix(
1849 i + (k + 2) * source_fe.
degree - k,
1850 j + (
l + 2) * this->degree -
l) =
1851 tmp * std::pow(0.5,
l) * factorial_l /
1852 (factorial_k * factorial_kl);
1858 for (
unsigned int k = 0; k < this->degree - 1;)
1860 interpolation_matrix(i + 2,
1861 k + (j + 2) * this->degree +
1863 interpolation_matrix(i + 2 * source_fe.
degree,
1864 j + (k + 4) * this->degree -
1882 const unsigned int codim = dim - 1;
1885 unsigned int n = this->degree + 1;
1886 for (
unsigned int i = 1; i < codim; ++i)
1887 n *= this->degree + 1;
1889 this->generalized_face_support_points.resize(n);
1894 for (
unsigned int iz = 0; iz <= ((codim > 2) ? this->degree : 0); ++iz)
1895 for (
unsigned int iy = 0; iy <= ((codim > 1) ? this->degree : 0); ++iy)
1896 for (
unsigned int ix = 0; ix <= this->degree; ++ix)
1922 this->generalized_face_support_points[face_renumber[k++]] = p;
1929 std::vector<unsigned int>
1932 std::vector<unsigned int> dpo(dim + 1, 1
U);
1933 for (
unsigned int i = 1; i < dpo.size(); ++i)
1934 dpo[i] = dpo[i - 1] * (deg - 1);
1941 std::vector<unsigned int>
1952 const unsigned int n = degree + 1;
1991 unsigned int next_index = 0;
1993 h2l[next_index++] = 0;
1994 h2l[next_index++] = 1;
1995 h2l[next_index++] = n;
1996 h2l[next_index++] = n + 1;
1999 h2l[next_index++] = (2 + i) * n;
2002 h2l[next_index++] = (2 + i) * n + 1;
2005 h2l[next_index++] = 2 + i;
2008 h2l[next_index++] = n + 2 + i;
2014 h2l[next_index++] = (2 + i) * n + 2 + j;
2023 unsigned int next_index = 0;
2024 const unsigned int n2 = n * n;
2027 h2l[next_index++] = 0;
2028 h2l[next_index++] = 1;
2029 h2l[next_index++] = n;
2030 h2l[next_index++] = n + 1;
2032 h2l[next_index++] = n2;
2033 h2l[next_index++] = n2 + 1;
2034 h2l[next_index++] = n2 + n;
2035 h2l[next_index++] = n2 + n + 1;
2040 h2l[next_index++] = (2 + i) * n;
2042 h2l[next_index++] = (2 + i) * n + 1;
2044 h2l[next_index++] = 2 + i;
2046 h2l[next_index++] = n + 2 + i;
2049 h2l[next_index++] = n2 + (2 + i) * n;
2051 h2l[next_index++] = n2 + (2 + i) * n + 1;
2053 h2l[next_index++] = n2 + 2 + i;
2055 h2l[next_index++] = n2 + n + 2 + i;
2058 h2l[next_index++] = (2 + i) * n2;
2060 h2l[next_index++] = (2 + i) * n2 + 1;
2062 h2l[next_index++] = (2 + i) * n2 + n;
2064 h2l[next_index++] = (2 + i) * n2 + n + 1;
2072 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n;
2076 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 1;
2080 h2l[next_index++] = (2 + i) * n2 + 2 + j;
2084 h2l[next_index++] = (2 + i) * n2 + n + 2 + j;
2088 h2l[next_index++] = (2 + i) * n + 2 + j;
2092 h2l[next_index++] = n2 + (2 + i) * n + 2 + j;
2100 h2l[next_index++] = (2 + i) * n2 + (2 + j) * n + 2 + k;
2115 std::vector<unsigned int>
2117 const unsigned int degree)
2121 return internal::FE_Q_Hierarchical::invert_numbering(
2129 std::vector<unsigned int>
2133 return std::vector<unsigned int>();
2140 const unsigned int face_index)
const
2152 return (((shape_index == 0) && (face_index == 0)) ||
2153 ((shape_index == 1) && (face_index == 1)));
2161 const unsigned int face_index)
const
2170 if (((dim == 2) && (shape_index >= this->first_quad_index)) ||
2171 ((dim == 3) && (shape_index >= this->first_hex_index)))
2176 if (shape_index < this->first_line_index)
2183 const unsigned int vertex_no = shape_index;
2186 for (
unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
2192 else if (shape_index < this->first_quad_index)
2195 const unsigned int line_index =
2196 (shape_index - this->first_line_index) / this->dofs_per_line;
2200 for (
unsigned int i = 0; i < GeometryInfo<dim>::lines_per_face; ++i)
2205 else if (shape_index < this->first_hex_index)
2208 const unsigned int quad_index =
2209 (shape_index - this->first_quad_index) / this->dofs_per_quad;
2210 Assert(
static_cast<signed int>(quad_index) <
2223 return (quad_index == face_index);
2245 std::vector<unsigned int>
2248 Assert((sub_degree > 0) && (sub_degree <= this->degree),
2253 std::vector<unsigned int> embedding_dofs(sub_degree + 1);
2254 for (
unsigned int i = 0; i < (sub_degree + 1); ++i)
2255 embedding_dofs[i] = i;
2257 return embedding_dofs;
2260 if (sub_degree == 1)
2262 std::vector<unsigned int> embedding_dofs(
2265 embedding_dofs[i] = i;
2267 return embedding_dofs;
2269 else if (sub_degree == this->degree)
2271 std::vector<unsigned int> embedding_dofs(this->dofs_per_cell);
2272 for (
unsigned int i = 0; i < this->dofs_per_cell; ++i)
2273 embedding_dofs[i] = i;
2275 return embedding_dofs;
2278 if ((dim == 2) || (dim == 3))
2280 std::vector<unsigned int> embedding_dofs(
2281 (dim == 2) ? (sub_degree + 1) * (sub_degree + 1) :
2282 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2284 for (
unsigned int i = 0;
2286 (sub_degree + 1) * (sub_degree + 1) :
2287 (sub_degree + 1) * (sub_degree + 1) * (sub_degree + 1));
2292 embedding_dofs[i] = i;
2297 const unsigned int j =
2299 const unsigned int line =
2304 line * (this->degree - 1) + j;
2312 const unsigned int j =
2316 const unsigned int k =
2321 const unsigned int face =
2324 k * (sub_degree - 1) - j) /
2325 ((sub_degree - 1) * (sub_degree - 1));
2330 face * (this->degree - 1) * (this->degree - 1) +
2331 k * (this->degree - 1) + j;
2339 (sub_degree - 1) * (sub_degree - 1))
2341 const unsigned int j =
2347 const unsigned int k =
2355 const unsigned int l =
2360 j - k * (sub_degree - 1)) /
2361 ((sub_degree - 1) * (sub_degree - 1));
2367 (this->degree - 1) +
2368 l * (this->degree - 1) * (this->degree - 1) +
2369 k * (this->degree - 1) + j;
2373 return embedding_dofs;
2378 return std::vector<unsigned int>();
2385 std::pair<Table<2, bool>, std::vector<unsigned int>>
2390 constant_modes(0, i) =
true;
2392 i < this->dofs_per_cell;
2394 constant_modes(0, i) =
false;
2395 return std::pair<Table<2, bool>, std::vector<unsigned int>>(
2396 constant_modes, std::vector<unsigned int>(1, 0));
2412 #include "fe_q_hierarchical.inst"