583 const QGauss<1> edge_quadrature(2 * this->degree);
584 const std::vector<Point<1>> &edge_quadrature_points =
586 const unsigned int n_edge_quadrature_points = edge_quadrature.
size();
597 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
598 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
601 const double weight = 2.0 * edge_quadrature.
weight(q_point);
603 if (edge_quadrature_points[q_point](0) < 0.5)
606 0.0, 2.0 * edge_quadrature_points[q_point](0));
608 this->restriction[index][0](0, dof) +=
610 this->shape_value_component(dof, quadrature_point, 1);
611 quadrature_point(0) = 1.0;
612 this->restriction[index][1](this->degree, dof) +=
614 this->shape_value_component(dof, quadrature_point, 1);
615 quadrature_point(0) = quadrature_point(1);
616 quadrature_point(1) = 0.0;
617 this->restriction[index][0](2 * this->degree, dof) +=
619 this->shape_value_component(dof, quadrature_point, 0);
620 quadrature_point(1) = 1.0;
621 this->restriction[index][2](3 * this->degree, dof) +=
623 this->shape_value_component(dof, quadrature_point, 0);
629 0.0, 2.0 * edge_quadrature_points[q_point](0) - 1.0);
631 this->restriction[index][2](0, dof) +=
633 this->shape_value_component(dof, quadrature_point, 1);
634 quadrature_point(0) = 1.0;
635 this->restriction[index][3](this->degree, dof) +=
637 this->shape_value_component(dof, quadrature_point, 1);
638 quadrature_point(0) = quadrature_point(1);
639 quadrature_point(1) = 0.0;
640 this->restriction[index][1](2 * this->degree, dof) +=
642 this->shape_value_component(dof, quadrature_point, 0);
643 quadrature_point(1) = 1.0;
644 this->restriction[index][3](3 * this->degree, dof) +=
646 this->shape_value_component(dof, quadrature_point, 0);
654 if (this->degree > 1)
656 const unsigned int deg = this->degree - 1;
657 const std::vector<Polynomials::Polynomial<double>>
658 &legendre_polynomials =
664 n_edge_quadrature_points);
666 for (
unsigned int q_point = 0;
667 q_point < n_edge_quadrature_points;
670 const double weight =
673 for (
unsigned int i = 0; i < deg; ++i)
674 assembling_matrix(i, q_point) =
675 weight * legendre_polynomials[i + 1].value(
676 edge_quadrature_points[q_point](0));
681 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
682 system_matrix_inv.
invert(system_matrix);
689 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
690 for (
unsigned int i = 0; i < 2; ++i)
694 for (
unsigned int q_point = 0;
695 q_point < n_edge_quadrature_points;
698 const double weight = edge_quadrature.
weight(q_point);
700 i, edge_quadrature_points[q_point](0));
702 edge_quadrature_points[q_point](0), i);
704 if (edge_quadrature_points[q_point](0) < 0.5)
707 i, 2.0 * edge_quadrature_points[q_point](0));
711 (2.0 * this->shape_value_component(
712 dof, quadrature_point_2, 1) -
713 this->restriction[index][i](i * this->degree,
715 this->shape_value_component(i * this->degree,
720 this->restriction[index][i + 2](i * this->degree,
722 this->shape_value_component(i * this->degree,
726 2.0 * edge_quadrature_points[q_point](0), i);
729 (2.0 * this->shape_value_component(
730 dof, quadrature_point_2, 0) -
731 this->restriction[index][2 * i]((i + 2) *
734 this->shape_value_component((i + 2) *
740 this->restriction[index][2 * i + 1](
741 (i + 2) * this->degree, dof) *
742 this->shape_value_component(
743 (i + 2) * this->degree, quadrature_point_1, 0);
750 this->restriction[index][i](i * this->degree,
752 this->shape_value_component(i * this->degree,
758 2.0 * edge_quadrature_points[q_point](0) - 1.0);
762 (2.0 * this->shape_value_component(
763 dof, quadrature_point_2, 1) -
764 this->restriction[index][i + 2](i * this->degree,
766 this->shape_value_component(i * this->degree,
771 this->restriction[index][2 * i]((i + 2) *
774 this->shape_value_component(
775 (i + 2) * this->degree, quadrature_point_1, 0);
777 2.0 * edge_quadrature_points[q_point](0) - 1.0,
781 (2.0 * this->shape_value_component(
782 dof, quadrature_point_2, 0) -
783 this->restriction[index][2 * i + 1](
784 (i + 2) * this->degree, dof) *
785 this->shape_value_component((i + 2) *
791 for (
unsigned int j = 0; j < this->degree - 1; ++j)
794 legendre_polynomials[j + 1].value(
795 edge_quadrature_points[q_point](0));
797 for (
unsigned int k = 0; k < tmp.
size(); ++k)
798 system_rhs(j, k) += tmp(k) * L_j;
802 system_matrix_inv.
mmult(solution, system_rhs);
804 for (
unsigned int j = 0; j < this->degree - 1; ++j)
805 for (
unsigned int k = 0; k < 2; ++k)
807 if (
std::abs(solution(j, k)) > 1e-14)
808 this->restriction[index][i + 2 * k](
809 i * this->degree + j + 1, dof) = solution(j, k);
811 if (
std::abs(solution(j, k + 2)) > 1e-14)
812 this->restriction[index][2 * i + k](
813 (i + 2) * this->degree + j + 1, dof) =
819 const std::vector<Point<dim>> &quadrature_points =
821 const std::vector<Polynomials::Polynomial<double>>
822 &lobatto_polynomials =
824 const unsigned int n_boundary_dofs =
826 const unsigned int n_quadrature_points = quadrature.
size();
831 n_quadrature_points);
833 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
838 for (
unsigned int i = 0; i < this->degree; ++i)
841 weight * legendre_polynomials[i].value(
842 quadrature_points[q_point](0));
844 for (
unsigned int j = 0; j < this->degree - 1; ++j)
845 assembling_matrix(i * (this->degree - 1) + j,
847 L_i * lobatto_polynomials[j + 2].value(
848 quadrature_points[q_point](1));
853 assembling_matrix.
m());
855 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
856 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
857 system_matrix_inv.
invert(system_matrix);
860 solution.reinit(system_matrix_inv.
m(), 8);
861 system_rhs.reinit(system_matrix_inv.
m(), 8);
864 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
868 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
873 if (quadrature_points[q_point](0) < 0.5)
875 if (quadrature_points[q_point](1) < 0.5)
878 2.0 * quadrature_points[q_point](0),
879 2.0 * quadrature_points[q_point](1));
881 tmp(0) += 2.0 * this->shape_value_component(
882 dof, quadrature_point, 0);
883 tmp(1) += 2.0 * this->shape_value_component(
884 dof, quadrature_point, 1);
890 2.0 * quadrature_points[q_point](0),
891 2.0 * quadrature_points[q_point](1) - 1.0);
893 tmp(4) += 2.0 * this->shape_value_component(
894 dof, quadrature_point, 0);
895 tmp(5) += 2.0 * this->shape_value_component(
896 dof, quadrature_point, 1);
900 else if (quadrature_points[q_point](1) < 0.5)
903 2.0 * quadrature_points[q_point](0) - 1.0,
904 2.0 * quadrature_points[q_point](1));
907 2.0 * this->shape_value_component(dof,
911 2.0 * this->shape_value_component(dof,
919 2.0 * quadrature_points[q_point](0) - 1.0,
920 2.0 * quadrature_points[q_point](1) - 1.0);
923 2.0 * this->shape_value_component(dof,
927 2.0 * this->shape_value_component(dof,
932 for (
unsigned int i = 0; i < 2; ++i)
933 for (
unsigned int j = 0; j < this->degree; ++j)
936 this->restriction[index][i](j + 2 * this->degree,
938 this->shape_value_component(
939 j + 2 * this->degree,
940 quadrature_points[q_point],
943 this->restriction[index][i](i * this->degree + j,
945 this->shape_value_component(
946 i * this->degree + j,
947 quadrature_points[q_point],
949 tmp(2 * (i + 2)) -= this->restriction[index][i + 2](
950 j + 3 * this->degree, dof) *
951 this->shape_value_component(
952 j + 3 * this->degree,
953 quadrature_points[q_point],
955 tmp(2 * i + 5) -= this->restriction[index][i + 2](
956 i * this->degree + j, dof) *
957 this->shape_value_component(
958 i * this->degree + j,
959 quadrature_points[q_point],
963 tmp *= quadrature.
weight(q_point);
965 for (
unsigned int i = 0; i < this->degree; ++i)
967 const double L_i_0 = legendre_polynomials[i].value(
968 quadrature_points[q_point](0));
969 const double L_i_1 = legendre_polynomials[i].value(
970 quadrature_points[q_point](1));
972 for (
unsigned int j = 0; j < this->degree - 1; ++j)
975 L_i_0 * lobatto_polynomials[j + 2].value(
976 quadrature_points[q_point](1));
978 L_i_1 * lobatto_polynomials[j + 2].value(
979 quadrature_points[q_point](0));
981 for (
unsigned int k = 0; k < 4; ++k)
983 system_rhs(i * (this->degree - 1) + j,
984 2 * k) += tmp(2 * k) * l_j_0;
985 system_rhs(i * (this->degree - 1) + j,
987 tmp(2 * k + 1) * l_j_1;
993 system_matrix_inv.
mmult(solution, system_rhs);
995 for (
unsigned int i = 0; i < this->degree; ++i)
996 for (
unsigned int j = 0; j < this->degree - 1; ++j)
997 for (
unsigned int k = 0; k < 4; ++k)
999 if (
std::abs(solution(i * (this->degree - 1) + j,
1001 this->restriction[index][k](i * (this->degree - 1) +
1002 j + n_boundary_dofs,
1004 solution(i * (this->degree - 1) + j, 2 * k);
1006 if (
std::abs(solution(i * (this->degree - 1) + j,
1007 2 * k + 1)) > 1e-14)
1008 this->restriction[index][k](
1009 i + (this->degree - 1 + j) * this->degree +
1012 solution(i * (this->degree - 1) + j, 2 * k + 1);
1026 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1027 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
1030 const double weight = 2.0 * edge_quadrature.
weight(q_point);
1032 if (edge_quadrature_points[q_point](0) < 0.5)
1033 for (
unsigned int i = 0; i < 2; ++i)
1034 for (
unsigned int j = 0; j < 2; ++j)
1037 i, 2.0 * edge_quadrature_points[q_point](0), j);
1039 this->restriction[index][i + 4 * j]((i + 4 * j) *
1043 this->shape_value_component(dof, quadrature_point, 1);
1045 Point<dim>(2.0 * edge_quadrature_points[q_point](0),
1048 this->restriction[index][2 * (i + 2 * j)](
1049 (i + 4 * j + 2) * this->degree, dof) +=
1051 this->shape_value_component(dof, quadrature_point, 0);
1055 2.0 * edge_quadrature_points[q_point](0));
1056 this->restriction[index][i + 2 * j]((i + 2 * (j + 4)) *
1060 this->shape_value_component(dof, quadrature_point, 2);
1064 for (
unsigned int i = 0; i < 2; ++i)
1065 for (
unsigned int j = 0; j < 2; ++j)
1068 i, 2.0 * edge_quadrature_points[q_point](0) - 1.0, j);
1070 this->restriction[index][i + 4 * j + 2]((i + 4 * j) *
1074 this->shape_value_component(dof, quadrature_point, 1);
1076 2.0 * edge_quadrature_points[q_point](0) - 1.0, i, j);
1077 this->restriction[index][2 * (i + 2 * j) + 1](
1078 (i + 4 * j + 2) * this->degree, dof) +=
1080 this->shape_value_component(dof, quadrature_point, 0);
1082 i, j, 2.0 * edge_quadrature_points[q_point](0) - 1.0);
1083 this->restriction[index][i + 2 * (j + 2)](
1084 (i + 2 * (j + 4)) * this->degree, dof) +=
1086 this->shape_value_component(dof, quadrature_point, 2);
1094 if (this->degree > 1)
1096 const unsigned int deg = this->degree - 1;
1097 const std::vector<Polynomials::Polynomial<double>>
1098 &legendre_polynomials =
1104 n_edge_quadrature_points);
1106 for (
unsigned int q_point = 0;
1107 q_point < n_edge_quadrature_points;
1110 const double weight =
1113 for (
unsigned int i = 0; i < deg; ++i)
1114 assembling_matrix(i, q_point) =
1115 weight * legendre_polynomials[i + 1].value(
1116 edge_quadrature_points[q_point](0));
1121 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1122 system_matrix_inv.
invert(system_matrix);
1129 for (
unsigned int i = 0; i < 2; ++i)
1130 for (
unsigned int j = 0; j < 2; ++j)
1131 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell();
1136 for (
unsigned int q_point = 0;
1137 q_point < n_edge_quadrature_points;
1140 const double weight = edge_quadrature.
weight(q_point);
1142 i, edge_quadrature_points[q_point](0), j);
1144 edge_quadrature_points[q_point](0), i, j);
1146 i, j, edge_quadrature_points[q_point](0));
1148 if (edge_quadrature_points[q_point](0) < 0.5)
1151 i, 2.0 * edge_quadrature_points[q_point](0), j);
1154 weight * (2.0 * this->shape_value_component(
1155 dof, quadrature_point_3, 1) -
1156 this->restriction[index][i + 4 * j](
1157 (i + 4 * j) * this->degree, dof) *
1158 this->shape_value_component(
1159 (i + 4 * j) * this->degree,
1164 this->restriction[index][i + 4 * j + 2](
1165 (i + 4 * j) * this->degree, dof) *
1166 this->shape_value_component((i + 4 * j) *
1171 2.0 * edge_quadrature_points[q_point](0), i, j);
1174 (2.0 * this->shape_value_component(
1175 dof, quadrature_point_3, 0) -
1176 this->restriction[index][2 * (i + 2 * j)](
1177 (i + 4 * j + 2) * this->degree, dof) *
1178 this->shape_value_component(
1179 (i + 4 * j + 2) * this->degree,
1184 this->restriction[index][2 * (i + 2 * j) + 1](
1185 (i + 4 * j + 2) * this->degree, dof) *
1186 this->shape_value_component((i + 4 * j + 2) *
1191 i, j, 2.0 * edge_quadrature_points[q_point](0));
1194 (2.0 * this->shape_value_component(
1195 dof, quadrature_point_3, 2) -
1196 this->restriction[index][i + 2 * j](
1197 (i + 2 * (j + 4)) * this->degree, dof) *
1198 this->shape_value_component(
1199 (i + 2 * (j + 4)) * this->degree,
1204 this->restriction[index][i + 2 * (j + 2)](
1205 (i + 2 * (j + 4)) * this->degree, dof) *
1206 this->shape_value_component((i + 2 * (j + 4)) *
1216 this->restriction[index][i + 4 * j](
1217 (i + 4 * j) * this->degree, dof) *
1218 this->shape_value_component((i + 4 * j) *
1225 2.0 * edge_quadrature_points[q_point](0) - 1.0,
1229 (2.0 * this->shape_value_component(
1230 dof, quadrature_point_3, 1) -
1231 this->restriction[index][i + 4 * j + 2](
1232 (i + 4 * j) * this->degree, dof) *
1233 this->shape_value_component(
1234 (i + 4 * j) * this->degree,
1239 this->restriction[index][2 * (i + 2 * j)](
1240 (i + 4 * j + 2) * this->degree, dof) *
1241 this->shape_value_component((i + 4 * j + 2) *
1246 2.0 * edge_quadrature_points[q_point](0) - 1.0,
1251 (2.0 * this->shape_value_component(
1252 dof, quadrature_point_3, 0) -
1253 this->restriction[index][2 * (i + 2 * j) + 1](
1254 (i + 4 * j + 2) * this->degree, dof) *
1255 this->shape_value_component(
1256 (i + 4 * j + 2) * this->degree,
1261 this->restriction[index][i + 2 * j](
1262 (i + 2 * (j + 4)) * this->degree, dof) *
1263 this->shape_value_component((i + 2 * (j + 4)) *
1270 2.0 * edge_quadrature_points[q_point](0) - 1.0);
1273 (2.0 * this->shape_value_component(
1274 dof, quadrature_point_3, 2) -
1275 this->restriction[index][i + 2 * (j + 2)](
1276 (i + 2 * (j + 4)) * this->degree, dof) *
1277 this->shape_value_component(
1278 (i + 2 * (j + 4)) * this->degree,
1283 for (
unsigned int k = 0; k < deg; ++k)
1286 legendre_polynomials[k + 1].value(
1287 edge_quadrature_points[q_point](0));
1289 for (
unsigned int l = 0; l < tmp.
size(); ++l)
1290 system_rhs(k, l) += tmp(l) * L_k;
1294 system_matrix_inv.
mmult(solution, system_rhs);
1296 for (
unsigned int k = 0; k < 2; ++k)
1297 for (
unsigned int l = 0; l < deg; ++l)
1299 if (
std::abs(solution(l, k)) > 1e-14)
1300 this->restriction[index][i + 2 * (2 * j + k)](
1301 (i + 4 * j) * this->degree + l + 1, dof) =
1304 if (
std::abs(solution(l, k + 2)) > 1e-14)
1305 this->restriction[index][2 * (i + 2 * j) + k](
1306 (i + 4 * j + 2) * this->degree + l + 1, dof) =
1309 if (
std::abs(solution(l, k + 4)) > 1e-14)
1310 this->restriction[index][i + 2 * (j + 2 * k)](
1311 (i + 2 * (j + 4)) * this->degree + l + 1, dof) =
1316 const QGauss<2> face_quadrature(2 * this->degree);
1317 const std::vector<Point<2>> &face_quadrature_points =
1319 const std::vector<Polynomials::Polynomial<double>>
1320 &lobatto_polynomials =
1322 const unsigned int n_edge_dofs =
1324 const unsigned int n_face_quadrature_points =
1325 face_quadrature.
size();
1329 n_face_quadrature_points);
1331 for (
unsigned int q_point = 0;
1332 q_point < n_face_quadrature_points;
1335 const double weight =
1338 for (
unsigned int i = 0; i <= deg; ++i)
1341 weight * legendre_polynomials[i].value(
1342 face_quadrature_points[q_point](0));
1344 for (
unsigned int j = 0; j < deg; ++j)
1345 assembling_matrix(i * deg + j, q_point) =
1346 L_i * lobatto_polynomials[j + 2].value(
1347 face_quadrature_points[q_point](1));
1352 assembling_matrix.
m());
1354 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1355 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1356 system_matrix_inv.
invert(system_matrix);
1359 solution.reinit(system_matrix_inv.
m(), 24);
1360 system_rhs.reinit(system_matrix_inv.
m(), 24);
1363 for (
unsigned int i = 0; i < 2; ++i)
1364 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1368 for (
unsigned int q_point = 0;
1369 q_point < n_face_quadrature_points;
1374 if (face_quadrature_points[q_point](0) < 0.5)
1376 if (face_quadrature_points[q_point](1) < 0.5)
1380 2.0 * face_quadrature_points[q_point](0),
1381 2.0 * face_quadrature_points[q_point](1));
1383 tmp(0) += 2.0 * this->shape_value_component(
1384 dof, quadrature_point_0, 1);
1385 tmp(1) += 2.0 * this->shape_value_component(
1386 dof, quadrature_point_0, 2);
1388 2.0 * face_quadrature_points[q_point](0),
1390 2.0 * face_quadrature_points[q_point](1));
1391 tmp(8) += 2.0 * this->shape_value_component(
1392 dof, quadrature_point_0, 2);
1393 tmp(9) += 2.0 * this->shape_value_component(
1394 dof, quadrature_point_0, 0);
1396 2.0 * face_quadrature_points[q_point](0),
1397 2.0 * face_quadrature_points[q_point](1),
1399 tmp(16) += 2.0 * this->shape_value_component(
1400 dof, quadrature_point_0, 0);
1401 tmp(17) += 2.0 * this->shape_value_component(
1402 dof, quadrature_point_0, 1);
1409 2.0 * face_quadrature_points[q_point](0),
1410 2.0 * face_quadrature_points[q_point](1) -
1413 tmp(2) += 2.0 * this->shape_value_component(
1414 dof, quadrature_point_0, 1);
1415 tmp(3) += 2.0 * this->shape_value_component(
1416 dof, quadrature_point_0, 2);
1418 2.0 * face_quadrature_points[q_point](0),
1420 2.0 * face_quadrature_points[q_point](1) -
1422 tmp(10) += 2.0 * this->shape_value_component(
1423 dof, quadrature_point_0, 2);
1424 tmp(11) += 2.0 * this->shape_value_component(
1425 dof, quadrature_point_0, 0);
1427 2.0 * face_quadrature_points[q_point](0),
1428 2.0 * face_quadrature_points[q_point](1) -
1431 tmp(18) += 2.0 * this->shape_value_component(
1432 dof, quadrature_point_0, 0);
1433 tmp(19) += 2.0 * this->shape_value_component(
1434 dof, quadrature_point_0, 1);
1438 else if (face_quadrature_points[q_point](1) < 0.5)
1442 2.0 * face_quadrature_points[q_point](0) - 1.0,
1443 2.0 * face_quadrature_points[q_point](1));
1445 tmp(4) += 2.0 * this->shape_value_component(
1446 dof, quadrature_point_0, 1);
1447 tmp(5) += 2.0 * this->shape_value_component(
1448 dof, quadrature_point_0, 2);
1450 2.0 * face_quadrature_points[q_point](0) - 1.0,
1452 2.0 * face_quadrature_points[q_point](1));
1453 tmp(12) += 2.0 * this->shape_value_component(
1454 dof, quadrature_point_0, 2);
1455 tmp(13) += 2.0 * this->shape_value_component(
1456 dof, quadrature_point_0, 0);
1458 2.0 * face_quadrature_points[q_point](0) - 1.0,
1459 2.0 * face_quadrature_points[q_point](1),
1461 tmp(20) += 2.0 * this->shape_value_component(
1462 dof, quadrature_point_0, 0);
1463 tmp(21) += 2.0 * this->shape_value_component(
1464 dof, quadrature_point_0, 1);
1471 2.0 * face_quadrature_points[q_point](0) - 1.0,
1472 2.0 * face_quadrature_points[q_point](1) - 1.0);
1474 tmp(6) += 2.0 * this->shape_value_component(
1475 dof, quadrature_point_0, 1);
1476 tmp(7) += 2.0 * this->shape_value_component(
1477 dof, quadrature_point_0, 2);
1479 2.0 * face_quadrature_points[q_point](0) - 1.0,
1481 2.0 * face_quadrature_points[q_point](1) - 1.0);
1482 tmp(14) += 2.0 * this->shape_value_component(
1483 dof, quadrature_point_0, 2);
1484 tmp(15) += 2.0 * this->shape_value_component(
1485 dof, quadrature_point_0, 0);
1487 2.0 * face_quadrature_points[q_point](0) - 1.0,
1488 2.0 * face_quadrature_points[q_point](1) - 1.0,
1490 tmp(22) += 2.0 * this->shape_value_component(
1491 dof, quadrature_point_0, 0);
1492 tmp(23) += 2.0 * this->shape_value_component(
1493 dof, quadrature_point_0, 1);
1498 face_quadrature_points[q_point](0),
1499 face_quadrature_points[q_point](1));
1501 face_quadrature_points[q_point](0),
1503 face_quadrature_points[q_point](1));
1505 face_quadrature_points[q_point](0),
1506 face_quadrature_points[q_point](1),
1509 for (
unsigned int j = 0; j < 2; ++j)
1510 for (
unsigned int k = 0; k < 2; ++k)
1511 for (
unsigned int l = 0; l <= deg; ++l)
1513 tmp(2 * (j + 2 * k)) -=
1514 this->restriction[index][i + 2 * (2 * j + k)](
1515 (i + 4 * j) * this->degree + l, dof) *
1516 this->shape_value_component(
1517 (i + 4 * j) * this->degree + l,
1520 tmp(2 * (j + 2 * k) + 1) -=
1521 this->restriction[index][i + 2 * (2 * j + k)](
1522 (i + 2 * (k + 4)) * this->degree + l, dof) *
1523 this->shape_value_component(
1524 (i + 2 * (k + 4)) * this->degree + l,
1527 tmp(2 * (j + 2 * (k + 2))) -=
1528 this->restriction[index][2 * (i + 2 * j) + k](
1529 (2 * (i + 4) + k) * this->degree + l, dof) *
1530 this->shape_value_component(
1531 (2 * (i + 4) + k) * this->degree + l,
1534 tmp(2 * (j + 2 * k) + 9) -=
1535 this->restriction[index][2 * (i + 2 * j) + k](
1536 (i + 4 * j + 2) * this->degree + l, dof) *
1537 this->shape_value_component(
1538 (i + 4 * j + 2) * this->degree + l,
1541 tmp(2 * (j + 2 * (k + 4))) -=
1542 this->restriction[index][2 * (2 * i + j) + k](
1543 (4 * i + j + 2) * this->degree + l, dof) *
1544 this->shape_value_component(
1545 (4 * i + j + 2) * this->degree + l,
1548 tmp(2 * (j + 2 * k) + 17) -=
1549 this->restriction[index][2 * (2 * i + j) + k](
1550 (4 * i + k) * this->degree + l, dof) *
1551 this->shape_value_component(
1552 (4 * i + k) * this->degree + l,
1557 tmp *= face_quadrature.
weight(q_point);
1559 for (
unsigned int j = 0; j <= deg; ++j)
1561 const double L_j_0 = legendre_polynomials[j].value(
1562 face_quadrature_points[q_point](0));
1563 const double L_j_1 = legendre_polynomials[j].value(
1564 face_quadrature_points[q_point](1));
1566 for (
unsigned int k = 0; k < deg; ++k)
1568 const double l_k_0 =
1569 L_j_0 * lobatto_polynomials[k + 2].value(
1570 face_quadrature_points[q_point](1));
1571 const double l_k_1 =
1572 L_j_1 * lobatto_polynomials[k + 2].value(
1573 face_quadrature_points[q_point](0));
1575 for (
unsigned int l = 0; l < 4; ++l)
1577 system_rhs(j * deg + k, 2 * l) +=
1579 system_rhs(j * deg + k, 2 * l + 1) +=
1580 tmp(2 * l + 1) * l_k_1;
1581 system_rhs(j * deg + k, 2 * (l + 4)) +=
1582 tmp(2 * (l + 4)) * l_k_1;
1583 system_rhs(j * deg + k, 2 * l + 9) +=
1584 tmp(2 * l + 9) * l_k_0;
1585 system_rhs(j * deg + k, 2 * (l + 8)) +=
1586 tmp(2 * (l + 8)) * l_k_0;
1587 system_rhs(j * deg + k, 2 * l + 17) +=
1588 tmp(2 * l + 17) * l_k_1;
1594 system_matrix_inv.
mmult(solution, system_rhs);
1596 for (
unsigned int j = 0; j < 2; ++j)
1597 for (
unsigned int k = 0; k < 2; ++k)
1598 for (
unsigned int l = 0; l <= deg; ++l)
1599 for (
unsigned int m = 0; m < deg; ++m)
1602 2 * (j + 2 * k))) > 1e-14)
1603 this->restriction[index][i + 2 * (2 * j + k)](
1604 (2 * i * this->degree + l) * deg + m +
1606 dof) = solution(l * deg + m, 2 * (j + 2 * k));
1609 2 * (j + 2 * k) + 1)) >
1611 this->restriction[index][i + 2 * (2 * j + k)](
1612 ((2 * i + 1) * deg + m) * this->degree + l +
1615 solution(l * deg + m, 2 * (j + 2 * k) + 1);
1618 2 * (j + 2 * (k + 2)))) >
1620 this->restriction[index][2 * (i + 2 * j) + k](
1621 (2 * (i + 2) * this->degree + l) * deg + m +
1624 solution(l * deg + m, 2 * (j + 2 * (k + 2)));
1627 2 * (j + 2 * k) + 9)) >
1629 this->restriction[index][2 * (i + 2 * j) + k](
1630 ((2 * i + 5) * deg + m) * this->degree + l +
1633 solution(l * deg + m, 2 * (j + 2 * k) + 9);
1636 2 * (j + 2 * (k + 4)))) >
1638 this->restriction[index][2 * (2 * i + j) + k](
1639 (2 * (i + 4) * this->degree + l) * deg + m +
1642 solution(l * deg + m, 2 * (j + 2 * (k + 4)));
1645 2 * (j + 2 * k) + 17)) >
1647 this->restriction[index][2 * (2 * i + j) + k](
1648 ((2 * i + 9) * deg + m) * this->degree + l +
1651 solution(l * deg + m, 2 * (j + 2 * k) + 17);
1656 const std::vector<Point<dim>> &quadrature_points =
1658 const unsigned int n_boundary_dofs =
1661 const unsigned int n_quadrature_points = quadrature.
size();
1665 n_quadrature_points);
1667 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1672 for (
unsigned int i = 0; i <= deg; ++i)
1675 weight * legendre_polynomials[i].value(
1676 quadrature_points[q_point](0));
1678 for (
unsigned int j = 0; j < deg; ++j)
1681 L_i * lobatto_polynomials[j + 2].value(
1682 quadrature_points[q_point](1));
1684 for (
unsigned int k = 0; k < deg; ++k)
1685 assembling_matrix((i * deg + j) * deg + k,
1687 l_j * lobatto_polynomials[k + 2].value(
1688 quadrature_points[q_point](2));
1694 assembling_matrix.
m());
1696 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1697 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1698 system_matrix_inv.
invert(system_matrix);
1701 solution.reinit(system_matrix_inv.
m(), 24);
1702 system_rhs.reinit(system_matrix_inv.
m(), 24);
1705 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1709 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1714 if (quadrature_points[q_point](0) < 0.5)
1716 if (quadrature_points[q_point](1) < 0.5)
1718 if (quadrature_points[q_point](2) < 0.5)
1721 2.0 * quadrature_points[q_point](0),
1722 2.0 * quadrature_points[q_point](1),
1723 2.0 * quadrature_points[q_point](2));
1725 tmp(0) += 2.0 * this->shape_value_component(
1726 dof, quadrature_point, 0);
1727 tmp(1) += 2.0 * this->shape_value_component(
1728 dof, quadrature_point, 1);
1729 tmp(2) += 2.0 * this->shape_value_component(
1730 dof, quadrature_point, 2);
1736 2.0 * quadrature_points[q_point](0),
1737 2.0 * quadrature_points[q_point](1),
1738 2.0 * quadrature_points[q_point](2) - 1.0);
1740 tmp(3) += 2.0 * this->shape_value_component(
1741 dof, quadrature_point, 0);
1742 tmp(4) += 2.0 * this->shape_value_component(
1743 dof, quadrature_point, 1);
1744 tmp(5) += 2.0 * this->shape_value_component(
1745 dof, quadrature_point, 2);
1749 else if (quadrature_points[q_point](2) < 0.5)
1752 2.0 * quadrature_points[q_point](0),
1753 2.0 * quadrature_points[q_point](1) - 1.0,
1754 2.0 * quadrature_points[q_point](2));
1756 tmp(6) += 2.0 * this->shape_value_component(
1757 dof, quadrature_point, 0);
1758 tmp(7) += 2.0 * this->shape_value_component(
1759 dof, quadrature_point, 1);
1760 tmp(8) += 2.0 * this->shape_value_component(
1761 dof, quadrature_point, 2);
1767 2.0 * quadrature_points[q_point](0),
1768 2.0 * quadrature_points[q_point](1) - 1.0,
1769 2.0 * quadrature_points[q_point](2) - 1.0);
1771 tmp(9) += 2.0 * this->shape_value_component(
1772 dof, quadrature_point, 0);
1773 tmp(10) += 2.0 * this->shape_value_component(
1774 dof, quadrature_point, 1);
1775 tmp(11) += 2.0 * this->shape_value_component(
1776 dof, quadrature_point, 2);
1780 else if (quadrature_points[q_point](1) < 0.5)
1782 if (quadrature_points[q_point](2) < 0.5)
1785 2.0 * quadrature_points[q_point](0) - 1.0,
1786 2.0 * quadrature_points[q_point](1),
1787 2.0 * quadrature_points[q_point](2));
1789 tmp(12) += 2.0 * this->shape_value_component(
1790 dof, quadrature_point, 0);
1791 tmp(13) += 2.0 * this->shape_value_component(
1792 dof, quadrature_point, 1);
1793 tmp(14) += 2.0 * this->shape_value_component(
1794 dof, quadrature_point, 2);
1800 2.0 * quadrature_points[q_point](0) - 1.0,
1801 2.0 * quadrature_points[q_point](1),
1802 2.0 * quadrature_points[q_point](2) - 1.0);
1804 tmp(15) += 2.0 * this->shape_value_component(
1805 dof, quadrature_point, 0);
1806 tmp(16) += 2.0 * this->shape_value_component(
1807 dof, quadrature_point, 1);
1808 tmp(17) += 2.0 * this->shape_value_component(
1809 dof, quadrature_point, 2);
1813 else if (quadrature_points[q_point](2) < 0.5)
1816 2.0 * quadrature_points[q_point](0) - 1.0,
1817 2.0 * quadrature_points[q_point](1) - 1.0,
1818 2.0 * quadrature_points[q_point](2));
1821 2.0 * this->shape_value_component(dof,
1825 2.0 * this->shape_value_component(dof,
1829 2.0 * this->shape_value_component(dof,
1837 2.0 * quadrature_points[q_point](0) - 1.0,
1838 2.0 * quadrature_points[q_point](1) - 1.0,
1839 2.0 * quadrature_points[q_point](2) - 1.0);
1842 2.0 * this->shape_value_component(dof,
1846 2.0 * this->shape_value_component(dof,
1850 2.0 * this->shape_value_component(dof,
1855 for (
unsigned int i = 0; i < 2; ++i)
1856 for (
unsigned int j = 0; j < 2; ++j)
1857 for (
unsigned int k = 0; k < 2; ++k)
1858 for (
unsigned int l = 0; l <= deg; ++l)
1860 tmp(3 * (i + 2 * (j + 2 * k))) -=
1861 this->restriction[index][2 * (2 * i + j) + k](
1862 (4 * i + j + 2) * this->degree + l, dof) *
1863 this->shape_value_component(
1864 (4 * i + j + 2) * this->degree + l,
1865 quadrature_points[q_point],
1867 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1868 this->restriction[index][2 * (2 * i + j) + k](
1869 (4 * i + k) * this->degree + l, dof) *
1870 this->shape_value_component(
1871 (4 * i + k) * this->degree + l,
1872 quadrature_points[q_point],
1874 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1875 this->restriction[index][2 * (2 * i + j) + k](
1876 (2 * (j + 4) + k) * this->degree + l, dof) *
1877 this->shape_value_component(
1878 (2 * (j + 4) + k) * this->degree + l,
1879 quadrature_points[q_point],
1882 for (
unsigned int m = 0; m < deg; ++m)
1884 tmp(3 * (i + 2 * (j + 2 * k))) -=
1885 this->restriction[index][2 * (2 * i + j) +
1887 ((2 * j + 5) * deg + m) * this->degree +
1890 this->shape_value_component(
1891 ((2 * j + 5) * deg + m) * this->degree +
1893 quadrature_points[q_point],
1895 tmp(3 * (i + 2 * (j + 2 * k))) -=
1896 this->restriction[index][2 * (2 * i + j) +
1898 (2 * (i + 4) * this->degree + l) * deg +
1901 this->shape_value_component(
1902 (2 * (i + 4) * this->degree + l) * deg +
1904 quadrature_points[q_point],
1906 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1907 this->restriction[index][2 * (2 * i + j) +
1909 (2 * k * this->degree + l) * deg + m +
1912 this->shape_value_component(
1913 (2 * k * this->degree + l) * deg + m +
1915 quadrature_points[q_point],
1917 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
1918 this->restriction[index][2 * (2 * i + j) +
1920 ((2 * i + 9) * deg + m) * this->degree +
1923 this->shape_value_component(
1924 ((2 * i + 9) * deg + m) * this->degree +
1926 quadrature_points[q_point],
1928 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1929 this->restriction[index][2 * (2 * i + j) +
1931 ((2 * k + 1) * deg + m) * this->degree +
1934 this->shape_value_component(
1935 ((2 * k + 1) * deg + m) * this->degree +
1937 quadrature_points[q_point],
1939 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
1940 this->restriction[index][2 * (2 * i + j) +
1942 (2 * (j + 2) * this->degree + l) * deg +
1945 this->shape_value_component(
1946 (2 * (j + 2) * this->degree + l) * deg +
1948 quadrature_points[q_point],
1953 tmp *= quadrature.
weight(q_point);
1955 for (
unsigned int i = 0; i <= deg; ++i)
1957 const double L_i_0 = legendre_polynomials[i].value(
1958 quadrature_points[q_point](0));
1959 const double L_i_1 = legendre_polynomials[i].value(
1960 quadrature_points[q_point](1));
1961 const double L_i_2 = legendre_polynomials[i].value(
1962 quadrature_points[q_point](2));
1964 for (
unsigned int j = 0; j < deg; ++j)
1966 const double l_j_0 =
1967 L_i_0 * lobatto_polynomials[j + 2].value(
1968 quadrature_points[q_point](1));
1969 const double l_j_1 =
1970 L_i_1 * lobatto_polynomials[j + 2].value(
1971 quadrature_points[q_point](0));
1972 const double l_j_2 =
1973 L_i_2 * lobatto_polynomials[j + 2].value(
1974 quadrature_points[q_point](0));
1976 for (
unsigned int k = 0; k < deg; ++k)
1978 const double l_k_0 =
1979 l_j_0 * lobatto_polynomials[k + 2].value(
1980 quadrature_points[q_point](2));
1981 const double l_k_1 =
1982 l_j_1 * lobatto_polynomials[k + 2].value(
1983 quadrature_points[q_point](2));
1984 const double l_k_2 =
1985 l_j_2 * lobatto_polynomials[k + 2].value(
1986 quadrature_points[q_point](1));
1988 for (
unsigned int l = 0; l < 8; ++l)
1990 system_rhs((i * deg + j) * deg + k,
1991 3 * l) += tmp(3 * l) * l_k_0;
1992 system_rhs((i * deg + j) * deg + k,
1994 tmp(3 * l + 1) * l_k_1;
1995 system_rhs((i * deg + j) * deg + k,
1997 tmp(3 * l + 2) * l_k_2;
2004 system_matrix_inv.
mmult(solution, system_rhs);
2006 for (
unsigned int i = 0; i < 2; ++i)
2007 for (
unsigned int j = 0; j < 2; ++j)
2008 for (
unsigned int k = 0; k < 2; ++k)
2009 for (
unsigned int l = 0; l <= deg; ++l)
2010 for (
unsigned int m = 0; m < deg; ++m)
2011 for (
unsigned int n = 0; n < deg; ++n)
2014 solution((l * deg + m) * deg + n,
2015 3 * (i + 2 * (j + 2 * k)))) >
2017 this->restriction[index][2 * (2 * i + j) + k](
2018 (l * deg + m) * deg + n + n_boundary_dofs,
2019 dof) = solution((l * deg + m) * deg + n,
2020 3 * (i + 2 * (j + 2 * k)));
2023 solution((l * deg + m) * deg + n,
2024 3 * (i + 2 * (j + 2 * k)) + 1)) >
2026 this->restriction[index][2 * (2 * i + j) + k](
2027 (l + (m + deg) * this->degree) * deg + n +
2030 solution((l * deg + m) * deg + n,
2031 3 * (i + 2 * (j + 2 * k)) + 1);
2034 solution((l * deg + m) * deg + n,
2035 3 * (i + 2 * (j + 2 * k)) + 2)) >
2037 this->restriction[index][2 * (2 * i + j) + k](
2039 ((m + 2 * deg) * deg + n) * this->degree +
2042 solution((l * deg + m) * deg + n,
2043 3 * (i + 2 * (j + 2 * k)) + 2);
3150 std::vector<double> & nodal_values)
const
3155 const unsigned int face_no = 0;
3157 const unsigned int deg = this->degree - 1;
3158 Assert(support_point_values.size() == this->generalized_support_points.size(),
3160 this->generalized_support_points.size()));
3161 Assert(support_point_values[0].size() == this->n_components(),
3163 this->n_components()));
3164 Assert(nodal_values.size() == this->n_dofs_per_cell(),
3166 std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3174 const QGauss<1> reference_edge_quadrature(this->degree);
3175 const unsigned int n_edge_points = reference_edge_quadrature.
size();
3177 for (
unsigned int i = 0; i < 2; ++i)
3178 for (
unsigned int j = 0; j < 2; ++j)
3180 for (
unsigned int q_point = 0; q_point < n_edge_points;
3182 nodal_values[(i + 2 * j) * this->degree] +=
3183 reference_edge_quadrature.
weight(q_point) *
3184 support_point_values[q_point + (i + 2 * j) * n_edge_points]
3190 if (
std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3191 nodal_values[(i + 2 * j) * this->degree] = 0.0;
3204 if (this->degree > 1)
3209 const std::vector<Polynomials::Polynomial<double>>
3210 &lobatto_polynomials =
3214 std::vector<Polynomials::Polynomial<double>>
3215 lobatto_polynomials_grad(this->degree);
3217 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3218 lobatto_polynomials_grad[i] =
3219 lobatto_polynomials[i + 1].derivative();
3224 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3225 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3226 for (
unsigned int q_point = 0; q_point < n_edge_points;
3228 system_matrix(i, j) +=
3229 boundary_weights(q_point, j) *
3230 lobatto_polynomials_grad[i + 1].value(
3231 this->generalized_face_support_points[face_no][q_point](
3237 system_matrix_inv.
invert(system_matrix);
3244 for (
unsigned int line = 0;
3245 line < GeometryInfo<dim>::lines_per_cell;
3251 for (
unsigned int q_point = 0; q_point < n_edge_points;
3255 support_point_values[line * n_edge_points + q_point]
3256 [line_coordinate[line]] -
3257 nodal_values[line * this->degree] *
3258 this->shape_value_component(
3259 line * this->degree,
3260 this->generalized_support_points[line *
3263 line_coordinate[line]);
3265 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3266 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3269 system_matrix_inv.
vmult(solution, system_rhs);
3275 for (
unsigned int i = 0; i < solution.
size(); ++i)
3277 nodal_values[line * this->degree + i + 1] = solution(i);
3289 const QGauss<dim> reference_quadrature(this->degree);
3290 const unsigned int n_interior_points =
3291 reference_quadrature.
size();
3292 const std::vector<Polynomials::Polynomial<double>>
3293 &legendre_polynomials =
3297 system_matrix.reinit((this->degree - 1) * this->degree,
3298 (this->degree - 1) * this->degree);
3301 for (
unsigned int i = 0; i < this->degree; ++i)
3302 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3303 for (
unsigned int k = 0; k < this->degree; ++k)
3304 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3305 for (
unsigned int q_point = 0;
3306 q_point < n_interior_points;
3308 system_matrix(i * (this->degree - 1) + j,
3309 k * (this->degree - 1) + l) +=
3310 reference_quadrature.
weight(q_point) *
3311 legendre_polynomials[i].value(
3312 this->generalized_support_points
3314 n_edge_points](0)) *
3315 lobatto_polynomials[j + 2].value(
3316 this->generalized_support_points
3318 n_edge_points](1)) *
3319 lobatto_polynomials_grad[k].value(
3320 this->generalized_support_points
3322 n_edge_points](0)) *
3323 lobatto_polynomials[l + 2].value(
3324 this->generalized_support_points
3328 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3329 system_matrix_inv.
invert(system_matrix);
3333 system_rhs.
reinit(system_matrix_inv.
m());
3336 for (
unsigned int q_point = 0; q_point < n_interior_points;
3340 support_point_values[q_point +
3344 for (
unsigned int i = 0; i < 2; ++i)
3345 for (
unsigned int j = 0; j <= deg; ++j)
3346 tmp -= nodal_values[(i + 2) * this->degree + j] *
3347 this->shape_value_component(
3348 (i + 2) * this->degree + j,
3349 this->generalized_support_points
3354 for (
unsigned int i = 0; i <= deg; ++i)
3355 for (
unsigned int j = 0; j < deg; ++j)
3356 system_rhs(i * deg + j) +=
3357 reference_quadrature.
weight(q_point) * tmp *
3358 lobatto_polynomials_grad[i].value(
3359 this->generalized_support_points
3361 n_edge_points](0)) *
3362 lobatto_polynomials[j + 2].value(
3363 this->generalized_support_points
3368 solution.
reinit(system_matrix.
m());
3369 system_matrix_inv.
vmult(solution, system_rhs);
3375 for (
unsigned int i = 0; i <= deg; ++i)
3376 for (
unsigned int j = 0; j < deg; ++j)
3377 if (
std::abs(solution(i * deg + j)) > 1e-14)
3380 solution(i * deg + j);
3387 for (
unsigned int q_point = 0; q_point < n_interior_points;
3391 support_point_values[q_point +
3395 for (
unsigned int i = 0; i < 2; ++i)
3396 for (
unsigned int j = 0; j <= deg; ++j)
3397 tmp -= nodal_values[i * this->degree + j] *
3398 this->shape_value_component(
3399 i * this->degree + j,
3400 this->generalized_support_points
3405 for (
unsigned int i = 0; i <= deg; ++i)
3406 for (
unsigned int j = 0; j < deg; ++j)
3407 system_rhs(i * deg + j) +=
3408 reference_quadrature.
weight(q_point) * tmp *
3409 lobatto_polynomials_grad[i].value(
3410 this->generalized_support_points
3412 n_edge_points](1)) *
3413 lobatto_polynomials[j + 2].value(
3414 this->generalized_support_points
3419 system_matrix_inv.
vmult(solution, system_rhs);
3425 for (
unsigned int i = 0; i <= deg; ++i)
3426 for (
unsigned int j = 0; j < deg; ++j)
3427 if (
std::abs(solution(i * deg + j)) > 1e-14)
3430 this->degree] = solution(i * deg + j);
3440 const QGauss<1> reference_edge_quadrature(this->degree);
3441 const unsigned int n_edge_points = reference_edge_quadrature.
size();
3443 for (
unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3445 for (
unsigned int i = 0; i < 4; ++i)
3446 nodal_values[(i + 8) * this->degree] +=
3447 reference_edge_quadrature.
weight(q_point) *
3448 support_point_values[q_point + (i + 8) * n_edge_points][2];
3450 for (
unsigned int i = 0; i < 2; ++i)
3451 for (
unsigned int j = 0; j < 2; ++j)
3452 for (
unsigned int k = 0; k < 2; ++k)
3453 nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3454 reference_edge_quadrature.
weight(q_point) *
3455 support_point_values[q_point + (i + 2 * (2 * j + k)) *
3456 n_edge_points][1 - k];
3463 for (
unsigned int i = 0; i < 4; ++i)
3464 if (
std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3465 nodal_values[(i + 8) * this->degree] = 0.0;
3467 for (
unsigned int i = 0; i < 2; ++i)
3468 for (
unsigned int j = 0; j < 2; ++j)
3469 for (
unsigned int k = 0; k < 2; ++k)
3471 nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3473 nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3484 if (this->degree > 1)
3489 const std::vector<Polynomials::Polynomial<double>>
3490 &lobatto_polynomials =
3494 std::vector<Polynomials::Polynomial<double>>
3495 lobatto_polynomials_grad(this->degree);
3497 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3498 lobatto_polynomials_grad[i] =
3499 lobatto_polynomials[i + 1].derivative();
3504 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3505 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3506 for (
unsigned int q_point = 0; q_point < n_edge_points;
3508 system_matrix(i, j) +=
3509 boundary_weights(q_point, j) *
3510 lobatto_polynomials_grad[i + 1].value(
3511 this->generalized_face_support_points[face_no][q_point](
3517 system_matrix_inv.
invert(system_matrix);
3521 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3525 for (
unsigned int line = 0;
3526 line < GeometryInfo<dim>::lines_per_cell;
3532 for (
unsigned int q_point = 0; q_point < this->degree;
3536 support_point_values[line * this->degree + q_point]
3537 [line_coordinate[line]] -
3538 nodal_values[line * this->degree] *
3539 this->shape_value_component(
3540 line * this->degree,
3542 ->generalized_support_points[line * this->degree +
3544 line_coordinate[line]);
3546 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3547 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3550 system_matrix_inv.
vmult(solution, system_rhs);
3556 for (
unsigned int i = 0; i < solution.
size(); ++i)
3558 nodal_values[line * this->degree + i + 1] = solution(i);
3569 const std::vector<Polynomials::Polynomial<double>>
3570 &legendre_polynomials =
3573 const unsigned int n_face_points = n_edge_points * n_edge_points;
3575 system_matrix.reinit((this->degree - 1) * this->degree,
3576 (this->degree - 1) * this->degree);
3579 for (
unsigned int i = 0; i < this->degree; ++i)
3580 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3581 for (
unsigned int k = 0; k < this->degree; ++k)
3582 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3583 for (
unsigned int q_point = 0; q_point < n_face_points;
3585 system_matrix(i * (this->degree - 1) + j,
3586 k * (this->degree - 1) + l) +=
3587 boundary_weights(q_point + n_edge_points,
3588 2 * (k * (this->degree - 1) + l)) *
3589 legendre_polynomials[i].value(
3590 this->generalized_face_support_points
3591 [face_no][q_point + 4 * n_edge_points](0)) *
3592 lobatto_polynomials[j + 2].value(
3593 this->generalized_face_support_points
3594 [face_no][q_point + 4 * n_edge_points](1));
3596 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3597 system_matrix_inv.
invert(system_matrix);
3598 solution.
reinit(system_matrix.
m());
3599 system_rhs.
reinit(system_matrix.
m());
3603 {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3620 for (
unsigned int q_point = 0; q_point < n_face_points;
3624 support_point_values[q_point +
3627 [face_coordinates[face][0]];
3629 for (
unsigned int i = 0; i < 2; ++i)
3630 for (
unsigned int j = 0; j <= deg; ++j)
3634 this->shape_value_component(
3636 this->generalized_support_points
3639 face_coordinates[face][0]);
3641 for (
unsigned int i = 0; i <= deg; ++i)
3642 for (
unsigned int j = 0; j < deg; ++j)
3643 system_rhs(i * deg + j) +=
3644 boundary_weights(q_point + n_edge_points,
3645 2 * (i * deg + j)) *
3649 system_matrix_inv.
vmult(solution, system_rhs);
3655 for (
unsigned int i = 0; i <= deg; ++i)
3656 for (
unsigned int j = 0; j < deg; ++j)
3657 if (
std::abs(solution(i * deg + j)) > 1e-14)
3658 nodal_values[(2 * face * this->degree + i +
3662 solution(i * deg + j);
3669 for (
unsigned int q_point = 0; q_point < n_face_points;
3673 support_point_values[q_point +
3676 [face_coordinates[face][1]];
3678 for (
unsigned int i = 2;
3679 i < GeometryInfo<dim>::lines_per_face;
3681 for (
unsigned int j = 0; j <= deg; ++j)
3685 this->shape_value_component(
3687 this->generalized_support_points
3690 face_coordinates[face][1]);
3692 for (
unsigned int i = 0; i <= deg; ++i)
3693 for (
unsigned int j = 0; j < deg; ++j)
3694 system_rhs(i * deg + j) +=
3695 boundary_weights(q_point + n_edge_points,
3696 2 * (i * deg + j) + 1) *
3700 system_matrix_inv.
vmult(solution, system_rhs);
3706 for (
unsigned int i = 0; i <= deg; ++i)
3707 for (
unsigned int j = 0; j < deg; ++j)
3708 if (
std::abs(solution(i * deg + j)) > 1e-14)
3709 nodal_values[((2 * face + 1) * deg + j +
3712 i] = solution(i * deg + j);
3720 const QGauss<dim> reference_quadrature(this->degree);
3721 const unsigned int n_interior_points =
3722 reference_quadrature.
size();
3726 system_matrix.reinit(this->degree * deg * deg,
3727 this->degree * deg * deg);
3730 for (
unsigned int i = 0; i <= deg; ++i)
3731 for (
unsigned int j = 0; j < deg; ++j)
3732 for (
unsigned int k = 0; k < deg; ++k)
3733 for (
unsigned int l = 0; l <= deg; ++l)
3734 for (
unsigned int m = 0; m < deg; ++m)
3735 for (
unsigned int n = 0; n < deg; ++n)
3736 for (
unsigned int q_point = 0;
3737 q_point < n_interior_points;
3739 system_matrix((i * deg + j) * deg + k,
3740 (l * deg + m) * deg + n) +=
3741 reference_quadrature.
weight(q_point) *
3742 legendre_polynomials[i].value(
3743 this->generalized_support_points
3748 n_face_points](0)) *
3749 lobatto_polynomials[j + 2].value(
3750 this->generalized_support_points
3755 n_face_points](1)) *
3756 lobatto_polynomials[k + 2].value(
3757 this->generalized_support_points
3762 n_face_points](2)) *
3763 lobatto_polynomials_grad[l].value(
3764 this->generalized_support_points
3769 n_face_points](0)) *
3770 lobatto_polynomials[m + 2].value(
3771 this->generalized_support_points
3776 n_face_points](1)) *
3777 lobatto_polynomials[n + 2].value(
3778 this->generalized_support_points
3785 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3786 system_matrix_inv.
invert(system_matrix);
3788 system_rhs.
reinit(system_matrix.
m());
3791 for (
unsigned int q_point = 0; q_point < n_interior_points;
3795 support_point_values[q_point +
3801 for (
unsigned int i = 0; i <= deg; ++i)
3803 for (
unsigned int j = 0; j < 2; ++j)
3804 for (
unsigned int k = 0; k < 2; ++k)
3806 nodal_values[i + (j + 4 * k + 2) * this->degree] *
3807 this->shape_value_component(
3808 i + (j + 4 * k + 2) * this->degree,
3809 this->generalized_support_points
3817 for (
unsigned int j = 0; j < deg; ++j)
3818 for (
unsigned int k = 0; k < 4; ++k)
3820 nodal_values[(i + 2 * (k + 2) * this->degree +
3825 this->shape_value_component(
3826 (i + 2 * (k + 2) * this->degree +
3830 this->generalized_support_points
3839 for (
unsigned int i = 0; i <= deg; ++i)
3840 for (
unsigned int j = 0; j < deg; ++j)
3841 for (
unsigned int k = 0; k < deg; ++k)
3842 system_rhs((i * deg + j) * deg + k) +=
3843 reference_quadrature.
weight(q_point) * tmp *
3844 lobatto_polynomials_grad[i].value(
3845 this->generalized_support_points
3850 n_face_points](0)) *
3851 lobatto_polynomials[j + 2].value(
3852 this->generalized_support_points
3857 n_face_points](1)) *
3858 lobatto_polynomials[k + 2].value(
3859 this->generalized_support_points
3868 system_matrix_inv.
vmult(solution, system_rhs);
3874 for (
unsigned int i = 0; i <= deg; ++i)
3875 for (
unsigned int j = 0; j < deg; ++j)
3876 for (
unsigned int k = 0; k < deg; ++k)
3877 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3884 solution((i * deg + j) * deg + k);
3889 for (
unsigned int q_point = 0; q_point < n_interior_points;
3893 support_point_values[q_point +
3899 for (
unsigned int i = 0; i <= deg; ++i)
3900 for (
unsigned int j = 0; j < 2; ++j)
3902 for (
unsigned int k = 0; k < 2; ++k)
3903 tmp -= nodal_values[i + (4 * j + k) * this->degree] *
3904 this->shape_value_component(
3905 i + (4 * j + k) * this->degree,
3906 this->generalized_support_points
3914 for (
unsigned int k = 0; k < deg; ++k)
3916 nodal_values[(i + 2 * j * this->degree +
3921 this->shape_value_component(
3922 (i + 2 * j * this->degree +
3926 this->generalized_support_points
3934 ((2 * j + 9) * deg + k +
3937 this->shape_value_component(
3938 i + ((2 * j + 9) * deg + k +
3941 this->generalized_support_points
3950 for (
unsigned int i = 0; i <= deg; ++i)
3951 for (
unsigned int j = 0; j < deg; ++j)
3952 for (
unsigned int k = 0; k < deg; ++k)
3953 system_rhs((i * deg + j) * deg + k) +=
3954 reference_quadrature.
weight(q_point) * tmp *
3955 lobatto_polynomials_grad[i].value(
3956 this->generalized_support_points
3961 n_face_points](1)) *
3962 lobatto_polynomials[j + 2].value(
3963 this->generalized_support_points
3968 n_face_points](0)) *
3969 lobatto_polynomials[k + 2].value(
3970 this->generalized_support_points
3978 system_matrix_inv.
vmult(solution, system_rhs);
3984 for (
unsigned int i = 0; i <= deg; ++i)
3985 for (
unsigned int j = 0; j < deg; ++j)
3986 for (
unsigned int k = 0; k < deg; ++k)
3987 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
3988 nodal_values[((i + this->degree +
3995 solution((i * deg + j) * deg + k);
4000 for (
unsigned int q_point = 0; q_point < n_interior_points;
4004 support_point_values[q_point +
4010 for (
unsigned int i = 0; i <= deg; ++i)
4011 for (
unsigned int j = 0; j < 4; ++j)
4013 tmp -= nodal_values[i + (j + 8) * this->degree] *
4014 this->shape_value_component(
4015 i + (j + 8) * this->degree,
4016 this->generalized_support_points
4024 for (
unsigned int k = 0; k < deg; ++k)
4027 ((2 * j + 1) * deg + k +
4030 this->shape_value_component(
4031 i + ((2 * j + 1) * deg + k +
4034 this->generalized_support_points
4043 for (
unsigned int i = 0; i <= deg; ++i)
4044 for (
unsigned int j = 0; j < deg; ++j)
4045 for (
unsigned int k = 0; k < deg; ++k)
4046 system_rhs((i * deg + j) * deg + k) +=
4047 reference_quadrature.
weight(q_point) * tmp *
4048 lobatto_polynomials_grad[i].value(
4049 this->generalized_support_points
4054 n_face_points](2)) *
4055 lobatto_polynomials[j + 2].value(
4056 this->generalized_support_points
4061 n_face_points](0)) *
4062 lobatto_polynomials[k + 2].value(
4063 this->generalized_support_points
4071 system_matrix_inv.
vmult(solution, system_rhs);
4077 for (
unsigned int i = 0; i <= deg; ++i)
4078 for (
unsigned int j = 0; j < deg; ++j)
4079 for (
unsigned int k = 0; k < deg; ++k)
4080 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4086 this->degree] = solution((i * deg + j) * deg + k);