512 const unsigned int dofs_1d =
513 2 * this->n_dofs_per_vertex() + this->n_dofs_per_line();
515 this->interface_constraints.TableBase<2, double>::reinit(
516 this->interface_constraints_size());
529 for (
unsigned int i = 0; i < dofs_1d; ++i)
530 this->interface_constraints(0, i) = dofs_subcell[0](1, i);
532 for (
unsigned int c = 0; c < GeometryInfo<1>::max_children_per_cell;
534 for (
unsigned int i = 0; i < dofs_1d; ++i)
535 for (
unsigned int j = 2; j < dofs_1d; ++j)
536 this->interface_constraints(1 + c * (this->degree - 1) + j - 2,
537 i) = dofs_subcell[c](j, i);
543 for (
unsigned int i = 0; i < dofs_1d * dofs_1d; ++i)
546 this->interface_constraints(0, face_renumber[i]) =
547 dofs_subcell[0](1, i % dofs_1d) *
548 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
551 this->interface_constraints(1, face_renumber[i]) =
552 dofs_subcell[0](0, i % dofs_1d) *
553 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
554 this->interface_constraints(2, face_renumber[i]) =
555 dofs_subcell[1](1, i % dofs_1d) *
556 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
557 this->interface_constraints(3, face_renumber[i]) =
558 dofs_subcell[0](1, i % dofs_1d) *
559 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
560 this->interface_constraints(4, face_renumber[i]) =
561 dofs_subcell[1](0, i % dofs_1d) *
562 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
565 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
567 this->interface_constraints(5 + j, face_renumber[i]) =
568 dofs_subcell[0](1, i % dofs_1d) *
569 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
570 this->interface_constraints(5 + (this->degree - 1) + j,
572 dofs_subcell[0](1, i % dofs_1d) *
573 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
574 this->interface_constraints(5 + 2 * (this->degree - 1) + j,
576 dofs_subcell[0](2 + j, i % dofs_1d) *
577 dofs_subcell[1](0, (i - (i % dofs_1d)) / dofs_1d);
578 this->interface_constraints(5 + 3 * (this->degree - 1) + j,
580 dofs_subcell[1](2 + j, i % dofs_1d) *
581 dofs_subcell[0](1, (i - (i % dofs_1d)) / dofs_1d);
585 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
588 this->interface_constraints(5 + 4 * (this->degree - 1) + j,
590 dofs_subcell[0](0, i % dofs_1d) *
591 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
592 this->interface_constraints(5 + 4 * (this->degree - 1) +
593 (this->degree - 1) + j,
595 dofs_subcell[0](0, 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 2 * (this->degree - 1) + j,
601 dofs_subcell[1](1, i % dofs_1d) *
602 dofs_subcell[0](2 + j, (i - (i % dofs_1d)) / dofs_1d);
603 this->interface_constraints(5 + 4 * (this->degree - 1) +
604 3 * (this->degree - 1) + j,
606 dofs_subcell[1](1, i % dofs_1d) *
607 dofs_subcell[1](2 + j, (i - (i % dofs_1d)) / dofs_1d);
609 this->interface_constraints(5 + 4 * (this->degree - 1) +
610 4 * (this->degree - 1) + j,
612 dofs_subcell[0](2 + j, i % dofs_1d) *
613 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
614 this->interface_constraints(5 + 4 * (this->degree - 1) +
615 5 * (this->degree - 1) + j,
617 dofs_subcell[1](2 + j, i % dofs_1d) *
618 dofs_subcell[0](0, (i - (i % dofs_1d)) / dofs_1d);
620 this->interface_constraints(5 + 4 * (this->degree - 1) +
621 6 * (this->degree - 1) + j,
623 dofs_subcell[0](2 + j, i % dofs_1d) *
624 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
625 this->interface_constraints(5 + 4 * (this->degree - 1) +
626 7 * (this->degree - 1) + j,
628 dofs_subcell[1](2 + j, i % dofs_1d) *
629 dofs_subcell[1](1, (i - (i % dofs_1d)) / dofs_1d);
633 for (
unsigned int j = 0; j < (this->degree - 1); ++j)
634 for (
unsigned int k = 0; k < (this->degree - 1); ++k)
637 this->interface_constraints(5 + 12 * (this->degree - 1) +
638 j + k * (this->degree - 1),
640 dofs_subcell[0](2 + j, i % dofs_1d) *
641 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
643 this->interface_constraints(5 + 12 * (this->degree - 1) +
644 j + k * (this->degree - 1) +
648 dofs_subcell[1](2 + j, i % dofs_1d) *
649 dofs_subcell[0](2 + k, (i - (i % dofs_1d)) / dofs_1d);
651 this->interface_constraints(5 + 12 * (this->degree - 1) +
652 j + k * (this->degree - 1) +
653 2 * (this->degree - 1) *
656 dofs_subcell[0](2 + j, i % dofs_1d) *
657 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
659 this->interface_constraints(5 + 12 * (this->degree - 1) +
660 j + k * (this->degree - 1) +
661 3 * (this->degree - 1) *
664 dofs_subcell[1](2 + j, i % dofs_1d) *
665 dofs_subcell[1](2 + k, (i - (i % dofs_1d)) / dofs_1d);
1020 const unsigned int subface,
1022 const unsigned int face_no)
const
1029 (
dynamic_cast<const FEQHierarchical *
>(&x_source_fe) !=
1033 Assert(interpolation_matrix.
n() == this->n_dofs_per_face(face_no),
1035 this->n_dofs_per_face(face_no)));
1067 interpolation_matrix(0, 0) = 1.0;
1068 interpolation_matrix(1, 0) = 0.5;
1069 interpolation_matrix(1, 1) = 0.5;
1071 for (
unsigned int dof = 2;
1072 dof < this->n_dofs_per_face(face_no);)
1074 interpolation_matrix(1, dof) = -1.0;
1078 int factorial_i = 1;
1080 for (
unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1085 int factorial_j = factorial_i;
1086 int factorial_ij = 1;
1088 for (
unsigned int j = i + 1;
1089 j < this->n_dofs_per_face(face_no);
1092 factorial_ij *= j - i;
1095 if (((i + j) & 1) != 0u)
1096 interpolation_matrix(i, j) =
1098 (factorial_i * factorial_ij);
1101 interpolation_matrix(i, j) =
1103 (factorial_i * factorial_ij);
1112 interpolation_matrix(0, 0) = 0.5;
1113 interpolation_matrix(0, 1) = 0.5;
1115 for (
unsigned int dof = 2;
1116 dof < this->n_dofs_per_face(face_no);)
1118 interpolation_matrix(0, dof) = -1.0;
1122 interpolation_matrix(1, 1) = 1.0;
1124 int factorial_i = 1;
1126 for (
unsigned int i = 2; i < this->n_dofs_per_face(face_no);
1131 int factorial_j = factorial_i;
1132 int factorial_ij = 1;
1134 for (
unsigned int j = i + 1;
1135 j < this->n_dofs_per_face(face_no);
1138 factorial_ij *= j - i;
1140 interpolation_matrix(i, j) =
1142 (factorial_i * factorial_ij);
1157 interpolation_matrix(0, 0) = 1.0;
1158 interpolation_matrix(1, 0) = 0.5;
1159 interpolation_matrix(1, 1) = 0.5;
1160 interpolation_matrix(2, 0) = 0.5;
1161 interpolation_matrix(2, 2) = 0.5;
1163 for (
unsigned int i = 0; i < this->degree - 1;)
1165 for (
unsigned int line = 0;
1166 line < GeometryInfo<3>::lines_per_face;
1168 interpolation_matrix(3,
1169 i + line * (this->degree - 1) +
1172 for (
unsigned int j = 0; j < this->degree - 1;)
1174 interpolation_matrix(3,
1175 i + (j + 4) * this->degree - j) =
1180 interpolation_matrix(1, i + 2 * (this->degree + 1)) =
1182 interpolation_matrix(2, i + 4) = -1.0;
1186 for (
unsigned int vertex = 0;
1187 vertex < GeometryInfo<3>::vertices_per_face;
1189 interpolation_matrix(3, vertex) = 0.25;
1191 int factorial_i = 1;
1193 for (
unsigned int i = 2; i <= this->degree; ++i)
1196 interpolation_matrix(i + 2, i + 2) = tmp;
1197 interpolation_matrix(i + 2 * source_fe.
degree,
1198 i + 2 * this->degree) = tmp;
1200 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
1202 interpolation_matrix(i + source_fe.
degree + 1,
1203 i + this->degree + 1) = tmp;
1204 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1205 i + 2 * this->degree) = tmp;
1206 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1207 i + 3 * this->degree - 1) = tmp;
1210 for (
unsigned int j = 0; j < this->degree - 1;)
1212 interpolation_matrix(i + source_fe.
degree + 1,
1213 (i + 2) * this->degree + j + 2 -
1215 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1216 i + (j + 4) * this->degree - j -
1221 int factorial_k = 1;
1223 for (
unsigned int j = 2; j <= this->degree; ++j)
1225 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1227 i + (j + 2) * this->degree - j) =
1230 int factorial_kl = 1;
1231 int factorial_l = factorial_k;
1233 for (
unsigned int k = j + 1; k < this->degree; ++k)
1235 factorial_kl *= k - j;
1238 if (((j + k) & 1) != 0u)
1239 interpolation_matrix(
1240 i + (j + 2) * source_fe.
degree - j,
1241 i + (k + 2) * this->degree - k) =
1243 factorial_l / (factorial_k * factorial_kl);
1246 interpolation_matrix(
1247 i + (j + 2) * source_fe.
degree - j,
1248 i + (k + 2) * this->degree - k) =
1250 (factorial_k * factorial_kl);
1255 int factorial_j = factorial_i;
1256 int factorial_ij = 1;
1258 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1260 factorial_ij *= j - i;
1263 if (((i + j) & 1) != 0u)
1266 factorial_j / (factorial_i * factorial_ij);
1267 interpolation_matrix(i + 2, j + 2) = tmp;
1268 interpolation_matrix(i + 2 * source_fe.
degree,
1269 j + 2 * this->degree) = tmp;
1272 for (
unsigned int k = 2; k <= this->degree; ++k)
1274 interpolation_matrix(
1275 i + (k + 2) * source_fe.
degree - k,
1276 j + (k + 2) * this->degree - k) =
1279 int factorial_l = factorial_k;
1280 int factorial_kl = 1;
1282 for (
unsigned int l = k + 1;
1286 factorial_kl *= l - k;
1289 if (((k + l) & 1) != 0u)
1290 interpolation_matrix(
1291 i + (k + 2) * source_fe.
degree - k,
1292 j + (l + 2) * this->degree - l) =
1295 (factorial_k * factorial_kl);
1298 interpolation_matrix(
1299 i + (k + 2) * source_fe.
degree - k,
1300 j + (l + 2) * this->degree - l) =
1303 (factorial_k * factorial_kl);
1308 interpolation_matrix(i + source_fe.
degree + 1,
1310 interpolation_matrix(i + source_fe.
degree + 1,
1311 j + this->degree + 1) = tmp;
1312 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1313 j + 2 * this->degree) = tmp;
1314 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1315 j + 3 * this->degree - 1) =
1319 for (
unsigned int k = 0; k < this->degree - 1;)
1321 interpolation_matrix(i + source_fe.
degree + 1,
1322 (j + 2) * this->degree +
1324 interpolation_matrix(
1325 i + 3 * source_fe.
degree - 1,
1326 j + (k + 4) * this->degree - k - 2) = tmp;
1333 (factorial_i * factorial_ij);
1334 interpolation_matrix(i + 2, j + 2) = tmp;
1335 interpolation_matrix(i + 2 * source_fe.
degree,
1336 j + 2 * this->degree) = tmp;
1339 for (
unsigned int k = 2; k <= this->degree; ++k)
1341 interpolation_matrix(
1342 i + (k + 2) * source_fe.
degree - k,
1343 j + (k + 2) * this->degree - k) =
1346 int factorial_l = factorial_k;
1347 int factorial_kl = 1;
1349 for (
unsigned int l = k + 1;
1353 factorial_kl *= l - k;
1356 if (((k + l) & 1) != 0u)
1357 interpolation_matrix(
1358 i + (k + 2) * source_fe.
degree - k,
1359 j + (l + 2) * this->degree - l) =
1362 (factorial_k * factorial_kl);
1365 interpolation_matrix(
1366 i + (k + 2) * source_fe.
degree - k,
1367 j + (l + 2) * this->degree - 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;
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) =
1512 (factorial_k * factorial_kl);
1518 for (
unsigned int k = 0; k < this->degree - 1;)
1520 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1521 j + (k + 4) * this->degree -
1527 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1528 j + 2 * this->degree) = tmp;
1529 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1530 j + 3 * this->degree - 1) = tmp;
1532 if (((i + j) & 1) != 0u)
1535 interpolation_matrix(i + 2, j + 2) = tmp;
1536 interpolation_matrix(i + 2, j + this->degree + 1) =
1538 interpolation_matrix(i + source_fe.
degree + 1,
1539 j + this->degree + 1) =
1543 for (
unsigned int k = 0; k < this->degree - 1;)
1545 interpolation_matrix(i + 2,
1546 k + (j + 2) * this->degree +
1552 int factorial_k = 1;
1554 for (
unsigned int j = 2; j <= this->degree; ++j)
1556 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1558 i + (j + 2) * this->degree - j) =
1561 int factorial_l = factorial_k;
1562 int factorial_kl = 1;
1564 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1566 factorial_kl *= k - j;
1569 if (((j + k) & 1) != 0u)
1570 interpolation_matrix(
1571 i + (j + 2) * source_fe.
degree - j,
1572 i + (k + 2) * this->degree - k) =
1574 factorial_l / (factorial_k * factorial_kl);
1577 interpolation_matrix(
1578 i + (j + 2) * source_fe.
degree - j,
1579 i + (k + 2) * this->degree - k) =
1581 (factorial_k * factorial_kl);
1591 interpolation_matrix(0, 0) = 0.5;
1592 interpolation_matrix(0, 2) = 0.5;
1593 interpolation_matrix(2, 2) = 1.0;
1594 interpolation_matrix(3, 2) = 0.5;
1595 interpolation_matrix(3, 3) = 0.5;
1597 for (
unsigned int i = 0; i < this->degree - 1;)
1599 for (
unsigned int line = 0;
1600 line < GeometryInfo<3>::lines_per_face;
1602 interpolation_matrix(1,
1603 i + line * (this->degree - 1) +
1606 for (
unsigned int j = 0; j < this->degree - 1;)
1608 interpolation_matrix(1,
1609 i + (j + 4) * this->degree - j) =
1614 interpolation_matrix(0, i + 4) = -1.0;
1615 interpolation_matrix(3, i + 3 * this->degree + 1) = -1.0;
1619 for (
unsigned int vertex = 0;
1620 vertex < GeometryInfo<3>::vertices_per_face;
1622 interpolation_matrix(1, vertex) = 0.25;
1624 int factorial_i = 1;
1626 for (
unsigned int i = 2; i <= this->degree; ++i)
1629 interpolation_matrix(i + 2, i + 2) = tmp;
1630 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1631 i + 3 * this->degree - 1) = tmp;
1633 interpolation_matrix(i + source_fe.
degree + 1, i + 2) =
1635 interpolation_matrix(i + source_fe.
degree + 1,
1636 i + this->degree + 1) = tmp;
1637 interpolation_matrix(i + 2 * source_fe.
degree,
1638 i + 2 * this->degree) = tmp;
1639 interpolation_matrix(i + 2 * source_fe.
degree,
1640 i + 3 * this->degree - 1) = tmp;
1643 for (
unsigned int j = 0; j < this->degree - 1;)
1645 interpolation_matrix(i + source_fe.
degree + 1,
1646 j + (i + 2) * this->degree + 2 -
1648 interpolation_matrix(i + 2 * source_fe.
degree,
1649 i + (j + 4) * this->degree - j -
1654 int factorial_k = 1;
1656 for (
unsigned int j = 2; j <= this->degree; ++j)
1658 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1660 i + (j + 2) * this->degree - j) =
1663 int factorial_kl = 1;
1664 int factorial_l = factorial_k;
1666 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1668 factorial_kl *= k - j;
1670 interpolation_matrix(
1671 i + (j + 2) * source_fe.
degree - j,
1672 i + (k + 2) * this->degree - k) =
1674 (factorial_k * factorial_kl);
1679 int factorial_j = factorial_i;
1680 int factorial_ij = 1;
1682 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1684 factorial_ij *= j - i;
1687 (factorial_i * factorial_ij);
1688 interpolation_matrix(i + 2, j + 2) = tmp;
1691 for (
unsigned int k = 0; k < this->degree - 1;)
1693 interpolation_matrix(i + source_fe.
degree + 1,
1694 k + (j + 2) * this->degree +
1700 interpolation_matrix(i + source_fe.
degree + 1,
1702 interpolation_matrix(i + source_fe.
degree + 1,
1703 j + this->degree + 1) = tmp;
1705 if (((i + j) & 1) != 0u)
1708 interpolation_matrix(i + 2 * source_fe.
degree,
1709 j + 2 * this->degree) = tmp;
1710 interpolation_matrix(i + 2 * source_fe.
degree,
1711 j + 3 * this->degree - 1) = tmp;
1713 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1714 j + 3 * this->degree - 1) = tmp;
1715 int factorial_k = 1;
1717 for (
unsigned int k = 2; k <= this->degree; ++k)
1719 interpolation_matrix(
1720 i + (k + 2) * source_fe.
degree - k,
1721 j + (k + 2) * this->degree - k) =
1724 int factorial_l = factorial_k;
1725 int factorial_kl = 1;
1727 for (
unsigned int l = k + 1; l <= this->degree;
1730 factorial_kl *= l - k;
1732 interpolation_matrix(
1733 i + (k + 2) * source_fe.
degree - k,
1734 j + (l + 2) * this->degree - l) =
1736 (factorial_k * factorial_kl);
1742 for (
unsigned int k = 0; k < this->degree - 1;)
1744 interpolation_matrix(i + 2 * source_fe.
degree,
1745 j + (k + 4) * this->degree -
1757 for (
unsigned int vertex = 0;
1758 vertex < GeometryInfo<3>::vertices_per_face;
1760 interpolation_matrix(0, vertex) = 0.25;
1762 for (
unsigned int i = 0; i < this->degree - 1;)
1764 for (
unsigned int line = 0;
1765 line < GeometryInfo<3>::lines_per_face;
1767 interpolation_matrix(0,
1768 i + line * (this->degree - 1) +
1771 for (
unsigned int j = 0; j < this->degree - 1;)
1773 interpolation_matrix(0,
1774 i + (j + 4) * this->degree - j) =
1779 interpolation_matrix(1, i + 4) = -1.0;
1780 interpolation_matrix(2, i + 3 * this->degree + 1) = -1.0;
1784 interpolation_matrix(1, 0) = 0.5;
1785 interpolation_matrix(1, 1) = 0.5;
1786 interpolation_matrix(2, 2) = 0.5;
1787 interpolation_matrix(2, 3) = 0.5;
1788 interpolation_matrix(3, 3) = 1.0;
1790 int factorial_i = 1;
1792 for (
unsigned int i = 2; i <= this->degree; ++i)
1795 interpolation_matrix(i + 2, i + 2) = tmp;
1796 interpolation_matrix(i + 2, i + this->degree + 1) = tmp;
1797 interpolation_matrix(i + 2 * source_fe.
degree,
1798 i + 2 * this->degree) = tmp;
1799 interpolation_matrix(i + 2 * source_fe.
degree,
1800 i + 3 * this->degree - 1) = tmp;
1803 for (
unsigned int j = 0; j < this->degree - 1;)
1805 interpolation_matrix(i + 2,
1806 j + (i + 2) * this->degree + 2 -
1808 interpolation_matrix(i + 2 * source_fe.
degree,
1809 i + (j + 4) * this->degree - 2) =
1815 interpolation_matrix(i + source_fe.
degree + 1,
1816 i + this->degree + 1) = tmp;
1817 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1818 i + 3 * this->degree - 1) = tmp;
1819 int factorial_k = 1;
1821 for (
unsigned int j = 2; j <= this->degree; ++j)
1823 interpolation_matrix(i + (j + 2) * source_fe.
degree -
1825 i + (j + 2) * this->degree - j) =
1828 int factorial_l = factorial_k;
1829 int factorial_kl = 1;
1831 for (
unsigned int k = j + 1; k <= this->degree; ++k)
1833 factorial_kl *= k - j;
1835 interpolation_matrix(
1836 i + (j + 2) * source_fe.
degree - j,
1837 i + (k + 2) * this->degree - k) =
1839 (factorial_k * factorial_kl);
1844 int factorial_j = factorial_i;
1845 int factorial_ij = 1;
1847 for (
unsigned int j = i + 1; j <= this->degree; ++j)
1849 factorial_ij *= j - i;
1852 (factorial_i * factorial_ij);
1853 interpolation_matrix(i + 2, j + 2) = tmp;
1854 interpolation_matrix(i + 2, j + this->degree + 1) =
1856 interpolation_matrix(i + 2 * source_fe.
degree,
1857 j + 2 * this->degree) = tmp;
1858 interpolation_matrix(i + 2 * source_fe.
degree,
1859 j + 3 * this->degree - 1) = tmp;
1861 interpolation_matrix(i + source_fe.
degree + 1,
1862 j + this->degree + 1) = tmp;
1863 interpolation_matrix(i + 3 * source_fe.
degree - 1,
1864 j + 3 * this->degree - 1) = tmp;
1865 int factorial_k = 1;
1867 for (
unsigned int k = 2; k <= this->degree; ++k)
1869 interpolation_matrix(
1870 i + (k + 2) * source_fe.
degree - k,
1871 j + (k + 2) * this->degree - k) =
1874 int factorial_l = factorial_k;
1875 int factorial_kl = 1;
1877 for (
unsigned int l = k + 1; l <= this->degree;
1880 factorial_kl *= l - k;
1882 interpolation_matrix(
1883 i + (k + 2) * source_fe.
degree - k,
1884 j + (l + 2) * this->degree - l) =
1886 (factorial_k * factorial_kl);
1892 for (
unsigned int k = 0; k < this->degree - 1;)
1894 interpolation_matrix(i + 2,
1895 k + (j + 2) * this->degree +
1897 interpolation_matrix(i + 2 * source_fe.
degree,
1898 j + (k + 4) * this->degree -