334 static const int c_swap_table_0 = 0;
336 static const int c_swap_table_1[8][3][2] = {
362 static const int c_swap_table_2[8][3][6] = {
363 {{-1, -1, -1, -1, -1, -1},
372 {{-1, -1, -1, -1, -1, -1},
375 {{-1, -1, -1, -1, -1, -1},
384 {{-1, -1, -1, -1, -1, -1},
386 {0, 0, 0, 1, 1, 1}}};
388 static const int c_swap_table_3[8][3][12] = {
390 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
391 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
392 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
393 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11},
394 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
395 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
396 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11},
397 {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
398 {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}},
399 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
400 {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
401 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}},
402 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
403 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
404 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
405 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11},
406 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
407 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
408 {{0, 4, 8, 1, 5, 9, 2, 6, 10, 3, 7, 11},
409 {0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0},
410 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0}},
411 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
412 {1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0},
413 {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}}};
415 static const int c_swap_table_4[8][3][20] = {
419 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
420 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
421 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
422 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
423 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
424 12, 17, 3, 8, 13, 18, 4, 9, 14, 19},
425 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
426 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
427 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
428 12, 17, 3, 8, 13, 18, 4, 9, 14, 19},
429 {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
430 {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}},
431 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
432 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
433 {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
434 {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}},
435 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
436 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
437 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
438 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
439 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
440 12, 17, 3, 8, 13, 18, 4, 9, 14, 19},
441 {1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0},
442 {1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0}},
443 {{0, 5, 10, 15, 1, 6, 11, 16, 2, 7,
444 12, 17, 3, 8, 13, 18, 4, 9, 14, 19},
445 {1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1},
446 {0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}},
447 {{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
448 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
449 {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
450 {1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1}}};
452 static const int *swap_table_array[5] = {&c_swap_table_0,
453 &c_swap_table_1[0][0][0],
454 &c_swap_table_2[0][0][0],
455 &c_swap_table_3[0][0][0],
456 &c_swap_table_4[0][0][0]};
458 static const int row_length[5] = {0, 2, 6, 12, 20};
459 static const int table_size[5] = {
460 0, 8 * 3 * 2, 8 * 3 * 6, 8 * 3 * 12, 8 * 3 * 20};
472 const unsigned int k = this->tensor_degree() - 1;
484 const unsigned int face_no = 0;
487 this->adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
488 this->reference_cell().n_face_orientations(face_no) *
489 this->n_dofs_per_quad(face_no),
493 this->adjust_quad_dof_sign_for_face_orientation_table[0].n_elements() ==
494 this->reference_cell().n_face_orientations(face_no) *
495 this->n_dofs_per_quad(face_no),
501 const int *swap_table = swap_table_array[k];
503 const unsigned int half_dofs = k * (k + 1);
505 const int rl = row_length[k];
507 int offset = table_size[0];
513 for (
const bool face_orientation : {
false,
true})
514 for (
const bool face_rotation : {
false,
true})
515 for (
const bool face_flip : {
false,
true})
530 for (
unsigned int indx_x = 0; indx_x < half_dofs; indx_x++)
532 offset = 3 * rl * case_no + 0 * rl + indx_x;
534 value = *(swap_table + offset);
540 const unsigned int indx_y =
541 half_dofs +
static_cast<unsigned int>(value);
545 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
546 indx_x, case_no) = indx_y - indx_x;
549 ->adjust_quad_dof_index_for_face_orientation_table[face_no](
550 indx_y, case_no) = indx_x - indx_y;
554 offset = 3 * rl * case_no + 1 * rl + indx_x;
556 value = *(swap_table + offset);
559 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
560 indx_x, case_no) =
static_cast<bool>(value);
563 offset = 3 * rl * case_no + 2 * rl + indx_x;
565 value = *(swap_table + offset);
568 this->adjust_quad_dof_sign_for_face_orientation_table[face_no](
569 indx_x + half_dofs, case_no) =
static_cast<bool>(value);
922 const QGauss<1> edge_quadrature(2 * this->degree);
923 const std::vector<Point<1>> &edge_quadrature_points =
924 edge_quadrature.get_points();
925 const unsigned int n_edge_quadrature_points = edge_quadrature.size();
936 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
937 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
940 const double weight = 2.0 * edge_quadrature.weight(q_point);
942 if (edge_quadrature_points[q_point][0] < 0.5)
945 0.0, 2.0 * edge_quadrature_points[q_point][0]);
947 this->restriction[index][0](0, dof) +=
949 this->shape_value_component(dof, quadrature_point, 1);
950 quadrature_point[0] = 1.0;
951 this->restriction[index][1](this->degree, dof) +=
953 this->shape_value_component(dof, quadrature_point, 1);
954 quadrature_point[0] = quadrature_point[1];
955 quadrature_point[1] = 0.0;
956 this->restriction[index][0](2 * this->degree, dof) +=
958 this->shape_value_component(dof, quadrature_point, 0);
959 quadrature_point[1] = 1.0;
960 this->restriction[index][2](3 * this->degree, dof) +=
962 this->shape_value_component(dof, quadrature_point, 0);
968 0.0, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
970 this->restriction[index][2](0, dof) +=
972 this->shape_value_component(dof, quadrature_point, 1);
973 quadrature_point[0] = 1.0;
974 this->restriction[index][3](this->degree, dof) +=
976 this->shape_value_component(dof, quadrature_point, 1);
977 quadrature_point[0] = quadrature_point[1];
978 quadrature_point[1] = 0.0;
979 this->restriction[index][1](2 * this->degree, dof) +=
981 this->shape_value_component(dof, quadrature_point, 0);
982 quadrature_point[1] = 1.0;
983 this->restriction[index][3](3 * this->degree, dof) +=
985 this->shape_value_component(dof, quadrature_point, 0);
993 if (this->degree > 1)
995 const unsigned int deg = this->degree - 1;
996 const std::vector<Polynomials::Polynomial<double>>
997 &legendre_polynomials =
1003 n_edge_quadrature_points);
1005 for (
unsigned int q_point = 0;
1006 q_point < n_edge_quadrature_points;
1009 const double weight =
1010 std::sqrt(edge_quadrature.weight(q_point));
1012 for (
unsigned int i = 0; i < deg; ++i)
1013 assembling_matrix(i, q_point) =
1014 weight * legendre_polynomials[i + 1].value(
1015 edge_quadrature_points[q_point][0]);
1020 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1021 system_matrix_inv.
invert(system_matrix);
1028 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1029 for (
unsigned int i = 0; i < 2; ++i)
1033 for (
unsigned int q_point = 0;
1034 q_point < n_edge_quadrature_points;
1037 const double weight = edge_quadrature.weight(q_point);
1039 i, edge_quadrature_points[q_point][0]);
1041 edge_quadrature_points[q_point][0], i);
1043 if (edge_quadrature_points[q_point][0] < 0.5)
1046 i, 2.0 * edge_quadrature_points[q_point][0]);
1050 (2.0 * this->shape_value_component(
1051 dof, quadrature_point_2, 1) -
1052 this->restriction[index][i](i * this->degree,
1054 this->shape_value_component(i * this->degree,
1059 this->restriction[index][i + 2](i * this->degree,
1061 this->shape_value_component(i * this->degree,
1065 2.0 * edge_quadrature_points[q_point][0], i);
1068 (2.0 * this->shape_value_component(
1069 dof, quadrature_point_2, 0) -
1070 this->restriction[index][2 * i]((i + 2) *
1073 this->shape_value_component((i + 2) *
1079 this->restriction[index][2 * i + 1](
1080 (i + 2) * this->degree, dof) *
1081 this->shape_value_component(
1082 (i + 2) * this->degree, quadrature_point_1, 0);
1089 this->restriction[index][i](i * this->degree,
1091 this->shape_value_component(i * this->degree,
1097 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1101 (2.0 * this->shape_value_component(
1102 dof, quadrature_point_2, 1) -
1103 this->restriction[index][i + 2](i * this->degree,
1105 this->shape_value_component(i * this->degree,
1110 this->restriction[index][2 * i]((i + 2) *
1113 this->shape_value_component(
1114 (i + 2) * this->degree, quadrature_point_1, 0);
1116 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1120 (2.0 * this->shape_value_component(
1121 dof, quadrature_point_2, 0) -
1122 this->restriction[index][2 * i + 1](
1123 (i + 2) * this->degree, dof) *
1124 this->shape_value_component((i + 2) *
1130 for (
unsigned int j = 0; j < this->degree - 1; ++j)
1133 legendre_polynomials[j + 1].value(
1134 edge_quadrature_points[q_point][0]);
1136 for (
unsigned int k = 0; k < tmp.
size(); ++k)
1137 system_rhs(j, k) += tmp(k) * L_j;
1141 system_matrix_inv.
mmult(solution, system_rhs);
1143 for (
unsigned int j = 0; j < this->degree - 1; ++j)
1144 for (
unsigned int k = 0; k < 2; ++k)
1146 if (
std::abs(solution(j, k)) > 1e-14)
1147 this->restriction[index][i + 2 * k](
1148 i * this->degree + j + 1, dof) = solution(j, k);
1150 if (
std::abs(solution(j, k + 2)) > 1e-14)
1151 this->restriction[index][2 * i + k](
1152 (i + 2) * this->degree + j + 1, dof) =
1158 const std::vector<Point<dim>> &quadrature_points =
1159 quadrature.get_points();
1160 const std::vector<Polynomials::Polynomial<double>>
1161 &lobatto_polynomials =
1163 const unsigned int n_boundary_dofs =
1165 const unsigned int n_quadrature_points = quadrature.size();
1170 n_quadrature_points);
1172 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1175 const double weight =
std::sqrt(quadrature.weight(q_point));
1177 for (
unsigned int i = 0; i < this->degree; ++i)
1180 weight * legendre_polynomials[i].value(
1181 quadrature_points[q_point][0]);
1183 for (
unsigned int j = 0; j < this->degree - 1; ++j)
1184 assembling_matrix(i * (this->degree - 1) + j,
1186 L_i * lobatto_polynomials[j + 2].value(
1187 quadrature_points[q_point][1]);
1192 assembling_matrix.
m());
1194 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1195 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1196 system_matrix_inv.
invert(system_matrix);
1199 solution.reinit(system_matrix_inv.
m(), 8);
1200 system_rhs.reinit(system_matrix_inv.
m(), 8);
1203 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1207 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
1212 if (quadrature_points[q_point][0] < 0.5)
1214 if (quadrature_points[q_point][1] < 0.5)
1217 2.0 * quadrature_points[q_point][0],
1218 2.0 * quadrature_points[q_point][1]);
1220 tmp(0) += 2.0 * this->shape_value_component(
1221 dof, quadrature_point, 0);
1222 tmp(1) += 2.0 * this->shape_value_component(
1223 dof, quadrature_point, 1);
1229 2.0 * quadrature_points[q_point][0],
1230 2.0 * quadrature_points[q_point][1] - 1.0);
1232 tmp(4) += 2.0 * this->shape_value_component(
1233 dof, quadrature_point, 0);
1234 tmp(5) += 2.0 * this->shape_value_component(
1235 dof, quadrature_point, 1);
1239 else if (quadrature_points[q_point][1] < 0.5)
1242 2.0 * quadrature_points[q_point][0] - 1.0,
1243 2.0 * quadrature_points[q_point][1]);
1246 2.0 * this->shape_value_component(dof,
1250 2.0 * this->shape_value_component(dof,
1258 2.0 * quadrature_points[q_point][0] - 1.0,
1259 2.0 * quadrature_points[q_point][1] - 1.0);
1262 2.0 * this->shape_value_component(dof,
1266 2.0 * this->shape_value_component(dof,
1271 for (
unsigned int i = 0; i < 2; ++i)
1272 for (
unsigned int j = 0; j < this->degree; ++j)
1275 this->restriction[index][i](j + 2 * this->degree,
1277 this->shape_value_component(
1278 j + 2 * this->degree,
1279 quadrature_points[q_point],
1282 this->restriction[index][i](i * this->degree + j,
1284 this->shape_value_component(
1285 i * this->degree + j,
1286 quadrature_points[q_point],
1288 tmp(2 * (i + 2)) -= this->restriction[index][i + 2](
1289 j + 3 * this->degree, dof) *
1290 this->shape_value_component(
1291 j + 3 * this->degree,
1292 quadrature_points[q_point],
1294 tmp(2 * i + 5) -= this->restriction[index][i + 2](
1295 i * this->degree + j, dof) *
1296 this->shape_value_component(
1297 i * this->degree + j,
1298 quadrature_points[q_point],
1302 tmp *= quadrature.weight(q_point);
1304 for (
unsigned int i = 0; i < this->degree; ++i)
1306 const double L_i_0 = legendre_polynomials[i].value(
1307 quadrature_points[q_point][0]);
1308 const double L_i_1 = legendre_polynomials[i].value(
1309 quadrature_points[q_point][1]);
1311 for (
unsigned int j = 0; j < this->degree - 1; ++j)
1313 const double l_j_0 =
1314 L_i_0 * lobatto_polynomials[j + 2].value(
1315 quadrature_points[q_point][1]);
1316 const double l_j_1 =
1317 L_i_1 * lobatto_polynomials[j + 2].value(
1318 quadrature_points[q_point][0]);
1320 for (
unsigned int k = 0; k < 4; ++k)
1322 system_rhs(i * (this->degree - 1) + j,
1323 2 * k) += tmp(2 * k) * l_j_0;
1324 system_rhs(i * (this->degree - 1) + j,
1326 tmp(2 * k + 1) * l_j_1;
1332 system_matrix_inv.
mmult(solution, system_rhs);
1334 for (
unsigned int i = 0; i < this->degree; ++i)
1335 for (
unsigned int j = 0; j < this->degree - 1; ++j)
1336 for (
unsigned int k = 0; k < 4; ++k)
1338 if (
std::abs(solution(i * (this->degree - 1) + j,
1340 this->restriction[index][k](i * (this->degree - 1) +
1341 j + n_boundary_dofs,
1343 solution(i * (this->degree - 1) + j, 2 * k);
1345 if (
std::abs(solution(i * (this->degree - 1) + j,
1346 2 * k + 1)) > 1e-14)
1347 this->restriction[index][k](
1348 i + (this->degree - 1 + j) * this->degree +
1351 solution(i * (this->degree - 1) + j, 2 * k + 1);
1365 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1366 for (
unsigned int q_point = 0; q_point < n_edge_quadrature_points;
1369 const double weight = 2.0 * edge_quadrature.weight(q_point);
1371 if (edge_quadrature_points[q_point][0] < 0.5)
1372 for (
unsigned int i = 0; i < 2; ++i)
1373 for (
unsigned int j = 0; j < 2; ++j)
1376 i, 2.0 * edge_quadrature_points[q_point][0], j);
1378 this->restriction[index][i + 4 * j]((i + 4 * j) *
1382 this->shape_value_component(dof, quadrature_point, 1);
1384 Point<dim>(2.0 * edge_quadrature_points[q_point][0],
1387 this->restriction[index][2 * (i + 2 * j)](
1388 (i + 4 * j + 2) * this->degree, dof) +=
1390 this->shape_value_component(dof, quadrature_point, 0);
1394 2.0 * edge_quadrature_points[q_point][0]);
1395 this->restriction[index][i + 2 * j]((i + 2 * (j + 4)) *
1399 this->shape_value_component(dof, quadrature_point, 2);
1403 for (
unsigned int i = 0; i < 2; ++i)
1404 for (
unsigned int j = 0; j < 2; ++j)
1407 i, 2.0 * edge_quadrature_points[q_point][0] - 1.0, j);
1409 this->restriction[index][i + 4 * j + 2]((i + 4 * j) *
1413 this->shape_value_component(dof, quadrature_point, 1);
1415 2.0 * edge_quadrature_points[q_point][0] - 1.0, i, j);
1416 this->restriction[index][2 * (i + 2 * j) + 1](
1417 (i + 4 * j + 2) * this->degree, dof) +=
1419 this->shape_value_component(dof, quadrature_point, 0);
1421 i, j, 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1422 this->restriction[index][i + 2 * (j + 2)](
1423 (i + 2 * (j + 4)) * this->degree, dof) +=
1425 this->shape_value_component(dof, quadrature_point, 2);
1433 if (this->degree > 1)
1435 const unsigned int deg = this->degree - 1;
1436 const std::vector<Polynomials::Polynomial<double>>
1437 &legendre_polynomials =
1443 n_edge_quadrature_points);
1445 for (
unsigned int q_point = 0;
1446 q_point < n_edge_quadrature_points;
1449 const double weight =
1450 std::sqrt(edge_quadrature.weight(q_point));
1452 for (
unsigned int i = 0; i < deg; ++i)
1453 assembling_matrix(i, q_point) =
1454 weight * legendre_polynomials[i + 1].value(
1455 edge_quadrature_points[q_point][0]);
1460 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1461 system_matrix_inv.
invert(system_matrix);
1468 for (
unsigned int i = 0; i < 2; ++i)
1469 for (
unsigned int j = 0; j < 2; ++j)
1470 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell();
1475 for (
unsigned int q_point = 0;
1476 q_point < n_edge_quadrature_points;
1479 const double weight = edge_quadrature.weight(q_point);
1481 i, edge_quadrature_points[q_point][0], j);
1483 edge_quadrature_points[q_point][0], i, j);
1485 i, j, edge_quadrature_points[q_point][0]);
1487 if (edge_quadrature_points[q_point][0] < 0.5)
1490 i, 2.0 * edge_quadrature_points[q_point][0], j);
1493 weight * (2.0 * this->shape_value_component(
1494 dof, quadrature_point_3, 1) -
1495 this->restriction[index][i + 4 * j](
1496 (i + 4 * j) * this->degree, dof) *
1497 this->shape_value_component(
1498 (i + 4 * j) * this->degree,
1503 this->restriction[index][i + 4 * j + 2](
1504 (i + 4 * j) * this->degree, dof) *
1505 this->shape_value_component((i + 4 * j) *
1510 2.0 * edge_quadrature_points[q_point][0], i, j);
1513 (2.0 * this->shape_value_component(
1514 dof, quadrature_point_3, 0) -
1515 this->restriction[index][2 * (i + 2 * j)](
1516 (i + 4 * j + 2) * this->degree, dof) *
1517 this->shape_value_component(
1518 (i + 4 * j + 2) * this->degree,
1523 this->restriction[index][2 * (i + 2 * j) + 1](
1524 (i + 4 * j + 2) * this->degree, dof) *
1525 this->shape_value_component((i + 4 * j + 2) *
1530 i, j, 2.0 * edge_quadrature_points[q_point][0]);
1533 (2.0 * this->shape_value_component(
1534 dof, quadrature_point_3, 2) -
1535 this->restriction[index][i + 2 * j](
1536 (i + 2 * (j + 4)) * this->degree, dof) *
1537 this->shape_value_component(
1538 (i + 2 * (j + 4)) * this->degree,
1543 this->restriction[index][i + 2 * (j + 2)](
1544 (i + 2 * (j + 4)) * this->degree, dof) *
1545 this->shape_value_component((i + 2 * (j + 4)) *
1555 this->restriction[index][i + 4 * j](
1556 (i + 4 * j) * this->degree, dof) *
1557 this->shape_value_component((i + 4 * j) *
1564 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1568 (2.0 * this->shape_value_component(
1569 dof, quadrature_point_3, 1) -
1570 this->restriction[index][i + 4 * j + 2](
1571 (i + 4 * j) * this->degree, dof) *
1572 this->shape_value_component(
1573 (i + 4 * j) * this->degree,
1578 this->restriction[index][2 * (i + 2 * j)](
1579 (i + 4 * j + 2) * this->degree, dof) *
1580 this->shape_value_component((i + 4 * j + 2) *
1585 2.0 * edge_quadrature_points[q_point][0] - 1.0,
1590 (2.0 * this->shape_value_component(
1591 dof, quadrature_point_3, 0) -
1592 this->restriction[index][2 * (i + 2 * j) + 1](
1593 (i + 4 * j + 2) * this->degree, dof) *
1594 this->shape_value_component(
1595 (i + 4 * j + 2) * this->degree,
1600 this->restriction[index][i + 2 * j](
1601 (i + 2 * (j + 4)) * this->degree, dof) *
1602 this->shape_value_component((i + 2 * (j + 4)) *
1609 2.0 * edge_quadrature_points[q_point][0] - 1.0);
1612 (2.0 * this->shape_value_component(
1613 dof, quadrature_point_3, 2) -
1614 this->restriction[index][i + 2 * (j + 2)](
1615 (i + 2 * (j + 4)) * this->degree, dof) *
1616 this->shape_value_component(
1617 (i + 2 * (j + 4)) * this->degree,
1622 for (
unsigned int k = 0; k < deg; ++k)
1625 legendre_polynomials[k + 1].value(
1626 edge_quadrature_points[q_point][0]);
1628 for (
unsigned int l = 0; l < tmp.
size(); ++l)
1629 system_rhs(k, l) += tmp(l) * L_k;
1633 system_matrix_inv.
mmult(solution, system_rhs);
1635 for (
unsigned int k = 0; k < 2; ++k)
1636 for (
unsigned int l = 0; l < deg; ++l)
1638 if (
std::abs(solution(l, k)) > 1e-14)
1639 this->restriction[index][i + 2 * (2 * j + k)](
1640 (i + 4 * j) * this->degree + l + 1, dof) =
1643 if (
std::abs(solution(l, k + 2)) > 1e-14)
1644 this->restriction[index][2 * (i + 2 * j) + k](
1645 (i + 4 * j + 2) * this->degree + l + 1, dof) =
1648 if (
std::abs(solution(l, k + 4)) > 1e-14)
1649 this->restriction[index][i + 2 * (j + 2 * k)](
1650 (i + 2 * (j + 4)) * this->degree + l + 1, dof) =
1655 const QGauss<2> face_quadrature(2 * this->degree);
1656 const std::vector<Point<2>> &face_quadrature_points =
1657 face_quadrature.get_points();
1658 const std::vector<Polynomials::Polynomial<double>>
1659 &lobatto_polynomials =
1661 const unsigned int n_edge_dofs =
1663 const unsigned int n_face_quadrature_points =
1664 face_quadrature.size();
1668 n_face_quadrature_points);
1670 for (
unsigned int q_point = 0;
1671 q_point < n_face_quadrature_points;
1674 const double weight =
1675 std::sqrt(face_quadrature.weight(q_point));
1677 for (
unsigned int i = 0; i <= deg; ++i)
1680 weight * legendre_polynomials[i].value(
1681 face_quadrature_points[q_point][0]);
1683 for (
unsigned int j = 0; j < deg; ++j)
1684 assembling_matrix(i * deg + j, q_point) =
1685 L_i * lobatto_polynomials[j + 2].value(
1686 face_quadrature_points[q_point][1]);
1691 assembling_matrix.
m());
1693 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
1694 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
1695 system_matrix_inv.
invert(system_matrix);
1698 solution.reinit(system_matrix_inv.
m(), 24);
1699 system_rhs.reinit(system_matrix_inv.
m(), 24);
1702 for (
unsigned int i = 0; i < 2; ++i)
1703 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
1707 for (
unsigned int q_point = 0;
1708 q_point < n_face_quadrature_points;
1713 if (face_quadrature_points[q_point][0] < 0.5)
1715 if (face_quadrature_points[q_point][1] < 0.5)
1719 2.0 * face_quadrature_points[q_point][0],
1720 2.0 * face_quadrature_points[q_point][1]);
1722 tmp(0) += 2.0 * this->shape_value_component(
1723 dof, quadrature_point_0, 1);
1724 tmp(1) += 2.0 * this->shape_value_component(
1725 dof, quadrature_point_0, 2);
1727 2.0 * face_quadrature_points[q_point][0],
1729 2.0 * face_quadrature_points[q_point][1]);
1730 tmp(8) += 2.0 * this->shape_value_component(
1731 dof, quadrature_point_0, 2);
1732 tmp(9) += 2.0 * this->shape_value_component(
1733 dof, quadrature_point_0, 0);
1735 2.0 * face_quadrature_points[q_point][0],
1736 2.0 * face_quadrature_points[q_point][1],
1738 tmp(16) += 2.0 * this->shape_value_component(
1739 dof, quadrature_point_0, 0);
1740 tmp(17) += 2.0 * this->shape_value_component(
1741 dof, quadrature_point_0, 1);
1748 2.0 * face_quadrature_points[q_point][0],
1749 2.0 * face_quadrature_points[q_point][1] -
1752 tmp(2) += 2.0 * this->shape_value_component(
1753 dof, quadrature_point_0, 1);
1754 tmp(3) += 2.0 * this->shape_value_component(
1755 dof, quadrature_point_0, 2);
1757 2.0 * face_quadrature_points[q_point][0],
1759 2.0 * face_quadrature_points[q_point][1] -
1761 tmp(10) += 2.0 * this->shape_value_component(
1762 dof, quadrature_point_0, 2);
1763 tmp(11) += 2.0 * this->shape_value_component(
1764 dof, quadrature_point_0, 0);
1766 2.0 * face_quadrature_points[q_point][0],
1767 2.0 * face_quadrature_points[q_point][1] -
1770 tmp(18) += 2.0 * this->shape_value_component(
1771 dof, quadrature_point_0, 0);
1772 tmp(19) += 2.0 * this->shape_value_component(
1773 dof, quadrature_point_0, 1);
1777 else if (face_quadrature_points[q_point][1] < 0.5)
1781 2.0 * face_quadrature_points[q_point][0] - 1.0,
1782 2.0 * face_quadrature_points[q_point][1]);
1784 tmp(4) += 2.0 * this->shape_value_component(
1785 dof, quadrature_point_0, 1);
1786 tmp(5) += 2.0 * this->shape_value_component(
1787 dof, quadrature_point_0, 2);
1789 2.0 * face_quadrature_points[q_point][0] - 1.0,
1791 2.0 * face_quadrature_points[q_point][1]);
1792 tmp(12) += 2.0 * this->shape_value_component(
1793 dof, quadrature_point_0, 2);
1794 tmp(13) += 2.0 * this->shape_value_component(
1795 dof, quadrature_point_0, 0);
1797 2.0 * face_quadrature_points[q_point][0] - 1.0,
1798 2.0 * face_quadrature_points[q_point][1],
1800 tmp(20) += 2.0 * this->shape_value_component(
1801 dof, quadrature_point_0, 0);
1802 tmp(21) += 2.0 * this->shape_value_component(
1803 dof, quadrature_point_0, 1);
1810 2.0 * face_quadrature_points[q_point][0] - 1.0,
1811 2.0 * face_quadrature_points[q_point][1] - 1.0);
1813 tmp(6) += 2.0 * this->shape_value_component(
1814 dof, quadrature_point_0, 1);
1815 tmp(7) += 2.0 * this->shape_value_component(
1816 dof, quadrature_point_0, 2);
1818 2.0 * face_quadrature_points[q_point][0] - 1.0,
1820 2.0 * face_quadrature_points[q_point][1] - 1.0);
1821 tmp(14) += 2.0 * this->shape_value_component(
1822 dof, quadrature_point_0, 2);
1823 tmp(15) += 2.0 * this->shape_value_component(
1824 dof, quadrature_point_0, 0);
1826 2.0 * face_quadrature_points[q_point][0] - 1.0,
1827 2.0 * face_quadrature_points[q_point][1] - 1.0,
1829 tmp(22) += 2.0 * this->shape_value_component(
1830 dof, quadrature_point_0, 0);
1831 tmp(23) += 2.0 * this->shape_value_component(
1832 dof, quadrature_point_0, 1);
1837 face_quadrature_points[q_point][0],
1838 face_quadrature_points[q_point][1]);
1840 face_quadrature_points[q_point][0],
1842 face_quadrature_points[q_point][1]);
1844 face_quadrature_points[q_point][0],
1845 face_quadrature_points[q_point][1],
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(2 * (j + 2 * k)) -=
1853 this->restriction[index][i + 2 * (2 * j + k)](
1854 (i + 4 * j) * this->degree + l, dof) *
1855 this->shape_value_component(
1856 (i + 4 * j) * this->degree + l,
1859 tmp(2 * (j + 2 * k) + 1) -=
1860 this->restriction[index][i + 2 * (2 * j + k)](
1861 (i + 2 * (k + 4)) * this->degree + l, dof) *
1862 this->shape_value_component(
1863 (i + 2 * (k + 4)) * this->degree + l,
1866 tmp(2 * (j + 2 * (k + 2))) -=
1867 this->restriction[index][2 * (i + 2 * j) + k](
1868 (2 * (i + 4) + k) * this->degree + l, dof) *
1869 this->shape_value_component(
1870 (2 * (i + 4) + k) * this->degree + l,
1873 tmp(2 * (j + 2 * k) + 9) -=
1874 this->restriction[index][2 * (i + 2 * j) + k](
1875 (i + 4 * j + 2) * this->degree + l, dof) *
1876 this->shape_value_component(
1877 (i + 4 * j + 2) * this->degree + l,
1880 tmp(2 * (j + 2 * (k + 4))) -=
1881 this->restriction[index][2 * (2 * i + j) + k](
1882 (4 * i + j + 2) * this->degree + l, dof) *
1883 this->shape_value_component(
1884 (4 * i + j + 2) * this->degree + l,
1887 tmp(2 * (j + 2 * k) + 17) -=
1888 this->restriction[index][2 * (2 * i + j) + k](
1889 (4 * i + k) * this->degree + l, dof) *
1890 this->shape_value_component(
1891 (4 * i + k) * this->degree + l,
1896 tmp *= face_quadrature.weight(q_point);
1898 for (
unsigned int j = 0; j <= deg; ++j)
1900 const double L_j_0 = legendre_polynomials[j].value(
1901 face_quadrature_points[q_point][0]);
1902 const double L_j_1 = legendre_polynomials[j].value(
1903 face_quadrature_points[q_point][1]);
1905 for (
unsigned int k = 0; k < deg; ++k)
1907 const double l_k_0 =
1908 L_j_0 * lobatto_polynomials[k + 2].value(
1909 face_quadrature_points[q_point][1]);
1910 const double l_k_1 =
1911 L_j_1 * lobatto_polynomials[k + 2].value(
1912 face_quadrature_points[q_point][0]);
1914 for (
unsigned int l = 0; l < 4; ++l)
1916 system_rhs(j * deg + k, 2 * l) +=
1918 system_rhs(j * deg + k, 2 * l + 1) +=
1919 tmp(2 * l + 1) * l_k_1;
1920 system_rhs(j * deg + k, 2 * (l + 4)) +=
1921 tmp(2 * (l + 4)) * l_k_1;
1922 system_rhs(j * deg + k, 2 * l + 9) +=
1923 tmp(2 * l + 9) * l_k_0;
1924 system_rhs(j * deg + k, 2 * (l + 8)) +=
1925 tmp(2 * (l + 8)) * l_k_0;
1926 system_rhs(j * deg + k, 2 * l + 17) +=
1927 tmp(2 * l + 17) * l_k_1;
1933 system_matrix_inv.
mmult(solution, system_rhs);
1935 for (
unsigned int j = 0; j < 2; ++j)
1936 for (
unsigned int k = 0; k < 2; ++k)
1937 for (
unsigned int l = 0; l <= deg; ++l)
1938 for (
unsigned int m = 0; m < deg; ++m)
1941 2 * (j + 2 * k))) > 1e-14)
1942 this->restriction[index][i + 2 * (2 * j + k)](
1943 (2 * i * this->degree + l) * deg + m +
1945 dof) = solution(l * deg + m, 2 * (j + 2 * k));
1948 2 * (j + 2 * k) + 1)) >
1950 this->restriction[index][i + 2 * (2 * j + k)](
1951 ((2 * i + 1) * deg + m) * this->degree + l +
1954 solution(l * deg + m, 2 * (j + 2 * k) + 1);
1957 2 * (j + 2 * (k + 2)))) >
1959 this->restriction[index][2 * (i + 2 * j) + k](
1960 (2 * (i + 2) * this->degree + l) * deg + m +
1963 solution(l * deg + m, 2 * (j + 2 * (k + 2)));
1966 2 * (j + 2 * k) + 9)) >
1968 this->restriction[index][2 * (i + 2 * j) + k](
1969 ((2 * i + 5) * deg + m) * this->degree + l +
1972 solution(l * deg + m, 2 * (j + 2 * k) + 9);
1975 2 * (j + 2 * (k + 4)))) >
1977 this->restriction[index][2 * (2 * i + j) + k](
1978 (2 * (i + 4) * this->degree + l) * deg + m +
1981 solution(l * deg + m, 2 * (j + 2 * (k + 4)));
1984 2 * (j + 2 * k) + 17)) >
1986 this->restriction[index][2 * (2 * i + j) + k](
1987 ((2 * i + 9) * deg + m) * this->degree + l +
1990 solution(l * deg + m, 2 * (j + 2 * k) + 17);
1995 const std::vector<Point<dim>> &quadrature_points =
1996 quadrature.get_points();
1997 const unsigned int n_boundary_dofs =
2000 const unsigned int n_quadrature_points = quadrature.size();
2004 n_quadrature_points);
2006 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
2009 const double weight =
std::sqrt(quadrature.weight(q_point));
2011 for (
unsigned int i = 0; i <= deg; ++i)
2014 weight * legendre_polynomials[i].value(
2015 quadrature_points[q_point][0]);
2017 for (
unsigned int j = 0; j < deg; ++j)
2020 L_i * lobatto_polynomials[j + 2].value(
2021 quadrature_points[q_point][1]);
2023 for (
unsigned int k = 0; k < deg; ++k)
2024 assembling_matrix((i * deg + j) * deg + k,
2026 l_j * lobatto_polynomials[k + 2].value(
2027 quadrature_points[q_point][2]);
2033 assembling_matrix.
m());
2035 assembling_matrix.
mTmult(system_matrix, assembling_matrix);
2036 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
2037 system_matrix_inv.
invert(system_matrix);
2040 solution.reinit(system_matrix_inv.
m(), 24);
2041 system_rhs.reinit(system_matrix_inv.
m(), 24);
2044 for (
unsigned int dof = 0; dof < this->n_dofs_per_cell(); ++dof)
2048 for (
unsigned int q_point = 0; q_point < n_quadrature_points;
2053 if (quadrature_points[q_point][0] < 0.5)
2055 if (quadrature_points[q_point][1] < 0.5)
2057 if (quadrature_points[q_point][2] < 0.5)
2060 2.0 * quadrature_points[q_point][0],
2061 2.0 * quadrature_points[q_point][1],
2062 2.0 * quadrature_points[q_point][2]);
2064 tmp(0) += 2.0 * this->shape_value_component(
2065 dof, quadrature_point, 0);
2066 tmp(1) += 2.0 * this->shape_value_component(
2067 dof, quadrature_point, 1);
2068 tmp(2) += 2.0 * this->shape_value_component(
2069 dof, quadrature_point, 2);
2075 2.0 * quadrature_points[q_point][0],
2076 2.0 * quadrature_points[q_point][1],
2077 2.0 * quadrature_points[q_point][2] - 1.0);
2079 tmp(3) += 2.0 * this->shape_value_component(
2080 dof, quadrature_point, 0);
2081 tmp(4) += 2.0 * this->shape_value_component(
2082 dof, quadrature_point, 1);
2083 tmp(5) += 2.0 * this->shape_value_component(
2084 dof, quadrature_point, 2);
2088 else if (quadrature_points[q_point][2] < 0.5)
2091 2.0 * quadrature_points[q_point][0],
2092 2.0 * quadrature_points[q_point][1] - 1.0,
2093 2.0 * quadrature_points[q_point][2]);
2095 tmp(6) += 2.0 * this->shape_value_component(
2096 dof, quadrature_point, 0);
2097 tmp(7) += 2.0 * this->shape_value_component(
2098 dof, quadrature_point, 1);
2099 tmp(8) += 2.0 * this->shape_value_component(
2100 dof, quadrature_point, 2);
2106 2.0 * quadrature_points[q_point][0],
2107 2.0 * quadrature_points[q_point][1] - 1.0,
2108 2.0 * quadrature_points[q_point][2] - 1.0);
2110 tmp(9) += 2.0 * this->shape_value_component(
2111 dof, quadrature_point, 0);
2112 tmp(10) += 2.0 * this->shape_value_component(
2113 dof, quadrature_point, 1);
2114 tmp(11) += 2.0 * this->shape_value_component(
2115 dof, quadrature_point, 2);
2119 else if (quadrature_points[q_point][1] < 0.5)
2121 if (quadrature_points[q_point][2] < 0.5)
2124 2.0 * quadrature_points[q_point][0] - 1.0,
2125 2.0 * quadrature_points[q_point][1],
2126 2.0 * quadrature_points[q_point][2]);
2128 tmp(12) += 2.0 * this->shape_value_component(
2129 dof, quadrature_point, 0);
2130 tmp(13) += 2.0 * this->shape_value_component(
2131 dof, quadrature_point, 1);
2132 tmp(14) += 2.0 * this->shape_value_component(
2133 dof, quadrature_point, 2);
2139 2.0 * quadrature_points[q_point][0] - 1.0,
2140 2.0 * quadrature_points[q_point][1],
2141 2.0 * quadrature_points[q_point][2] - 1.0);
2143 tmp(15) += 2.0 * this->shape_value_component(
2144 dof, quadrature_point, 0);
2145 tmp(16) += 2.0 * this->shape_value_component(
2146 dof, quadrature_point, 1);
2147 tmp(17) += 2.0 * this->shape_value_component(
2148 dof, quadrature_point, 2);
2152 else if (quadrature_points[q_point][2] < 0.5)
2155 2.0 * quadrature_points[q_point][0] - 1.0,
2156 2.0 * quadrature_points[q_point][1] - 1.0,
2157 2.0 * quadrature_points[q_point][2]);
2160 2.0 * this->shape_value_component(dof,
2164 2.0 * this->shape_value_component(dof,
2168 2.0 * this->shape_value_component(dof,
2176 2.0 * quadrature_points[q_point][0] - 1.0,
2177 2.0 * quadrature_points[q_point][1] - 1.0,
2178 2.0 * quadrature_points[q_point][2] - 1.0);
2181 2.0 * this->shape_value_component(dof,
2185 2.0 * this->shape_value_component(dof,
2189 2.0 * this->shape_value_component(dof,
2194 for (
unsigned int i = 0; i < 2; ++i)
2195 for (
unsigned int j = 0; j < 2; ++j)
2196 for (
unsigned int k = 0; k < 2; ++k)
2197 for (
unsigned int l = 0; l <= deg; ++l)
2199 tmp(3 * (i + 2 * (j + 2 * k))) -=
2200 this->restriction[index][2 * (2 * i + j) + k](
2201 (4 * i + j + 2) * this->degree + l, dof) *
2202 this->shape_value_component(
2203 (4 * i + j + 2) * this->degree + l,
2204 quadrature_points[q_point],
2206 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2207 this->restriction[index][2 * (2 * i + j) + k](
2208 (4 * i + k) * this->degree + l, dof) *
2209 this->shape_value_component(
2210 (4 * i + k) * this->degree + l,
2211 quadrature_points[q_point],
2213 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2214 this->restriction[index][2 * (2 * i + j) + k](
2215 (2 * (j + 4) + k) * this->degree + l, dof) *
2216 this->shape_value_component(
2217 (2 * (j + 4) + k) * this->degree + l,
2218 quadrature_points[q_point],
2221 for (
unsigned int m = 0; m < deg; ++m)
2223 tmp(3 * (i + 2 * (j + 2 * k))) -=
2224 this->restriction[index][2 * (2 * i + j) +
2226 ((2 * j + 5) * deg + m) * this->degree +
2229 this->shape_value_component(
2230 ((2 * j + 5) * deg + m) * this->degree +
2232 quadrature_points[q_point],
2234 tmp(3 * (i + 2 * (j + 2 * k))) -=
2235 this->restriction[index][2 * (2 * i + j) +
2237 (2 * (i + 4) * this->degree + l) * deg +
2240 this->shape_value_component(
2241 (2 * (i + 4) * this->degree + l) * deg +
2243 quadrature_points[q_point],
2245 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2246 this->restriction[index][2 * (2 * i + j) +
2248 (2 * k * this->degree + l) * deg + m +
2251 this->shape_value_component(
2252 (2 * k * this->degree + l) * deg + m +
2254 quadrature_points[q_point],
2256 tmp(3 * (i + 2 * (j + 2 * k)) + 1) -=
2257 this->restriction[index][2 * (2 * i + j) +
2259 ((2 * i + 9) * deg + m) * this->degree +
2262 this->shape_value_component(
2263 ((2 * i + 9) * deg + m) * this->degree +
2265 quadrature_points[q_point],
2267 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2268 this->restriction[index][2 * (2 * i + j) +
2270 ((2 * k + 1) * deg + m) * this->degree +
2273 this->shape_value_component(
2274 ((2 * k + 1) * deg + m) * this->degree +
2276 quadrature_points[q_point],
2278 tmp(3 * (i + 2 * (j + 2 * k)) + 2) -=
2279 this->restriction[index][2 * (2 * i + j) +
2281 (2 * (j + 2) * this->degree + l) * deg +
2284 this->shape_value_component(
2285 (2 * (j + 2) * this->degree + l) * deg +
2287 quadrature_points[q_point],
2292 tmp *= quadrature.weight(q_point);
2294 for (
unsigned int i = 0; i <= deg; ++i)
2296 const double L_i_0 = legendre_polynomials[i].value(
2297 quadrature_points[q_point][0]);
2298 const double L_i_1 = legendre_polynomials[i].value(
2299 quadrature_points[q_point][1]);
2300 const double L_i_2 = legendre_polynomials[i].value(
2301 quadrature_points[q_point][2]);
2303 for (
unsigned int j = 0; j < deg; ++j)
2305 const double l_j_0 =
2306 L_i_0 * lobatto_polynomials[j + 2].value(
2307 quadrature_points[q_point][1]);
2308 const double l_j_1 =
2309 L_i_1 * lobatto_polynomials[j + 2].value(
2310 quadrature_points[q_point][0]);
2311 const double l_j_2 =
2312 L_i_2 * lobatto_polynomials[j + 2].value(
2313 quadrature_points[q_point][0]);
2315 for (
unsigned int k = 0; k < deg; ++k)
2317 const double l_k_0 =
2318 l_j_0 * lobatto_polynomials[k + 2].value(
2319 quadrature_points[q_point][2]);
2320 const double l_k_1 =
2321 l_j_1 * lobatto_polynomials[k + 2].value(
2322 quadrature_points[q_point][2]);
2323 const double l_k_2 =
2324 l_j_2 * lobatto_polynomials[k + 2].value(
2325 quadrature_points[q_point][1]);
2327 for (
unsigned int l = 0; l < 8; ++l)
2329 system_rhs((i * deg + j) * deg + k,
2330 3 * l) += tmp(3 * l) * l_k_0;
2331 system_rhs((i * deg + j) * deg + k,
2333 tmp(3 * l + 1) * l_k_1;
2334 system_rhs((i * deg + j) * deg + k,
2336 tmp(3 * l + 2) * l_k_2;
2343 system_matrix_inv.
mmult(solution, system_rhs);
2345 for (
unsigned int i = 0; i < 2; ++i)
2346 for (
unsigned int j = 0; j < 2; ++j)
2347 for (
unsigned int k = 0; k < 2; ++k)
2348 for (
unsigned int l = 0; l <= deg; ++l)
2349 for (
unsigned int m = 0; m < deg; ++m)
2350 for (
unsigned int n = 0; n < deg; ++n)
2353 solution((l * deg + m) * deg + n,
2354 3 * (i + 2 * (j + 2 * k)))) >
2356 this->restriction[index][2 * (2 * i + j) + k](
2357 (l * deg + m) * deg + n + n_boundary_dofs,
2358 dof) = solution((l * deg + m) * deg + n,
2359 3 * (i + 2 * (j + 2 * k)));
2362 solution((l * deg + m) * deg + n,
2363 3 * (i + 2 * (j + 2 * k)) + 1)) >
2365 this->restriction[index][2 * (2 * i + j) + k](
2366 (l + (m + deg) * this->degree) * deg + n +
2369 solution((l * deg + m) * deg + n,
2370 3 * (i + 2 * (j + 2 * k)) + 1);
2373 solution((l * deg + m) * deg + n,
2374 3 * (i + 2 * (j + 2 * k)) + 2)) >
2376 this->restriction[index][2 * (2 * i + j) + k](
2378 ((m + 2 * deg) * deg + n) * this->degree +
2381 solution((l * deg + m) * deg + n,
2382 3 * (i + 2 * (j + 2 * k)) + 2);
3492 std::vector<double> &nodal_values)
const
3497 const unsigned int face_no = 0;
3499 const unsigned int deg = this->degree - 1;
3500 Assert(support_point_values.size() == this->generalized_support_points.size(),
3502 this->generalized_support_points.size()));
3503 Assert(support_point_values[0].
size() == this->n_components(),
3505 this->n_components()));
3506 Assert(nodal_values.size() == this->n_dofs_per_cell(),
3508 std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
3516 const QGauss<1> reference_edge_quadrature(this->degree);
3517 const unsigned int n_edge_points = reference_edge_quadrature.size();
3519 for (
unsigned int i = 0; i < 2; ++i)
3520 for (
unsigned int j = 0; j < 2; ++j)
3522 for (
unsigned int q_point = 0; q_point < n_edge_points;
3524 nodal_values[(i + 2 * j) * this->degree] +=
3525 reference_edge_quadrature.weight(q_point) *
3526 support_point_values[q_point + (i + 2 * j) * n_edge_points]
3532 if (
std::abs(nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
3533 nodal_values[(i + 2 * j) * this->degree] = 0.0;
3546 if (this->degree > 1)
3551 const std::vector<Polynomials::Polynomial<double>>
3552 &lobatto_polynomials =
3556 std::vector<Polynomials::Polynomial<double>>
3557 lobatto_polynomials_grad(this->degree);
3559 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3560 lobatto_polynomials_grad[i] =
3561 lobatto_polynomials[i + 1].derivative();
3566 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3567 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3568 for (
unsigned int q_point = 0; q_point < n_edge_points;
3570 system_matrix(i, j) +=
3571 boundary_weights(q_point, j) *
3572 lobatto_polynomials_grad[i + 1].value(
3573 this->generalized_face_support_points[face_no][q_point]
3579 system_matrix_inv.
invert(system_matrix);
3586 for (
unsigned int line = 0;
3587 line < GeometryInfo<dim>::lines_per_cell;
3593 for (
unsigned int q_point = 0; q_point < n_edge_points;
3597 support_point_values[line * n_edge_points + q_point]
3598 [line_coordinate[line]] -
3599 nodal_values[line * this->degree] *
3600 this->shape_value_component(
3601 line * this->degree,
3602 this->generalized_support_points[line *
3605 line_coordinate[line]);
3607 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3608 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3611 system_matrix_inv.
vmult(solution, system_rhs);
3617 for (
unsigned int i = 0; i < solution.
size(); ++i)
3619 nodal_values[line * this->degree + i + 1] = solution(i);
3631 const QGauss<dim> reference_quadrature(this->degree);
3632 const unsigned int n_interior_points =
3633 reference_quadrature.size();
3634 const std::vector<Polynomials::Polynomial<double>>
3635 &legendre_polynomials =
3639 system_matrix.reinit((this->degree - 1) * this->degree,
3640 (this->degree - 1) * this->degree);
3643 for (
unsigned int i = 0; i < this->degree; ++i)
3644 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3645 for (
unsigned int k = 0; k < this->degree; ++k)
3646 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3647 for (
unsigned int q_point = 0;
3648 q_point < n_interior_points;
3650 system_matrix(i * (this->degree - 1) + j,
3651 k * (this->degree - 1) + l) +=
3652 reference_quadrature.weight(q_point) *
3653 legendre_polynomials[i].value(
3654 this->generalized_support_points
3656 n_edge_points][0]) *
3657 lobatto_polynomials[j + 2].value(
3658 this->generalized_support_points
3660 n_edge_points][1]) *
3661 lobatto_polynomials_grad[k].value(
3662 this->generalized_support_points
3664 n_edge_points][0]) *
3665 lobatto_polynomials[l + 2].value(
3666 this->generalized_support_points
3670 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3671 system_matrix_inv.
invert(system_matrix);
3675 system_rhs.
reinit(system_matrix_inv.
m());
3678 for (
unsigned int q_point = 0; q_point < n_interior_points;
3682 support_point_values[q_point +
3686 for (
unsigned int i = 0; i < 2; ++i)
3687 for (
unsigned int j = 0; j <= deg; ++j)
3688 tmp -= nodal_values[(i + 2) * this->degree + j] *
3689 this->shape_value_component(
3690 (i + 2) * this->degree + j,
3691 this->generalized_support_points
3696 for (
unsigned int i = 0; i <= deg; ++i)
3697 for (
unsigned int j = 0; j < deg; ++j)
3698 system_rhs(i * deg + j) +=
3699 reference_quadrature.weight(q_point) * tmp *
3700 lobatto_polynomials_grad[i].value(
3701 this->generalized_support_points
3703 n_edge_points][0]) *
3704 lobatto_polynomials[j + 2].value(
3705 this->generalized_support_points
3710 solution.
reinit(system_matrix.
m());
3711 system_matrix_inv.
vmult(solution, system_rhs);
3717 for (
unsigned int i = 0; i <= deg; ++i)
3718 for (
unsigned int j = 0; j < deg; ++j)
3719 if (
std::abs(solution(i * deg + j)) > 1e-14)
3722 solution(i * deg + j);
3729 for (
unsigned int q_point = 0; q_point < n_interior_points;
3733 support_point_values[q_point +
3737 for (
unsigned int i = 0; i < 2; ++i)
3738 for (
unsigned int j = 0; j <= deg; ++j)
3739 tmp -= nodal_values[i * this->degree + j] *
3740 this->shape_value_component(
3741 i * this->degree + j,
3742 this->generalized_support_points
3747 for (
unsigned int i = 0; i <= deg; ++i)
3748 for (
unsigned int j = 0; j < deg; ++j)
3749 system_rhs(i * deg + j) +=
3750 reference_quadrature.weight(q_point) * tmp *
3751 lobatto_polynomials_grad[i].value(
3752 this->generalized_support_points
3754 n_edge_points][1]) *
3755 lobatto_polynomials[j + 2].value(
3756 this->generalized_support_points
3761 system_matrix_inv.
vmult(solution, system_rhs);
3767 for (
unsigned int i = 0; i <= deg; ++i)
3768 for (
unsigned int j = 0; j < deg; ++j)
3769 if (
std::abs(solution(i * deg + j)) > 1e-14)
3772 this->degree] = solution(i * deg + j);
3782 const QGauss<1> reference_edge_quadrature(this->degree);
3783 const unsigned int n_edge_points = reference_edge_quadrature.size();
3785 for (
unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
3787 for (
unsigned int i = 0; i < 4; ++i)
3788 nodal_values[(i + 8) * this->degree] +=
3789 reference_edge_quadrature.weight(q_point) *
3790 support_point_values[q_point + (i + 8) * n_edge_points][2];
3792 for (
unsigned int i = 0; i < 2; ++i)
3793 for (
unsigned int j = 0; j < 2; ++j)
3794 for (
unsigned int k = 0; k < 2; ++k)
3795 nodal_values[(i + 2 * (2 * j + k)) * this->degree] +=
3796 reference_edge_quadrature.weight(q_point) *
3797 support_point_values[q_point + (i + 2 * (2 * j + k)) *
3798 n_edge_points][1 - k];
3805 for (
unsigned int i = 0; i < 4; ++i)
3806 if (
std::abs(nodal_values[(i + 8) * this->degree]) < 1e-14)
3807 nodal_values[(i + 8) * this->degree] = 0.0;
3809 for (
unsigned int i = 0; i < 2; ++i)
3810 for (
unsigned int j = 0; j < 2; ++j)
3811 for (
unsigned int k = 0; k < 2; ++k)
3813 nodal_values[(i + 2 * (2 * j + k)) * this->degree]) <
3815 nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
3826 if (this->degree > 1)
3831 const std::vector<Polynomials::Polynomial<double>>
3832 &lobatto_polynomials =
3836 std::vector<Polynomials::Polynomial<double>>
3837 lobatto_polynomials_grad(this->degree);
3839 for (
unsigned int i = 0; i < lobatto_polynomials_grad.size(); ++i)
3840 lobatto_polynomials_grad[i] =
3841 lobatto_polynomials[i + 1].derivative();
3846 for (
unsigned int i = 0; i < system_matrix.
m(); ++i)
3847 for (
unsigned int j = 0; j < system_matrix.
n(); ++j)
3848 for (
unsigned int q_point = 0; q_point < n_edge_points;
3850 system_matrix(i, j) +=
3851 boundary_weights(q_point, j) *
3852 lobatto_polynomials_grad[i + 1].value(
3853 this->generalized_face_support_points[face_no][q_point]
3859 system_matrix_inv.
invert(system_matrix);
3863 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
3867 for (
unsigned int line = 0;
3868 line < GeometryInfo<dim>::lines_per_cell;
3874 for (
unsigned int q_point = 0; q_point < this->degree;
3878 support_point_values[line * this->degree + q_point]
3879 [line_coordinate[line]] -
3880 nodal_values[line * this->degree] *
3881 this->shape_value_component(
3882 line * this->degree,
3884 ->generalized_support_points[line * this->degree +
3886 line_coordinate[line]);
3888 for (
unsigned int i = 0; i < system_rhs.
size(); ++i)
3889 system_rhs(i) += boundary_weights(q_point, i) * tmp;
3892 system_matrix_inv.
vmult(solution, system_rhs);
3898 for (
unsigned int i = 0; i < solution.
size(); ++i)
3900 nodal_values[line * this->degree + i + 1] = solution(i);
3911 const std::vector<Polynomials::Polynomial<double>>
3912 &legendre_polynomials =
3915 const unsigned int n_face_points = n_edge_points * n_edge_points;
3917 system_matrix.reinit((this->degree - 1) * this->degree,
3918 (this->degree - 1) * this->degree);
3921 for (
unsigned int i = 0; i < this->degree; ++i)
3922 for (
unsigned int j = 0; j < this->degree - 1; ++j)
3923 for (
unsigned int k = 0; k < this->degree; ++k)
3924 for (
unsigned int l = 0; l < this->degree - 1; ++l)
3925 for (
unsigned int q_point = 0; q_point < n_face_points;
3927 system_matrix(i * (this->degree - 1) + j,
3928 k * (this->degree - 1) + l) +=
3929 boundary_weights(q_point + n_edge_points,
3930 2 * (k * (this->degree - 1) + l)) *
3931 legendre_polynomials[i].value(
3932 this->generalized_face_support_points
3933 [face_no][q_point + 4 * n_edge_points][0]) *
3934 lobatto_polynomials[j + 2].value(
3935 this->generalized_face_support_points
3936 [face_no][q_point + 4 * n_edge_points][1]);
3938 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
3939 system_matrix_inv.
invert(system_matrix);
3940 solution.
reinit(system_matrix.
m());
3941 system_rhs.
reinit(system_matrix.
m());
3945 {1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
3962 for (
unsigned int q_point = 0; q_point < n_face_points;
3966 support_point_values[q_point +
3969 [face_coordinates[face][0]];
3971 for (
unsigned int i = 0; i < 2; ++i)
3972 for (
unsigned int j = 0; j <= deg; ++j)
3976 this->shape_value_component(
3978 this->generalized_support_points
3981 face_coordinates[face][0]);
3983 for (
unsigned int i = 0; i <= deg; ++i)
3984 for (
unsigned int j = 0; j < deg; ++j)
3985 system_rhs(i * deg + j) +=
3986 boundary_weights(q_point + n_edge_points,
3987 2 * (i * deg + j)) *
3991 system_matrix_inv.
vmult(solution, system_rhs);
3997 for (
unsigned int i = 0; i <= deg; ++i)
3998 for (
unsigned int j = 0; j < deg; ++j)
3999 if (
std::abs(solution(i * deg + j)) > 1e-14)
4000 nodal_values[(2 * face * this->degree + i +
4004 solution(i * deg + j);
4011 for (
unsigned int q_point = 0; q_point < n_face_points;
4015 support_point_values[q_point +
4018 [face_coordinates[face][1]];
4020 for (
unsigned int i = 2;
4021 i < GeometryInfo<dim>::lines_per_face;
4023 for (
unsigned int j = 0; j <= deg; ++j)
4027 this->shape_value_component(
4029 this->generalized_support_points
4032 face_coordinates[face][1]);
4034 for (
unsigned int i = 0; i <= deg; ++i)
4035 for (
unsigned int j = 0; j < deg; ++j)
4036 system_rhs(i * deg + j) +=
4037 boundary_weights(q_point + n_edge_points,
4038 2 * (i * deg + j) + 1) *
4042 system_matrix_inv.
vmult(solution, system_rhs);
4048 for (
unsigned int i = 0; i <= deg; ++i)
4049 for (
unsigned int j = 0; j < deg; ++j)
4050 if (
std::abs(solution(i * deg + j)) > 1e-14)
4051 nodal_values[((2 * face + 1) * deg + j +
4054 i] = solution(i * deg + j);
4062 const QGauss<dim> reference_quadrature(this->degree);
4063 const unsigned int n_interior_points =
4064 reference_quadrature.size();
4068 system_matrix.reinit(this->degree * deg * deg,
4069 this->degree * deg * deg);
4072 for (
unsigned int i = 0; i <= deg; ++i)
4073 for (
unsigned int j = 0; j < deg; ++j)
4074 for (
unsigned int k = 0; k < deg; ++k)
4075 for (
unsigned int l = 0; l <= deg; ++l)
4076 for (
unsigned int m = 0; m < deg; ++m)
4077 for (
unsigned int n = 0; n < deg; ++n)
4078 for (
unsigned int q_point = 0;
4079 q_point < n_interior_points;
4081 system_matrix((i * deg + j) * deg + k,
4082 (l * deg + m) * deg + n) +=
4083 reference_quadrature.weight(q_point) *
4084 legendre_polynomials[i].value(
4085 this->generalized_support_points
4090 n_face_points][0]) *
4091 lobatto_polynomials[j + 2].value(
4092 this->generalized_support_points
4097 n_face_points][1]) *
4098 lobatto_polynomials[k + 2].value(
4099 this->generalized_support_points
4104 n_face_points][2]) *
4105 lobatto_polynomials_grad[l].value(
4106 this->generalized_support_points
4111 n_face_points][0]) *
4112 lobatto_polynomials[m + 2].value(
4113 this->generalized_support_points
4118 n_face_points][1]) *
4119 lobatto_polynomials[n + 2].value(
4120 this->generalized_support_points
4127 system_matrix_inv.reinit(system_matrix.
m(), system_matrix.
m());
4128 system_matrix_inv.
invert(system_matrix);
4130 system_rhs.
reinit(system_matrix.
m());
4133 for (
unsigned int q_point = 0; q_point < n_interior_points;
4137 support_point_values[q_point +
4143 for (
unsigned int i = 0; i <= deg; ++i)
4145 for (
unsigned int j = 0; j < 2; ++j)
4146 for (
unsigned int k = 0; k < 2; ++k)
4148 nodal_values[i + (j + 4 * k + 2) * this->degree] *
4149 this->shape_value_component(
4150 i + (j + 4 * k + 2) * this->degree,
4151 this->generalized_support_points
4159 for (
unsigned int j = 0; j < deg; ++j)
4160 for (
unsigned int k = 0; k < 4; ++k)
4162 nodal_values[(i + 2 * (k + 2) * this->degree +
4167 this->shape_value_component(
4168 (i + 2 * (k + 2) * this->degree +
4172 this->generalized_support_points
4181 for (
unsigned int i = 0; i <= deg; ++i)
4182 for (
unsigned int j = 0; j < deg; ++j)
4183 for (
unsigned int k = 0; k < deg; ++k)
4184 system_rhs((i * deg + j) * deg + k) +=
4185 reference_quadrature.weight(q_point) * tmp *
4186 lobatto_polynomials_grad[i].value(
4187 this->generalized_support_points
4192 n_face_points][0]) *
4193 lobatto_polynomials[j + 2].value(
4194 this->generalized_support_points
4199 n_face_points][1]) *
4200 lobatto_polynomials[k + 2].value(
4201 this->generalized_support_points
4210 system_matrix_inv.
vmult(solution, system_rhs);
4216 for (
unsigned int i = 0; i <= deg; ++i)
4217 for (
unsigned int j = 0; j < deg; ++j)
4218 for (
unsigned int k = 0; k < deg; ++k)
4219 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4226 solution((i * deg + j) * deg + k);
4231 for (
unsigned int q_point = 0; q_point < n_interior_points;
4235 support_point_values[q_point +
4241 for (
unsigned int i = 0; i <= deg; ++i)
4242 for (
unsigned int j = 0; j < 2; ++j)
4244 for (
unsigned int k = 0; k < 2; ++k)
4245 tmp -= nodal_values[i + (4 * j + k) * this->degree] *
4246 this->shape_value_component(
4247 i + (4 * j + k) * this->degree,
4248 this->generalized_support_points
4256 for (
unsigned int k = 0; k < deg; ++k)
4258 nodal_values[(i + 2 * j * this->degree +
4263 this->shape_value_component(
4264 (i + 2 * j * this->degree +
4268 this->generalized_support_points
4276 ((2 * j + 9) * deg + k +
4279 this->shape_value_component(
4280 i + ((2 * j + 9) * deg + k +
4283 this->generalized_support_points
4292 for (
unsigned int i = 0; i <= deg; ++i)
4293 for (
unsigned int j = 0; j < deg; ++j)
4294 for (
unsigned int k = 0; k < deg; ++k)
4295 system_rhs((i * deg + j) * deg + k) +=
4296 reference_quadrature.weight(q_point) * tmp *
4297 lobatto_polynomials_grad[i].value(
4298 this->generalized_support_points
4303 n_face_points][1]) *
4304 lobatto_polynomials[j + 2].value(
4305 this->generalized_support_points
4310 n_face_points][0]) *
4311 lobatto_polynomials[k + 2].value(
4312 this->generalized_support_points
4320 system_matrix_inv.
vmult(solution, system_rhs);
4326 for (
unsigned int i = 0; i <= deg; ++i)
4327 for (
unsigned int j = 0; j < deg; ++j)
4328 for (
unsigned int k = 0; k < deg; ++k)
4329 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4330 nodal_values[((i + this->degree +
4337 solution((i * deg + j) * deg + k);
4342 for (
unsigned int q_point = 0; q_point < n_interior_points;
4346 support_point_values[q_point +
4352 for (
unsigned int i = 0; i <= deg; ++i)
4353 for (
unsigned int j = 0; j < 4; ++j)
4355 tmp -= nodal_values[i + (j + 8) * this->degree] *
4356 this->shape_value_component(
4357 i + (j + 8) * this->degree,
4358 this->generalized_support_points
4366 for (
unsigned int k = 0; k < deg; ++k)
4369 ((2 * j + 1) * deg + k +
4372 this->shape_value_component(
4373 i + ((2 * j + 1) * deg + k +
4376 this->generalized_support_points
4385 for (
unsigned int i = 0; i <= deg; ++i)
4386 for (
unsigned int j = 0; j < deg; ++j)
4387 for (
unsigned int k = 0; k < deg; ++k)
4388 system_rhs((i * deg + j) * deg + k) +=
4389 reference_quadrature.weight(q_point) * tmp *
4390 lobatto_polynomials_grad[i].value(
4391 this->generalized_support_points
4396 n_face_points][2]) *
4397 lobatto_polynomials[j + 2].value(
4398 this->generalized_support_points
4403 n_face_points][0]) *
4404 lobatto_polynomials[k + 2].value(
4405 this->generalized_support_points
4413 system_matrix_inv.
vmult(solution, system_rhs);
4419 for (
unsigned int i = 0; i <= deg; ++i)
4420 for (
unsigned int j = 0; j < deg; ++j)
4421 for (
unsigned int k = 0; k < deg; ++k)
4422 if (
std::abs(solution((i * deg + j) * deg + k)) > 1e-14)
4428 this->degree] = solution((i * deg + j) * deg + k);