514 const unsigned int dofs_1d =
515 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
517 this->interface_constraints.TableBase<2, double>::reinit(
518 this->interface_constraints_size());
531 for (
unsigned int i = 0; i < dofs_1d; ++i)
532 this->interface_constraints(0, i) = dofs_subcell[0](1, i);
534 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
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);
545 for (
unsigned int i = 0; i < dofs_1d * dofs_1d; ++i)
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);
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);
567 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
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,
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,
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,
582 dofs_subcell[1](2 + j, i % dofs_1d) *
583 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
587 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
590 this->interface_constraints(5 + 4 * (this->degree - 1) + j,
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,
597 dofs_subcell[0](0, i % dofs_1d) *
598 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
600 this->interface_constraints(5 + 4 * (this->degree - 1) +
601 2 * (this->degree - 1) + j,
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,
608 dofs_subcell[1](1, i % dofs_1d) *
609 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
611 this->interface_constraints(5 + 4 * (this->degree - 1) +
612 4 * (this->degree - 1) + j,
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,
619 dofs_subcell[1](2 + j, i % dofs_1d) *
620 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
622 this->interface_constraints(5 + 4 * (this->degree - 1) +
623 6 * (this->degree - 1) + j,
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,
630 dofs_subcell[1](2 + j, i % dofs_1d) *
631 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
635 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
636 for (
unsigned int k = 0; k < (this->degree - 1); ++k)
639 this->interface_constraints(5 + 12 * (this->degree - 1) +
640 j + k * (this->degree - 1),
642 dofs_subcell[0](2 + j, i % dofs_1d) *
643 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
645 this->interface_constraints(5 + 12 * (this->degree - 1) +
646 j + k * (this->degree - 1) +
650 dofs_subcell[1](2 + j, i % dofs_1d) *
651 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
653 this->interface_constraints(5 + 12 * (this->degree - 1) +
654 j + k * (this->degree - 1) +
655 2 * (this->degree - 1) *
658 dofs_subcell[0](2 + j, i % dofs_1d) *
659 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
661 this->interface_constraints(5 + 12 * (this->degree - 1) +
662 j + k * (this->degree - 1) +
663 3 * (this->degree - 1) *
666 dofs_subcell[1](2 + j, i % dofs_1d) *
667 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
1022 const unsigned int subface,
1024 const unsigned int face_no)
const
1031 (
dynamic_cast<const FEQHierarchical *
>(&x_source_fe) !=
1035 Assert(interpolation_matrix.
n() == this->n_dofs_per_face(face_no),
1037 this->n_dofs_per_face(face_no)));
1069 interpolation_matrix(0, 0) = 1.0;
1070 interpolation_matrix(1, 0) = 0.5;
1071 interpolation_matrix(1, 1) = 0.5;
1073 for (
unsigned int dof = 2;
1074 dof < this->n_dofs_per_face(face_no);)
1076 interpolation_matrix(1, dof) = -1.0;
1080 int factorial_i = 1;
1082 for (
unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1085 interpolation_matrix(i, i) =
std::pow(0.5, i);
1087 int factorial_j = factorial_i;
1088 int factorial_ij = 1;
1090 for (
unsigned int j = i + 1;
1091 j < this->n_dofs_per_face(face_no);
1094 factorial_ij *= j - i;
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);
1103 interpolation_matrix(i, j) =
1105 (factorial_i * factorial_ij);
1114 interpolation_matrix(0, 0) = 0.5;
1115 interpolation_matrix(0, 1) = 0.5;
1117 for (
unsigned int dof = 2;
1118 dof < this->n_dofs_per_face(face_no);)
1120 interpolation_matrix(0, dof) = -1.0;
1124 interpolation_matrix(1, 1) = 1.0;
1126 int factorial_i = 1;
1128 for (
unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1131 interpolation_matrix(i, i) =
std::pow(0.5, i);
1133 int factorial_j = factorial_i;
1134 int factorial_ij = 1;
1136 for (
unsigned int j = i + 1;
1137 j < this->n_dofs_per_face(face_no);
1140 factorial_ij *= j - i;
1142 interpolation_matrix(i, j) =
1144 (factorial_i * factorial_ij);
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;
1165 for (
unsigned int i = 0; i < this->degree - 1;)
1167 for (
unsigned int line = 0;
1168 line < GeometryInfo<3>::lines_per_face;
1170 interpolation_matrix(3,
1171 i + line * (this->degree - 1) +
1174 for (
unsigned int j = 0; j < this->degree - 1;)
1176 interpolation_matrix(3,
1177 i + (j + 4) * this->degree - j) =
1182 interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1184 interpolation_matrix(2, i + 4) = -1.0;
1188 for (
unsigned int vertex = 0;
1189 vertex < GeometryInfo<3>::vertices_per_face;
1191 interpolation_matrix(3, vertex) = 0.25;
1193 int factorial_i = 1;
1195 for (
unsigned int i = 2; i <= this->degree; ++i)
1198 interpolation_matrix(i + 2, i + 2) = tmp;
1199 interpolation_matrix(i + 2 * source_fe.
degree,
1200 i + 2 * this->degree) = tmp;
1202 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
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;
1212 for (
unsigned int j = 0; j < this->degree - 1;)
1214 interpolation_matrix(i + source_fe.
degree + 1,
1215 (i + 2) * this->degree + j + 2 -
1217 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1218 i + (j + 4) * this->degree - j -
1223 int factorial_k = 1;
1225 for (
unsigned int j = 2; j <= this->degree; ++j)
1227 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1229 i + (j + 2) * this->degree - j) =
1232 int factorial_kl = 1;
1233 int factorial_l = factorial_k;
1235 for (
unsigned int k = j + 1; k < this->degree; ++k)
1237 factorial_kl *= k - j;
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);
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);
1257 int factorial_j = factorial_i;
1258 int factorial_ij = 1;
1260 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1262 factorial_ij *= j - i;
1265 if (((i + j) & 1) != 0u)
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;
1274 for (
unsigned int k = 2; k <= this->degree; ++k)
1276 interpolation_matrix(
1277 i + (k + 2) * source_fe.
degree - k,
1278 j + (k + 2) * this->degree - k) =
1281 int factorial_l = factorial_k;
1282 int factorial_kl = 1;
1284 for (
unsigned int l = k + 1;
1288 factorial_kl *= l - k;
1291 if (((k + l) & 1) != 0u)
1292 interpolation_matrix(
1293 i + (k + 2) * source_fe.
degree - k,
1294 j + (l + 2) * this->degree - l) =
1297 (factorial_k * factorial_kl);
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);
1309 interpolation_matrix(i + source_fe.
degree + 1,
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) =
1320 for (
unsigned int k = 0; k < this->degree - 1;)
1322 interpolation_matrix(i + source_fe.
degree + 1,
1323 (j + 2) * this->degree +
1325 interpolation_matrix(
1326 i + 3 * source_fe.
degree - 1,
1327 j + (k + 4) * this->degree - k - 2) = tmp;
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;
1340 for (
unsigned int k = 2; k <= this->degree; ++k)
1342 interpolation_matrix(
1343 i + (k + 2) * source_fe.
degree - k,
1344 j + (k + 2) * this->degree - k) =
1347 int factorial_l = factorial_k;
1348 int factorial_kl = 1;
1350 for (
unsigned int l = k + 1;
1354 factorial_kl *= l - k;
1357 if (((k + l) & 1) != 0u)
1358 interpolation_matrix(
1359 i + (k + 2) * source_fe.
degree - k,
1360 j + (l + 2) * this->degree - l) =
1363 (factorial_k * factorial_kl);
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);
1375 interpolation_matrix(i + source_fe.
degree + 1,
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) =
1386 for (
unsigned int k = 0; k < this->degree - 1;)
1388 interpolation_matrix(i + source_fe.
degree + 1,
1389 (j + 2) * this->degree +
1391 interpolation_matrix(
1392 i + 3 * source_fe.
degree - 1,
1393 j + (k + 4) * this->degree - k - 2) = tmp;
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;
1411 for (
unsigned int i = 0; i < this->degree - 1;)
1413 for (
unsigned int line = 0;
1414 line < GeometryInfo<3>::lines_per_face;
1416 interpolation_matrix(2,
1417 i + line * (this->degree - 1) +
1420 for (
unsigned int j = 0; j < this->degree - 1;)
1422 interpolation_matrix(2,
1423 i + (j + 4) * this->degree - j) =
1428 interpolation_matrix(0, i + 2 * (this->degree + 1)) =
1430 interpolation_matrix(3, i + this->degree + 3) = -1.0;
1434 for (
unsigned int vertex = 0;
1435 vertex < GeometryInfo<3>::vertices_per_face;
1437 interpolation_matrix(2, vertex) = 0.25;
1439 int factorial_i = 1;
1441 for (
unsigned int i = 2; i <= this->degree; ++i)
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;
1452 for (
unsigned int j = 0; j < this->degree - 1;)
1454 interpolation_matrix(i + 2,
1455 j + (i + 2) * this->degree + 2 -
1457 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1458 i + (j + 4) * this->degree - j -
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;
1469 int factorial_j = factorial_i;
1470 int factorial_ij = 1;
1472 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1474 factorial_ij *= j - i;
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;
1482 for (
unsigned int k = 2; k <= this->degree; ++k)
1484 interpolation_matrix(
1485 i + (k + 2) * source_fe.
degree - k,
1486 j + (k + 2) * this->degree - k) =
1489 int factorial_l = factorial_k;
1490 int factorial_kl = 1;
1492 for (
unsigned int l = k + 1; l <= this->degree;
1495 factorial_kl *= l - k;
1498 if (((k + l) & 1) != 0u)
1499 interpolation_matrix(
1500 i + (k + 2) * source_fe.
degree - k,
1501 j + (l + 2) * this->degree - l) =
1504 (factorial_k * factorial_kl);
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);
1517 for (
unsigned int k = 0; k < this->degree - 1;)
1519 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1520 j + (k + 4) * this->degree -
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;
1531 if (((i + j) & 1) != 0u)
1534 interpolation_matrix(i + 2, j + 2) = tmp;
1535 interpolation_matrix(i + 2, j + this->degree + 1) =
1537 interpolation_matrix(i + source_fe.
degree + 1,
1538 j + this->degree + 1) =
1542 for (
unsigned int k = 0; k < this->degree - 1;)
1544 interpolation_matrix(i + 2,
1545 k + (j + 2) * this->degree +
1551 int factorial_k = 1;
1553 for (
unsigned int j = 2; j <= this->degree; ++j)
1555 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1557 i + (j + 2) * this->degree - j) =
1560 int factorial_l = factorial_k;
1561 int factorial_kl = 1;
1563 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1565 factorial_kl *= k - j;
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);
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);
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;
1596 for (
unsigned int i = 0; i < this->degree - 1;)
1598 for (
unsigned int line = 0;
1599 line < GeometryInfo<3>::lines_per_face;
1601 interpolation_matrix(1,
1602 i + line * (this->degree - 1) +
1605 for (
unsigned int j = 0; j < this->degree - 1;)
1607 interpolation_matrix(1,
1608 i + (j + 4) * this->degree - j) =
1613 interpolation_matrix(0, i + 4) = -1.0;
1614 interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1618 for (
unsigned int vertex = 0;
1619 vertex < GeometryInfo<3>::vertices_per_face;
1621 interpolation_matrix(1, vertex) = 0.25;
1623 int factorial_i = 1;
1625 for (
unsigned int i = 2; i <= this->degree; ++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;
1632 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
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;
1642 for (
unsigned int j = 0; j < this->degree - 1;)
1644 interpolation_matrix(i + source_fe.
degree + 1,
1645 j + (i + 2) * this->degree + 2 -
1647 interpolation_matrix(i + 2 * source_fe.
degree,
1648 i + (j + 4) * this->degree - j -
1653 int factorial_k = 1;
1655 for (
unsigned int j = 2; j <= this->degree; ++j)
1657 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1659 i + (j + 2) * this->degree - j) =
1662 int factorial_kl = 1;
1663 int factorial_l = factorial_k;
1665 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1667 factorial_kl *= k - j;
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);
1678 int factorial_j = factorial_i;
1679 int factorial_ij = 1;
1681 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1683 factorial_ij *= j - i;
1685 tmp =
std::pow(0.5, j) * factorial_j /
1686 (factorial_i * factorial_ij);
1687 interpolation_matrix(i + 2, j + 2) = tmp;
1690 for (
unsigned int k = 0; k < this->degree - 1;)
1692 interpolation_matrix(i + source_fe.
degree + 1,
1693 k + (j + 2) * this->degree +
1699 interpolation_matrix(i + source_fe.
degree + 1,
1701 interpolation_matrix(i + source_fe.
degree + 1,
1702 j + this->degree + 1) = tmp;
1704 if (((i + j) & 1) != 0u)
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;
1712 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1713 j + 3 * this->degree - 1) = tmp;
1714 int factorial_k = 1;
1716 for (
unsigned int k = 2; k <= this->degree; ++k)
1718 interpolation_matrix(
1719 i + (k + 2) * source_fe.
degree - k,
1720 j + (k + 2) * this->degree - k) =
1723 int factorial_l = factorial_k;
1724 int factorial_kl = 1;
1726 for (
unsigned int l = k + 1; l <= this->degree;
1729 factorial_kl *= l - k;
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);
1741 for (
unsigned int k = 0; k < this->degree - 1;)
1743 interpolation_matrix(i + 2 * source_fe.
degree,
1744 j + (k + 4) * this->degree -
1756 for (
unsigned int vertex = 0;
1757 vertex < GeometryInfo<3>::vertices_per_face;
1759 interpolation_matrix(0, vertex) = 0.25;
1761 for (
unsigned int i = 0; i < this->degree - 1;)
1763 for (
unsigned int line = 0;
1764 line < GeometryInfo<3>::lines_per_face;
1766 interpolation_matrix(0,
1767 i + line * (this->degree - 1) +
1770 for (
unsigned int j = 0; j < this->degree - 1;)
1772 interpolation_matrix(0,
1773 i + (j + 4) * this->degree - j) =
1778 interpolation_matrix(1, i + 4) = -1.0;
1779 interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
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;
1789 int factorial_i = 1;
1791 for (
unsigned int i = 2; i <= this->degree; ++i)
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;
1802 for (
unsigned int j = 0; j < this->degree - 1;)
1804 interpolation_matrix(i + 2,
1805 j + (i + 2) * this->degree + 2 -
1807 interpolation_matrix(i + 2 * source_fe.
degree,
1808 i + (j + 4) * this->degree - 2) =
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;
1820 for (
unsigned int j = 2; j <= this->degree; ++j)
1822 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1824 i + (j + 2) * this->degree - j) =
1827 int factorial_l = factorial_k;
1828 int factorial_kl = 1;
1830 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1832 factorial_kl *= k - j;
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);
1843 int factorial_j = factorial_i;
1844 int factorial_ij = 1;
1846 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1848 factorial_ij *= j - i;
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) =
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;
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;
1866 for (
unsigned int k = 2; k <= this->degree; ++k)
1868 interpolation_matrix(
1869 i + (k + 2) * source_fe.
degree - k,
1870 j + (k + 2) * this->degree - k) =
1873 int factorial_l = factorial_k;
1874 int factorial_kl = 1;
1876 for (
unsigned int l = k + 1; l <= this->degree;
1879 factorial_kl *= l - k;
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);
1891 for (
unsigned int k = 0; k < this->degree - 1;)
1893 interpolation_matrix(i + 2,
1894 k + (j + 2) * this->degree +
1896 interpolation_matrix(i + 2 * source_fe.
degree,
1897 j + (k + 4) * this->degree -