576 const std::vector<Point<1>> &edge_quadrature_points =
577 edge_quadrature.get_points();
578 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
590 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
593 const double weight = 2.0 * edge_quadrature.weight(q_point);
595 if (edge_quadrature_points[q_point][0] < 0.5)
598 0.0, 2.0 * edge_quadrature_points[q_point][0]);
603 quadrature_point[0] = 1.0;
607 quadrature_point[0] = quadrature_point[1];
608 quadrature_point[1] = 0.0;
609 this->
restriction[index][0](2 * this->degree, dof) +=
612 quadrature_point[1] = 1.0;
613 this->
restriction[index][2](3 * this->degree, dof) +=
621 0.0, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
626 quadrature_point[0] = 1.0;
630 quadrature_point[0] = quadrature_point[1];
631 quadrature_point[1] = 0.0;
632 this->
restriction[index][1](2 * this->degree, dof) +=
635 quadrature_point[1] = 1.0;
636 this->
restriction[index][3](3 * this->degree, dof) +=
648 const unsigned int deg = this->
degree - 1;
649 const std::vector<Polynomials::Polynomial<double>>
650 &legendre_polynomials =
656 n_edge_quadrature_points);
658 for (
unsigned int q_point = 0;
659 q_point < n_edge_quadrature_points;
662 const double weight =
663 std::sqrt(edge_quadrature.weight(q_point));
665 for (
unsigned int i = 0; i < deg; ++i)
666 assembling_matrix(i, q_point) =
667 weight * legendre_polynomials[i + 1].value(
668 edge_quadrature_points[q_point][0]);
673 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
674 system_matrix_inv.
invert(system_matrix);
681 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
682 for (
unsigned int i = 0; i < 2; ++i)
686 for (
unsigned int q_point = 0;
687 q_point < n_edge_quadrature_points;
690 const double weight = edge_quadrature.weight(q_point);
692 i, edge_quadrature_points[q_point][0]);
694 edge_quadrature_points[q_point][0], i);
696 if (edge_quadrature_points[q_point][0] < 0.5)
699 i, 2.0 * edge_quadrature_points[q_point][0]);
704 dof, quadrature_point_2, 1) -
718 2.0 * edge_quadrature_points[q_point][0], i);
722 dof, quadrature_point_2, 0) -
723 this->restriction[index][2 * i]((i + 2) *
732 this->restriction[index][2 * i + 1](
733 (i + 2) * this->degree, dof) *
735 (i + 2) * this->degree, quadrature_point_1, 0);
750 2.0 * edge_quadrature_points[q_point][0] - 1.0);
755 dof, quadrature_point_2, 1) -
756 this->restriction[index][i + 2](i * this->degree,
763 this->restriction[index][2 * i]((i + 2) *
767 (i + 2) * this->degree, quadrature_point_1, 0);
769 2.0 * edge_quadrature_points[q_point][0] - 1.0,
774 dof, quadrature_point_2, 0) -
775 this->restriction[index][2 * i + 1](
776 (i + 2) * this->degree, dof) *
783 for (
unsigned int j = 0; j < this->
degree - 1; ++j)
786 legendre_polynomials[j + 1].value(
787 edge_quadrature_points[q_point][0]);
789 for (
unsigned int k = 0; k < tmp.
size(); ++k)
790 system_rhs(j, k) += tmp(k) * L_j;
794 system_matrix_inv.
mmult(solution, system_rhs);
796 for (
unsigned int j = 0; j < this->
degree - 1; ++j)
797 for (
unsigned int k = 0; k < 2; ++k)
799 if (
std::abs(solution(j, k)) > 1e-14)
801 i * this->degree + j + 1, dof) = solution(j, k);
803 if (
std::abs(solution(j, k + 2)) > 1e-14)
805 (i + 2) * this->degree + j + 1, dof) =
811 const std::vector<Point<dim>> &quadrature_points =
812 quadrature.get_points();
813 const std::vector<Polynomials::Polynomial<double>>
814 &lobatto_polynomials =
816 const unsigned int n_boundary_dofs =
818 const unsigned int n_quadrature_points = quadrature.size();
823 n_quadrature_points);
825 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
828 const double weight =
std::sqrt(quadrature.weight(q_point));
830 for (
unsigned int i = 0; i < this->degree; ++i)
833 weight * legendre_polynomials[i].value(
834 quadrature_points[q_point][0]);
836 for (
unsigned int j = 0; j < this->degree - 1; ++j)
837 assembling_matrix(i * (this->degree - 1) + j,
839 L_i * lobatto_polynomials[j + 2].value(
840 quadrature_points[q_point][1]);
845 assembling_matrix.
m());
847 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
848 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
849 system_matrix_inv.
invert(system_matrix);
852 solution.reinit(system_matrix_inv.
m(), 8);
853 system_rhs.reinit(system_matrix_inv.
m(), 8);
856 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
860 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
865 if (quadrature_points[q_point][0] < 0.5)
867 if (quadrature_points[q_point][1] < 0.5)
870 2.0 * quadrature_points[q_point][0],
871 2.0 * quadrature_points[q_point][1]);
874 dof, quadrature_point, 0);
876 dof, quadrature_point, 1);
882 2.0 * quadrature_points[q_point][0],
883 2.0 * quadrature_points[q_point][1] - 1.0);
886 dof, quadrature_point, 0);
888 dof, quadrature_point, 1);
892 else if (quadrature_points[q_point][1] < 0.5)
895 2.0 * quadrature_points[q_point][0] - 1.0,
896 2.0 * quadrature_points[q_point][1]);
911 2.0 * quadrature_points[q_point][0] - 1.0,
912 2.0 * quadrature_points[q_point][1] - 1.0);
924 for (
unsigned int i = 0; i < 2; ++i)
925 for (
unsigned int j = 0; j < this->degree; ++j)
931 j + 2 * this->degree,
932 quadrature_points[q_point],
938 i * this->degree + j,
939 quadrature_points[q_point],
941 tmp(2 * (i + 2)) -= this->
restriction[index][i + 2](
942 j + 3 * this->degree, dof) *
944 j + 3 * this->degree,
945 quadrature_points[q_point],
948 i * this->degree + j, dof) *
950 i * this->degree + j,
951 quadrature_points[q_point],
955 tmp *= quadrature.weight(q_point);
957 for (
unsigned int i = 0; i < this->degree; ++i)
959 const double L_i_0 = legendre_polynomials[i].value(
960 quadrature_points[q_point][0]);
961 const double L_i_1 = legendre_polynomials[i].value(
962 quadrature_points[q_point][1]);
964 for (
unsigned int j = 0; j < this->degree - 1; ++j)
967 L_i_0 * lobatto_polynomials[j + 2].value(
968 quadrature_points[q_point][1]);
970 L_i_1 * lobatto_polynomials[j + 2].value(
971 quadrature_points[q_point][0]);
973 for (
unsigned int k = 0; k < 4; ++k)
975 system_rhs(i * (this->degree - 1) + j,
976 2 * k) += tmp(2 * k) * l_j_0;
977 system_rhs(i * (this->degree - 1) + j,
979 tmp(2 * k + 1) * l_j_1;
985 system_matrix_inv.
mmult(solution, system_rhs);
987 for (
unsigned int i = 0; i < this->degree; ++i)
988 for (
unsigned int j = 0; j < this->degree - 1; ++j)
989 for (
unsigned int k = 0; k < 4; ++k)
991 if (
std::abs(solution(i * (this->degree - 1) + j,
993 this->
restriction[index][k](i * (this->degree - 1) +
996 solution(i * (this->degree - 1) + j, 2 * k);
998 if (
std::abs(solution(i * (this->degree - 1) + j,
1001 i + (this->degree - 1 + j) * this->degree +
1004 solution(i * (this->degree - 1) + j, 2 * k + 1);
1019 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
1022 const double weight = 2.0 * edge_quadrature.weight(q_point);
1024 if (edge_quadrature_points[q_point][0] < 0.5)
1025 for (
unsigned int i = 0; i < 2; ++i)
1026 for (
unsigned int j = 0; j < 2; ++j)
1029 i, 2.0 * edge_quadrature_points[q_point][0], j);
1037 Point<dim>(2.0 * edge_quadrature_points[q_point][0],
1041 (i + 4 * j + 2) * this->degree, dof) +=
1047 2.0 * edge_quadrature_points[q_point][0]);
1048 this->
restriction[index][i + 2 * j]((i + 2 * (j + 4)) *
1052 this->shape_value_component(dof, quadrature_point, 2);
1056 for (
unsigned int i = 0; i < 2; ++i)
1057 for (
unsigned int j = 0; j < 2; ++j)
1060 i, 2.0 * edge_quadrature_points[q_point][0] - 1.0, j);
1062 this->
restriction[index][i + 4 * j + 2]((i + 4 * j) *
1068 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j);
1070 (i + 4 * j + 2) * this->degree, dof) +=
1074 i, j, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1076 (i + 2 * (j + 4)) * this->degree, dof) +=
1078 this->shape_value_component(dof, quadrature_point, 2);
1088 const unsigned int deg = this->
degree - 1;
1089 const std::vector<Polynomials::Polynomial<double>>
1090 &legendre_polynomials =
1096 n_edge_quadrature_points);
1098 for (
unsigned int q_point = 0;
1099 q_point < n_edge_quadrature_points;
1102 const double weight =
1103 std::sqrt(edge_quadrature.weight(q_point));
1105 for (
unsigned int i = 0; i < deg; ++i)
1106 assembling_matrix(i, q_point) =
1107 weight * legendre_polynomials[i + 1].value(
1108 edge_quadrature_points[q_point][0]);
1113 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1114 system_matrix_inv.
invert(system_matrix);
1121 for (
unsigned int i = 0; i < 2; ++i)
1122 for (
unsigned int j = 0; j < 2; ++j)
1123 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell();
1128 for (
unsigned int q_point = 0;
1129 q_point < n_edge_quadrature_points;
1132 const double weight = edge_quadrature.weight(q_point);
1134 i, edge_quadrature_points[q_point][0], j);
1136 edge_quadrature_points[q_point][0], i, j);
1138 i, j, edge_quadrature_points[q_point][0]);
1140 if (edge_quadrature_points[q_point][0] < 0.5)
1143 i, 2.0 * edge_quadrature_points[q_point][0], j);
1147 dof, quadrature_point_3, 1) -
1149 (i + 4 * j) * this->
degree, dof) *
1151 (i + 4 * j) * this->degree,
1157 (i + 4 * j) * this->degree, dof) *
1163 2.0 * edge_quadrature_points[q_point][0], i, j);
1167 dof, quadrature_point_3, 0) -
1168 this->restriction[index][2 * (i + 2 * j)](
1169 (i + 4 * j + 2) * this->degree, dof) *
1171 (i + 4 * j + 2) * this->degree,
1176 this->restriction[index][2 * (i + 2 * j) + 1](
1177 (i + 4 * j + 2) * this->degree, dof) *
1183 i, j, 2.0 * edge_quadrature_points[q_point][0]);
1186 (2.0 * this->shape_value_component(
1187 dof, quadrature_point_3, 2) -
1188 this->restriction[index][i + 2 * j](
1189 (i + 2 * (j + 4)) * this->degree, dof) *
1190 this->shape_value_component(
1191 (i + 2 * (j + 4)) * this->degree,
1196 this->restriction[index][i + 2 * (j + 2)](
1197 (i + 2 * (j + 4)) * this->degree, dof) *
1198 this->shape_value_component((i + 2 * (j + 4)) *
1209 (i + 4 * j) * this->
degree, dof) *
1217 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1222 dof, quadrature_point_3, 1) -
1223 this->restriction[index][i + 4 * j + 2](
1224 (i + 4 * j) * this->degree, dof) *
1226 (i + 4 * j) * this->degree,
1231 this->restriction[index][2 * (i + 2 * j)](
1232 (i + 4 * j + 2) * this->degree, dof) *
1238 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1243 (2.0 * this->shape_value_component(
1244 dof, quadrature_point_3, 0) -
1245 this->restriction[index][2 * (i + 2 * j) + 1](
1246 (i + 4 * j + 2) * this->degree, dof) *
1247 this->shape_value_component(
1248 (i + 4 * j + 2) * this->degree,
1253 this->restriction[index][i + 2 * j](
1254 (i + 2 * (j + 4)) * this->degree, dof) *
1255 this->shape_value_component((i + 2 * (j + 4)) *
1262 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1265 (2.0 * this->shape_value_component(
1266 dof, quadrature_point_3, 2) -
1267 this->restriction[index][i + 2 * (j + 2)](
1268 (i + 2 * (j + 4)) * this->degree, dof) *
1269 this->shape_value_component(
1270 (i + 2 * (j + 4)) * this->degree,
1275 for (
unsigned int k = 0; k < deg; ++k)
1278 legendre_polynomials[k + 1].value(
1279 edge_quadrature_points[q_point][0]);
1281 for (
unsigned int l = 0; l < tmp.
size(); ++l)
1282 system_rhs(k, l) += tmp(l) * L_k;
1286 system_matrix_inv.
mmult(solution, system_rhs);
1288 for (
unsigned int k = 0; k < 2; ++k)
1289 for (
unsigned int l = 0; l < deg; ++l)
1291 if (
std::abs(solution(l, k)) > 1e-14)
1293 (i + 4 * j) * this->
degree + l + 1, dof) =
1296 if (
std::abs(solution(l, k + 2)) > 1e-14)
1298 (i + 4 * j + 2) * this->
degree + l + 1, dof) =
1301 if (
std::abs(solution(l, k + 4)) > 1e-14)
1303 (i + 2 * (j + 4)) * this->
degree + l + 1, dof) =
1309 const std::vector<Point<2>> &face_quadrature_points =
1310 face_quadrature.get_points();
1311 const std::vector<Polynomials::Polynomial<double>>
1312 &lobatto_polynomials =
1314 const unsigned int n_edge_dofs =
1316 const unsigned int n_face_quadrature_points =
1317 face_quadrature.size();
1321 n_face_quadrature_points);
1323 for (
unsigned int q_point = 0;
1324 q_point < n_face_quadrature_points;
1327 const double weight =
1328 std::sqrt(face_quadrature.weight(q_point));
1330 for (
unsigned int i = 0; i <= deg; ++i)
1333 weight * legendre_polynomials[i].value(
1334 face_quadrature_points[q_point][0]);
1336 for (
unsigned int j = 0; j < deg; ++j)
1337 assembling_matrix(i * deg + j, q_point) =
1338 L_i * lobatto_polynomials[j + 2].value(
1339 face_quadrature_points[q_point][1]);
1344 assembling_matrix.
m());
1346 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1347 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1348 system_matrix_inv.
invert(system_matrix);
1351 solution.reinit(system_matrix_inv.
m(), 24);
1352 system_rhs.reinit(system_matrix_inv.
m(), 24);
1355 for (
unsigned int i = 0; i < 2; ++i)
1356 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1360 for (
unsigned int q_point = 0;
1361 q_point < n_face_quadrature_points;
1366 if (face_quadrature_points[q_point][0] < 0.5)
1368 if (face_quadrature_points[q_point][1] < 0.5)
1372 2.0 * face_quadrature_points[q_point][0],
1373 2.0 * face_quadrature_points[q_point][1]);
1376 dof, quadrature_point_0, 1);
1378 dof, quadrature_point_0, 2);
1380 2.0 * face_quadrature_points[q_point][0],
1382 2.0 * face_quadrature_points[q_point][1]);
1384 dof, quadrature_point_0, 2);
1386 dof, quadrature_point_0, 0);
1388 2.0 * face_quadrature_points[q_point][0],
1389 2.0 * face_quadrature_points[q_point][1],
1392 dof, quadrature_point_0, 0);
1394 dof, quadrature_point_0, 1);
1401 2.0 * face_quadrature_points[q_point][0],
1402 2.0 * face_quadrature_points[q_point][1] -
1406 dof, quadrature_point_0, 1);
1408 dof, quadrature_point_0, 2);
1410 2.0 * face_quadrature_points[q_point][0],
1412 2.0 * face_quadrature_points[q_point][1] -
1415 dof, quadrature_point_0, 2);
1417 dof, quadrature_point_0, 0);
1419 2.0 * face_quadrature_points[q_point][0],
1420 2.0 * face_quadrature_points[q_point][1] -
1424 dof, quadrature_point_0, 0);
1426 dof, quadrature_point_0, 1);
1430 else if (face_quadrature_points[q_point][1] < 0.5)
1434 2.0 * face_quadrature_points[q_point][0] - 1.0,
1435 2.0 * face_quadrature_points[q_point][1]);
1438 dof, quadrature_point_0, 1);
1440 dof, quadrature_point_0, 2);
1442 2.0 * face_quadrature_points[q_point][0] - 1.0,
1444 2.0 * face_quadrature_points[q_point][1]);
1446 dof, quadrature_point_0, 2);
1448 dof, quadrature_point_0, 0);
1450 2.0 * face_quadrature_points[q_point][0] - 1.0,
1451 2.0 * face_quadrature_points[q_point][1],
1454 dof, quadrature_point_0, 0);
1456 dof, quadrature_point_0, 1);
1463 2.0 * face_quadrature_points[q_point][0] - 1.0,
1464 2.0 * face_quadrature_points[q_point][1] - 1.0);
1467 dof, quadrature_point_0, 1);
1469 dof, quadrature_point_0, 2);
1471 2.0 * face_quadrature_points[q_point][0] - 1.0,
1473 2.0 * face_quadrature_points[q_point][1] - 1.0);
1475 dof, quadrature_point_0, 2);
1477 dof, quadrature_point_0, 0);
1479 2.0 * face_quadrature_points[q_point][0] - 1.0,
1480 2.0 * face_quadrature_points[q_point][1] - 1.0,
1483 dof, quadrature_point_0, 0);
1485 dof, quadrature_point_0, 1);
1490 face_quadrature_points[q_point][0],
1491 face_quadrature_points[q_point][1]);
1493 face_quadrature_points[q_point][0],
1495 face_quadrature_points[q_point][1]);
1497 face_quadrature_points[q_point][0],
1498 face_quadrature_points[q_point][1],
1501 for (
unsigned int j = 0; j < 2; ++j)
1502 for (
unsigned int k = 0; k < 2; ++k)
1503 for (
unsigned int l = 0; l <= deg; ++l)
1505 tmp(2 * (j + 2 * k)) -=
1507 (i + 4 * j) * this->degree + l, dof) *
1509 (i + 4 * j) * this->degree + l,
1512 tmp(2 * (j + 2 * k) + 1) -=
1514 (i + 2 * (k + 4)) * this->degree + l, dof) *
1516 (i + 2 * (k + 4)) * this->degree + l,
1519 tmp(2 * (j + 2 * (k + 2))) -=
1521 (2 * (i + 4) + k) * this->degree + l, dof) *
1523 (2 * (i + 4) + k) * this->degree + l,
1526 tmp(2 * (j + 2 * k) + 9) -=
1528 (i + 4 * j + 2) * this->degree + l, dof) *
1530 (i + 4 * j + 2) * this->degree + l,
1533 tmp(2 * (j + 2 * (k + 4))) -=
1535 (4 * i + j + 2) * this->degree + l, dof) *
1537 (4 * i + j + 2) * this->degree + l,
1540 tmp(2 * (j + 2 * k) + 17) -=
1542 (4 * i + k) * this->degree + l, dof) *
1544 (4 * i + k) * this->degree + l,
1549 tmp *= face_quadrature.weight(q_point);
1551 for (
unsigned int j = 0; j <= deg; ++j)
1553 const double L_j_0 = legendre_polynomials[j].value(
1554 face_quadrature_points[q_point][0]);
1555 const double L_j_1 = legendre_polynomials[j].value(
1556 face_quadrature_points[q_point][1]);
1558 for (
unsigned int k = 0; k < deg; ++k)
1560 const double l_k_0 =
1561 L_j_0 * lobatto_polynomials[k + 2].value(
1562 face_quadrature_points[q_point][1]);
1563 const double l_k_1 =
1564 L_j_1 * lobatto_polynomials[k + 2].value(
1565 face_quadrature_points[q_point][0]);
1567 for (
unsigned int l = 0; l < 4; ++l)
1569 system_rhs(j * deg + k, 2 * l) +=
1571 system_rhs(j * deg + k, 2 * l + 1) +=
1572 tmp(2 * l + 1) * l_k_1;
1573 system_rhs(j * deg + k, 2 * (l + 4)) +=
1574 tmp(2 * (l + 4)) * l_k_1;
1575 system_rhs(j * deg + k, 2 * l + 9) +=
1576 tmp(2 * l + 9) * l_k_0;
1577 system_rhs(j * deg + k, 2 * (l + 8)) +=
1578 tmp(2 * (l + 8)) * l_k_0;
1579 system_rhs(j * deg + k, 2 * l + 17) +=
1580 tmp(2 * l + 17) * l_k_1;
1586 system_matrix_inv.
mmult(solution, system_rhs);
1588 for (
unsigned int j = 0; j < 2; ++j)
1589 for (
unsigned int k = 0; k < 2; ++k)
1590 for (
unsigned int l = 0; l <= deg; ++l)
1591 for (
unsigned int m = 0; m < deg; ++m)
1594 2 * (j + 2 * k))) > 1e-14)
1596 (2 * i * this->degree + l) * deg + m +
1598 dof) = solution(l * deg + m, 2 * (j + 2 * k));
1601 2 * (j + 2 * k) + 1)) >
1604 ((2 * i + 1) * deg + m) * this->degree + l +
1607 solution(l * deg + m, 2 * (j + 2 * k) + 1);
1610 2 * (j + 2 * (k + 2)))) >
1613 (2 * (i + 2) * this->degree + l) * deg + m +
1616 solution(l * deg + m, 2 * (j + 2 * (k + 2)));
1619 2 * (j + 2 * k) + 9)) >
1622 ((2 * i + 5) * deg + m) * this->degree + l +
1625 solution(l * deg + m, 2 * (j + 2 * k) + 9);
1628 2 * (j + 2 * (k + 4)))) >
1631 (2 * (i + 4) * this->degree + l) * deg + m +
1634 solution(l * deg + m, 2 * (j + 2 * (k + 4)));
1637 2 * (j + 2 * k) + 17)) >
1640 ((2 * i + 9) * deg + m) * this->degree + l +
1643 solution(l * deg + m, 2 * (j + 2 * k) + 17);
1648 const std::vector<Point<dim>> &quadrature_points =
1649 quadrature.get_points();
1650 const unsigned int n_boundary_dofs =
1653 const unsigned int n_quadrature_points = quadrature.size();
1657 n_quadrature_points);
1659 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1662 const double weight =
std::sqrt(quadrature.weight(q_point));
1664 for (
unsigned int i = 0; i <= deg; ++i)
1667 weight * legendre_polynomials[i].value(
1668 quadrature_points[q_point][0]);
1670 for (
unsigned int j = 0; j < deg; ++j)
1673 L_i * lobatto_polynomials[j + 2].value(
1674 quadrature_points[q_point][1]);
1676 for (
unsigned int k = 0; k < deg; ++k)
1677 assembling_matrix((i * deg + j) * deg + k,
1679 l_j * lobatto_polynomials[k + 2].value(
1680 quadrature_points[q_point][2]);
1686 assembling_matrix.
m());
1688 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1689 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1690 system_matrix_inv.
invert(system_matrix);
1693 solution.reinit(system_matrix_inv.
m(), 24);
1694 system_rhs.reinit(system_matrix_inv.
m(), 24);
1697 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1701 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1706 if (quadrature_points[q_point][0] < 0.5)
1708 if (quadrature_points[q_point][1] < 0.5)
1710 if (quadrature_points[q_point][2] < 0.5)
1713 2.0 * quadrature_points[q_point][0],
1714 2.0 * quadrature_points[q_point][1],
1715 2.0 * quadrature_points[q_point][2]);
1718 dof, quadrature_point, 0);
1720 dof, quadrature_point, 1);
1722 dof, quadrature_point, 2);
1728 2.0 * quadrature_points[q_point][0],
1729 2.0 * quadrature_points[q_point][1],
1730 2.0 * quadrature_points[q_point][2] - 1.0);
1733 dof, quadrature_point, 0);
1735 dof, quadrature_point, 1);
1737 dof, quadrature_point, 2);
1741 else if (quadrature_points[q_point][2] < 0.5)
1744 2.0 * quadrature_points[q_point][0],
1745 2.0 * quadrature_points[q_point][1] - 1.0,
1746 2.0 * quadrature_points[q_point][2]);
1749 dof, quadrature_point, 0);
1751 dof, quadrature_point, 1);
1753 dof, quadrature_point, 2);
1759 2.0 * quadrature_points[q_point][0],
1760 2.0 * quadrature_points[q_point][1] - 1.0,
1761 2.0 * quadrature_points[q_point][2] - 1.0);
1764 dof, quadrature_point, 0);
1766 dof, quadrature_point, 1);
1768 dof, quadrature_point, 2);
1772 else if (quadrature_points[q_point][1] < 0.5)
1774 if (quadrature_points[q_point][2] < 0.5)
1777 2.0 * quadrature_points[q_point][0] - 1.0,
1778 2.0 * quadrature_points[q_point][1],
1779 2.0 * quadrature_points[q_point][2]);
1782 dof, quadrature_point, 0);
1784 dof, quadrature_point, 1);
1786 dof, quadrature_point, 2);
1792 2.0 * quadrature_points[q_point][0] - 1.0,
1793 2.0 * quadrature_points[q_point][1],
1794 2.0 * quadrature_points[q_point][2] - 1.0);
1797 dof, quadrature_point, 0);
1799 dof, quadrature_point, 1);
1801 dof, quadrature_point, 2);
1805 else if (quadrature_points[q_point][2] < 0.5)
1808 2.0 * quadrature_points[q_point][0] - 1.0,
1809 2.0 * quadrature_points[q_point][1] - 1.0,
1810 2.0 * quadrature_points[q_point][2]);
1829 2.0 * quadrature_points[q_point][0] - 1.0,
1830 2.0 * quadrature_points[q_point][1] - 1.0,
1831 2.0 * quadrature_points[q_point][2] - 1.0);
1847 for (
unsigned int i = 0; i < 2; ++i)
1848 for (
unsigned int j = 0; j < 2; ++j)
1849 for (
unsigned int k = 0; k < 2; ++k)
1850 for (
unsigned int l = 0; l <= deg; ++l)
1852 tmp(3 * (i + 2 * (j + 2 * k))) -=
1854 (4 * i + j + 2) * this->degree + l, dof) *
1856 (4 * i + j + 2) * this->degree + l,
1857 quadrature_points[q_point],
1859 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1861 (4 * i + k) * this->degree + l, dof) *
1863 (4 * i + k) * this->degree + l,
1864 quadrature_points[q_point],
1866 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1868 (2 * (j + 4) + k) * this->degree + l, dof) *
1870 (2 * (j + 4) + k) * this->degree + l,
1871 quadrature_points[q_point],
1874 for (
unsigned int m = 0; m < deg; ++m)
1876 tmp(3 * (i + 2 * (j + 2 * k))) -=
1879 ((2 * j + 5) * deg + m) * this->degree +
1883 ((2 * j + 5) * deg + m) * this->degree +
1885 quadrature_points[q_point],
1887 tmp(3 * (i + 2 * (j + 2 * k))) -=
1890 (2 * (i + 4) * this->degree + l) * deg +
1894 (2 * (i + 4) * this->degree + l) * deg +
1896 quadrature_points[q_point],
1898 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1901 (2 * k * this->degree + l) * deg + m +
1905 (2 * k * this->degree + l) * deg + m +
1907 quadrature_points[q_point],
1909 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1912 ((2 * i + 9) * deg + m) * this->degree +
1916 ((2 * i + 9) * deg + m) * this->degree +
1918 quadrature_points[q_point],
1920 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1923 ((2 * k + 1) * deg + m) * this->degree +
1927 ((2 * k + 1) * deg + m) * this->degree +
1929 quadrature_points[q_point],
1931 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1934 (2 * (j + 2) * this->degree + l) * deg +
1938 (2 * (j + 2) * this->degree + l) * deg +
1940 quadrature_points[q_point],
1945 tmp *= quadrature.weight(q_point);
1947 for (
unsigned int i = 0; i <= deg; ++i)
1949 const double L_i_0 = legendre_polynomials[i].value(
1950 quadrature_points[q_point][0]);
1951 const double L_i_1 = legendre_polynomials[i].value(
1952 quadrature_points[q_point][1]);
1953 const double L_i_2 = legendre_polynomials[i].value(
1954 quadrature_points[q_point][2]);
1956 for (
unsigned int j = 0; j < deg; ++j)
1958 const double l_j_0 =
1959 L_i_0 * lobatto_polynomials[j + 2].value(
1960 quadrature_points[q_point][1]);
1961 const double l_j_1 =
1962 L_i_1 * lobatto_polynomials[j + 2].value(
1963 quadrature_points[q_point][0]);
1964 const double l_j_2 =
1965 L_i_2 * lobatto_polynomials[j + 2].value(
1966 quadrature_points[q_point][0]);
1968 for (
unsigned int k = 0; k < deg; ++k)
1970 const double l_k_0 =
1971 l_j_0 * lobatto_polynomials[k + 2].value(
1972 quadrature_points[q_point][2]);
1973 const double l_k_1 =
1974 l_j_1 * lobatto_polynomials[k + 2].value(
1975 quadrature_points[q_point][2]);
1976 const double l_k_2 =
1977 l_j_2 * lobatto_polynomials[k + 2].value(
1978 quadrature_points[q_point][1]);
1980 for (
unsigned int l = 0; l < 8; ++l)
1982 system_rhs((i * deg + j) * deg + k,
1983 3 * l) += tmp(3 * l) * l_k_0;
1984 system_rhs((i * deg + j) * deg + k,
1986 tmp(3 * l + 1) * l_k_1;
1987 system_rhs((i * deg + j) * deg + k,
1989 tmp(3 * l + 2) * l_k_2;
1996 system_matrix_inv.
mmult(solution, system_rhs);
1998 for (
unsigned int i = 0; i < 2; ++i)
1999 for (
unsigned int j = 0; j < 2; ++j)
2000 for (
unsigned int k = 0; k < 2; ++k)
2001 for (
unsigned int l = 0; l <= deg; ++l)
2002 for (
unsigned int m = 0; m < deg; ++m)
2003 for (
unsigned int n = 0; n < deg; ++n)
2006 solution((l * deg + m) * deg + n,
2007 3 * (i + 2 * (j + 2 * k)))) >
2010 (l * deg + m) * deg + n + n_boundary_dofs,
2011 dof) = solution((l * deg + m) * deg + n,
2012 3 * (i + 2 * (j + 2 * k)));
2015 solution((l * deg + m) * deg + n,
2016 3 * (i + 2 * (j + 2 * k)) + 1)) >
2019 (l + (m + deg) * this->degree) * deg + n +
2022 solution((l * deg + m) * deg + n,
2023 3 * (i + 2 * (j + 2 * k)) + 1);
2026 solution((l * deg + m) * deg + n,
2027 3 * (i + 2 * (j + 2 * k)) + 2)) >
2031 ((m + 2 * deg) * deg + n) * this->degree +
2034 solution((l * deg + m) * deg + n,
2035 3 * (i + 2 * (j + 2 * k)) + 2);
3142 std::vector<double> &nodal_values)
const
3147 const unsigned int face_no = 0;
3149 const unsigned int deg = this->
degree - 1;
3150 Assert(support_point_values.size() == this->generalized_support_points.size(),
3152 this->generalized_support_points.size()));
3156 Assert(nodal_values.size() == this->n_dofs_per_cell(),
3158 std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3167 const unsigned int n_edge_points = reference_edge_quadrature.size();
3169 for (
unsigned int i = 0; i < 2; ++i)
3170 for (
unsigned int j = 0; j < 2; ++j)
3172 for (
unsigned int q_point = 0; q_point < n_edge_points;
3174 nodal_values[(i + 2 * j) * this->
degree] +=
3175 reference_edge_quadrature.weight(q_point) *
3176 support_point_values[q_point + (i + 2 * j) * n_edge_points]
3183 nodal_values[(i + 2 * j) * this->degree] = 0.0;
3201 const std::vector<Polynomials::Polynomial<double>>
3202 &lobatto_polynomials =
3206 std::vector<Polynomials::Polynomial<double>>
3207 lobatto_polynomials_grad(this->
degree);
3209 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3210 lobatto_polynomials_grad[i] =
3211 lobatto_polynomials[i + 1].derivative();
3216 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3217 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3218 for (
unsigned int q_point = 0; q_point < n_edge_points;
3220 system_matrix(i, j) +=
3222 lobatto_polynomials_grad[i + 1].value(
3229 system_matrix_inv.
invert(system_matrix);
3236 for (
unsigned int line = 0;
3237 line < GeometryInfo<dim>::lines_per_cell;
3243 for (
unsigned int q_point = 0; q_point < n_edge_points;
3247 support_point_values[line * n_edge_points + q_point]
3248 [line_coordinate[line]] -
3249 nodal_values[line * this->
degree] *
3255 line_coordinate[line]);
3257 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3261 system_matrix_inv.
vmult(solution, system_rhs);
3267 for (
unsigned int i = 0; i < solution.
size(); ++i)
3269 nodal_values[line * this->
degree + i + 1] = solution(i);
3282 const unsigned int n_interior_points =
3283 reference_quadrature.size();
3284 const std::vector<Polynomials::Polynomial<double>>
3285 &legendre_polynomials =
3289 system_matrix.reinit((this->
degree - 1) * this->
degree,
3293 for (
unsigned int i = 0; i < this->
degree; ++i)
3294 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3295 for (
unsigned int k = 0; k < this->degree; ++k)
3296 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3297 for (
unsigned int q_point = 0;
3298 q_point < n_interior_points;
3300 system_matrix(i * (this->degree - 1) + j,
3301 k * (this->degree - 1) + l) +=
3302 reference_quadrature.weight(q_point) *
3303 legendre_polynomials[i].value(
3306 n_edge_points][0]) *
3307 lobatto_polynomials[j + 2].value(
3310 n_edge_points][1]) *
3311 lobatto_polynomials_grad[k].value(
3314 n_edge_points][0]) *
3315 lobatto_polynomials[l + 2].value(
3320 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3321 system_matrix_inv.
invert(system_matrix);
3325 system_rhs.
reinit(system_matrix_inv.
m());
3328 for (
unsigned int q_point = 0; q_point < n_interior_points;
3332 support_point_values[q_point +
3336 for (
unsigned int i = 0; i < 2; ++i)
3337 for (
unsigned int j = 0; j <= deg; ++j)
3338 tmp -= nodal_values[(i + 2) * this->degree + j] *
3340 (i + 2) * this->degree + j,
3346 for (
unsigned int i = 0; i <= deg; ++i)
3347 for (
unsigned int j = 0; j < deg; ++j)
3348 system_rhs(i * deg + j) +=
3349 reference_quadrature.weight(q_point) * tmp *
3350 lobatto_polynomials_grad[i].value(
3353 n_edge_points][0]) *
3354 lobatto_polynomials[j + 2].value(
3360 solution.
reinit(system_matrix.
m());
3361 system_matrix_inv.
vmult(solution, system_rhs);
3367 for (
unsigned int i = 0; i <= deg; ++i)
3368 for (
unsigned int j = 0; j < deg; ++j)
3369 if (
std::abs(solution(i * deg + j)) > 1e-14)
3372 solution(i * deg + j);
3379 for (
unsigned int q_point = 0; q_point < n_interior_points;
3383 support_point_values[q_point +
3387 for (
unsigned int i = 0; i < 2; ++i)
3388 for (
unsigned int j = 0; j <= deg; ++j)
3389 tmp -= nodal_values[i * this->degree + j] *
3391 i * this->degree + j,
3397 for (
unsigned int i = 0; i <= deg; ++i)
3398 for (
unsigned int j = 0; j < deg; ++j)
3399 system_rhs(i * deg + j) +=
3400 reference_quadrature.weight(q_point) * tmp *
3401 lobatto_polynomials_grad[i].value(
3404 n_edge_points][1]) *
3405 lobatto_polynomials[j + 2].value(
3411 system_matrix_inv.
vmult(solution, system_rhs);
3417 for (
unsigned int i = 0; i <= deg; ++i)
3418 for (
unsigned int j = 0; j < deg; ++j)
3419 if (
std::abs(solution(i * deg + j)) > 1e-14)
3422 this->degree] = solution(i * deg + j);
3433 const unsigned int n_edge_points = reference_edge_quadrature.size();
3435 for (
unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3437 for (
unsigned int i = 0; i < 4; ++i)
3438 nodal_values[(i + 8) * this->
degree] +=
3439 reference_edge_quadrature.weight(q_point) *
3440 support_point_values[q_point + (i + 8) * n_edge_points][2];
3442 for (
unsigned int i = 0; i < 2; ++i)
3443 for (
unsigned int j = 0; j < 2; ++j)
3444 for (
unsigned int k = 0; k < 2; ++k)
3445 nodal_values[(i + 2 * (2 * j + k)) * this->
degree] +=
3446 reference_edge_quadrature.weight(q_point) *
3447 support_point_values[q_point + (i + 2 * (2 * j + k)) *
3448 n_edge_points][1 - k];
3455 for (
unsigned int i = 0; i < 4; ++i)
3457 nodal_values[(i + 8) * this->
degree] = 0.0;
3459 for (
unsigned int i = 0; i < 2; ++i)
3460 for (
unsigned int j = 0; j < 2; ++j)
3461 for (
unsigned int k = 0; k < 2; ++k)
3463 nodal_values[(i + 2 * (2 * j + k)) * this->
degree]) <
3465 nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3476 if (this->degree > 1)
3481 const std::vector<Polynomials::Polynomial<double>>
3482 &lobatto_polynomials =
3486 std::vector<Polynomials::Polynomial<double>>
3487 lobatto_polynomials_grad(this->degree);
3489 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3490 lobatto_polynomials_grad[i] =
3491 lobatto_polynomials[i + 1].derivative();
3496 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3497 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3498 for (
unsigned int q_point = 0; q_point < n_edge_points;
3500 system_matrix(i, j) +=
3502 lobatto_polynomials_grad[i + 1].value(
3509 system_matrix_inv.
invert(system_matrix);
3513 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3517 for (
unsigned int line = 0;
3518 line < GeometryInfo<dim>::lines_per_cell;
3524 for (
unsigned int q_point = 0; q_point < this->degree;
3528 support_point_values[line * this->degree + q_point]
3529 [line_coordinate[line]] -
3530 nodal_values[line * this->degree] *
3532 line * this->degree,
3536 line_coordinate[line]);
3538 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3542 system_matrix_inv.
vmult(solution, system_rhs);
3548 for (
unsigned int i = 0; i < solution.
size(); ++i)
3550 nodal_values[line * this->degree + i + 1] = solution(i);
3561 const std::vector<Polynomials::Polynomial<double>>
3562 &legendre_polynomials =
3565 const unsigned int n_face_points = n_edge_points * n_edge_points;
3567 system_matrix.reinit((this->degree - 1) * this->degree,
3568 (this->degree - 1) * this->degree);
3571 for (
unsigned int i = 0; i < this->degree; ++i)
3572 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3573 for (
unsigned int k = 0; k < this->degree; ++k)
3574 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3575 for (
unsigned int q_point = 0; q_point < n_face_points;
3577 system_matrix(i * (this->degree - 1) + j,
3578 k * (this->degree - 1) + l) +=
3580 2 * (k * (this->degree - 1) + l)) *
3581 legendre_polynomials[i].value(
3583 [face_no][q_point + 4 * n_edge_points][0]) *
3584 lobatto_polynomials[j + 2].value(
3586 [face_no][q_point + 4 * n_edge_points][1]);
3588 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3589 system_matrix_inv.
invert(system_matrix);
3590 solution.
reinit(system_matrix.
m());
3591 system_rhs.
reinit(system_matrix.
m());
3595 {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3612 for (
unsigned int q_point = 0; q_point < n_face_points;
3616 support_point_values[q_point +
3619 [face_coordinates[face][0]];
3621 for (
unsigned int i = 0; i < 2; ++i)
3622 for (
unsigned int j = 0; j <= deg; ++j)
3624 nodal_values[edge_indices[face][i] * this->degree +
3627 edge_indices[face][i] * this->degree + j,
3631 face_coordinates[face][0]);
3633 for (
unsigned int i = 0; i <= deg; ++i)
3634 for (
unsigned int j = 0; j < deg; ++j)
3635 system_rhs(i * deg + j) +=
3637 2 * (i * deg + j)) *
3641 system_matrix_inv.
vmult(solution, system_rhs);
3647 for (
unsigned int i = 0; i <= deg; ++i)
3648 for (
unsigned int j = 0; j < deg; ++j)
3649 if (
std::abs(solution(i * deg + j)) > 1e-14)
3650 nodal_values[(2 * face * this->degree + i +
3654 solution(i * deg + j);
3661 for (
unsigned int q_point = 0; q_point < n_face_points;
3665 support_point_values[q_point +
3668 [face_coordinates[face][1]];
3670 for (
unsigned int i = 2;
3671 i < GeometryInfo<dim>::lines_per_face;
3673 for (
unsigned int j = 0; j <= deg; ++j)
3675 nodal_values[edge_indices[face][i] * this->degree +
3678 edge_indices[face][i] * this->degree + j,
3682 face_coordinates[face][1]);
3684 for (
unsigned int i = 0; i <= deg; ++i)
3685 for (
unsigned int j = 0; j < deg; ++j)
3686 system_rhs(i * deg + j) +=
3688 2 * (i * deg + j) + 1) *
3692 system_matrix_inv.
vmult(solution, system_rhs);
3698 for (
unsigned int i = 0; i <= deg; ++i)
3699 for (
unsigned int j = 0; j < deg; ++j)
3700 if (
std::abs(solution(i * deg + j)) > 1e-14)
3701 nodal_values[((2 * face + 1) * deg + j +
3704 i] = solution(i * deg + j);
3712 const QGauss<dim> reference_quadrature(this->degree);
3713 const unsigned int n_interior_points =
3714 reference_quadrature.size();
3718 system_matrix.reinit(this->degree * deg * deg,
3719 this->degree * deg * deg);
3722 for (
unsigned int i = 0; i <= deg; ++i)
3723 for (
unsigned int j = 0; j < deg; ++j)
3724 for (
unsigned int k = 0; k < deg; ++k)
3725 for (
unsigned int l = 0; l <= deg; ++l)
3726 for (
unsigned int m = 0; m < deg; ++m)
3727 for (
unsigned int n = 0; n < deg; ++n)
3728 for (
unsigned int q_point = 0;
3729 q_point < n_interior_points;
3731 system_matrix((i * deg + j) * deg + k,
3732 (l * deg + m) * deg + n) +=
3733 reference_quadrature.weight(q_point) *
3734 legendre_polynomials[i].value(
3740 n_face_points][0]) *
3741 lobatto_polynomials[j + 2].value(
3747 n_face_points][1]) *
3748 lobatto_polynomials[k + 2].value(
3754 n_face_points][2]) *
3755 lobatto_polynomials_grad[l].value(
3761 n_face_points][0]) *
3762 lobatto_polynomials[m + 2].value(
3768 n_face_points][1]) *
3769 lobatto_polynomials[n + 2].value(
3777 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3778 system_matrix_inv.
invert(system_matrix);
3780 system_rhs.
reinit(system_matrix.
m());
3783 for (
unsigned int q_point = 0; q_point < n_interior_points;
3787 support_point_values[q_point +
3793 for (
unsigned int i = 0; i <= deg; ++i)
3795 for (
unsigned int j = 0; j < 2; ++j)
3796 for (
unsigned int k = 0; k < 2; ++k)
3798 nodal_values[i + (j + 4 * k + 2) * this->degree] *
3800 i + (j + 4 * k + 2) * this->degree,
3809 for (
unsigned int j = 0; j < deg; ++j)
3810 for (
unsigned int k = 0; k < 4; ++k)
3812 nodal_values[(i + 2 * (k + 2) * this->degree +
3817 this->shape_value_component(
3818 (i + 2 * (k + 2) * this->degree +
3831 for (
unsigned int i = 0; i <= deg; ++i)
3832 for (
unsigned int j = 0; j < deg; ++j)
3833 for (
unsigned int k = 0; k < deg; ++k)
3834 system_rhs((i * deg + j) * deg + k) +=
3835 reference_quadrature.weight(q_point) * tmp *
3836 lobatto_polynomials_grad[i].value(
3842 n_face_points][0]) *
3843 lobatto_polynomials[j + 2].value(
3849 n_face_points][1]) *
3850 lobatto_polynomials[k + 2].value(
3860 system_matrix_inv.
vmult(solution, system_rhs);
3866 for (
unsigned int i = 0; i <= deg; ++i)
3867 for (
unsigned int j = 0; j < deg; ++j)
3868 for (
unsigned int k = 0; k < deg; ++k)
3869 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3876 solution((i * deg + j) * deg + k);
3881 for (
unsigned int q_point = 0; q_point < n_interior_points;
3885 support_point_values[q_point +
3891 for (
unsigned int i = 0; i <= deg; ++i)
3892 for (
unsigned int j = 0; j < 2; ++j)
3894 for (
unsigned int k = 0; k < 2; ++k)
3895 tmp -= nodal_values[i + (4 * j + k) * this->degree] *
3897 i + (4 * j + k) * this->degree,
3906 for (
unsigned int k = 0; k < deg; ++k)
3908 nodal_values[(i + 2 * j * this->degree +
3913 this->shape_value_component(
3914 (i + 2 * j * this->degree +
3926 ((2 * j + 9) * deg + k +
3929 this->shape_value_component(
3930 i + ((2 * j + 9) * deg + k +
3942 for (
unsigned int i = 0; i <= deg; ++i)
3943 for (
unsigned int j = 0; j < deg; ++j)
3944 for (
unsigned int k = 0; k < deg; ++k)
3945 system_rhs((i * deg + j) * deg + k) +=
3946 reference_quadrature.weight(q_point) * tmp *
3947 lobatto_polynomials_grad[i].value(
3953 n_face_points][1]) *
3954 lobatto_polynomials[j + 2].value(
3960 n_face_points][0]) *
3961 lobatto_polynomials[k + 2].value(
3970 system_matrix_inv.
vmult(solution, system_rhs);
3976 for (
unsigned int i = 0; i <= deg; ++i)
3977 for (
unsigned int j = 0; j < deg; ++j)
3978 for (
unsigned int k = 0; k < deg; ++k)
3979 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3980 nodal_values[((i + this->degree +
3987 solution((i * deg + j) * deg + k);
3992 for (
unsigned int q_point = 0; q_point < n_interior_points;
3996 support_point_values[q_point +
4002 for (
unsigned int i = 0; i <= deg; ++i)
4003 for (
unsigned int j = 0; j < 4; ++j)
4005 tmp -= nodal_values[i + (j + 8) * this->degree] *
4007 i + (j + 8) * this->degree,
4016 for (
unsigned int k = 0; k < deg; ++k)
4019 ((2 * j + 1) * deg + k +
4023 i + ((2 * j + 1) * deg + k +
4035 for (
unsigned int i = 0; i <= deg; ++i)
4036 for (
unsigned int j = 0; j < deg; ++j)
4037 for (
unsigned int k = 0; k < deg; ++k)
4038 system_rhs((i * deg + j) * deg + k) +=
4039 reference_quadrature.weight(q_point) * tmp *
4040 lobatto_polynomials_grad[i].value(
4046 n_face_points][2]) *
4047 lobatto_polynomials[j + 2].value(
4053 n_face_points][0]) *
4054 lobatto_polynomials[k + 2].value(
4063 system_matrix_inv.
vmult(solution, system_rhs);
4069 for (
unsigned int i = 0; i <= deg; ++i)
4070 for (
unsigned int j = 0; j < deg; ++j)
4071 for (
unsigned int k = 0; k < deg; ++k)
4072 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4078 this->degree] = solution((i * deg + j) * deg + k);