235 typename ::FiniteElement<dim, spacedim>::InternalDataBase>
236 data_ptr = std::make_unique<InternalData>();
238 data.update_each = requires_update_flags(update_flags);
241 const unsigned int degree(this->degree - 1);
247 const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
250 const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
253 const unsigned int n_q_points = quadrature.
size();
256 data.sigma_imj_values.resize(
258 std::vector<std::vector<double>>(vertices_per_cell,
259 std::vector<double>(vertices_per_cell)));
261 data.sigma_imj_grads.resize(vertices_per_cell,
262 std::vector<std::vector<double>>(
263 vertices_per_cell, std::vector<double>(dim)));
268 data.shape_values.resize(this->n_dofs_per_cell(),
274 data.shape_grads.resize(this->n_dofs_per_cell(),
281 data.shape_hessians.resize(this->n_dofs_per_cell(),
286 std::vector<Point<dim>> p_list(n_q_points);
296 std::vector<std::vector<double>> sigma(
297 n_q_points, std::vector<double>(lines_per_cell));
298 std::vector<std::vector<double>> lambda(
299 n_q_points, std::vector<double>(lines_per_cell));
301 for (
unsigned int q = 0; q < n_q_points; ++q)
303 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
304 sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
305 sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
306 sigma[q][3] = p_list[q][0] + p_list[q][1];
308 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
309 lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
310 lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
311 lambda[q][3] = p_list[q][0] * p_list[q][1];
312 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
314 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
316 data.sigma_imj_values[q][i][j] =
317 sigma[q][i] - sigma[q][j];
330 {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
331 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
332 unsigned int sigma_imj_component[vertices_per_cell]
335 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
337 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
344 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
345 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
349 sigma_imj_component[i][j] = 0;
350 for (
unsigned int d = 0; d < dim; ++d)
353 sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
358 sigma_imj_component[i][j] = d;
365 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
366 2.0 * sigma_imj_sign[i][j];
374 data.edge_sigma_values.resize(lines_per_cell);
375 data.edge_sigma_grads.resize(lines_per_cell);
376 for (
unsigned int m = 0; m < lines_per_cell; ++m)
378 data.edge_sigma_values[m].resize(n_q_points);
381 data.edge_sigma_grads[m].resize(dim);
391 data.edge_lambda_values.resize(lines_per_cell,
392 std::vector<double>(n_q_points));
393 data.edge_lambda_grads_2d.resize(lines_per_cell,
394 std::vector<double>(dim));
395 for (
unsigned int m = 0; m < lines_per_cell; ++m)
400 const unsigned int e1(
402 const unsigned int e2(
404 for (
unsigned int q = 0; q < n_q_points; ++q)
406 data.edge_sigma_values[m][q] =
407 data.sigma_imj_values[q][e2][e1];
408 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
411 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
413 data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
414 data.edge_lambda_grads_2d[1] = {1.0, 0.0};
415 data.edge_lambda_grads_2d[2] = {0.0, -1.0};
416 data.edge_lambda_grads_2d[3] = {0.0, 1.0};
484 const unsigned int cell_type1_offset = n_line_dofs;
491 const unsigned int cell_type2_offset =
492 cell_type1_offset + degree * degree;
500 const unsigned int cell_type3_offset1 =
501 cell_type2_offset + degree * degree;
502 const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
507 std::vector<Point<dim>> cell_points(n_q_points);
508 for (
unsigned int q = 0; q < n_q_points; ++q)
510 for (
unsigned int d = 0; d < dim; ++d)
512 cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
517 for (
unsigned int q = 0; q < n_q_points; ++q)
529 const unsigned int poly_length =
534 std::vector<std::vector<double>> polyx(
535 degree, std::vector<double>(poly_length));
536 std::vector<std::vector<double>> polyy(
537 degree, std::vector<double>(poly_length));
538 for (
unsigned int i = 0; i < degree; ++i)
543 IntegratedLegendrePolynomials[i + 2].value(
544 cell_points[q][0], polyx[i]);
545 IntegratedLegendrePolynomials[i + 2].value(
546 cell_points[q][1], polyy[i]);
551 for (
unsigned int j = 0; j < degree; ++j)
553 const unsigned int shift_j(j * degree);
554 for (
unsigned int i = 0; i < degree; ++i)
556 const unsigned int shift_ij(i + shift_j);
559 const unsigned int dof_index1(cell_type1_offset +
561 data.shape_values[dof_index1][q][0] =
562 2.0 * polyx[i][1] * polyy[j][0];
563 data.shape_values[dof_index1][q][1] =
564 2.0 * polyx[i][0] * polyy[j][1];
567 const unsigned int dof_index2(cell_type2_offset +
569 data.shape_values[dof_index2][q][0] =
570 data.shape_values[dof_index1][q][0];
571 data.shape_values[dof_index2][q][1] =
572 -1.0 * data.shape_values[dof_index1][q][1];
575 const unsigned int dof_index3_1(cell_type3_offset1 +
577 data.shape_values[dof_index3_1][q][0] = polyy[j][0];
578 data.shape_values[dof_index3_1][q][1] = 0.0;
580 const unsigned int dof_index3_2(cell_type3_offset2 +
582 data.shape_values[dof_index3_2][q][0] = 0.0;
583 data.shape_values[dof_index3_2][q][1] = polyx[j][0];
588 for (
unsigned int j = 0; j < degree; ++j)
590 const unsigned int shift_j(j * degree);
591 for (
unsigned int i = 0; i < degree; ++i)
593 const unsigned int shift_ij(i + shift_j);
596 const unsigned int dof_index1(cell_type1_offset +
598 data.shape_grads[dof_index1][q][0][0] =
599 4.0 * polyx[i][2] * polyy[j][0];
600 data.shape_grads[dof_index1][q][0][1] =
601 4.0 * polyx[i][1] * polyy[j][1];
602 data.shape_grads[dof_index1][q][1][0] =
603 data.shape_grads[dof_index1][q][0][1];
604 data.shape_grads[dof_index1][q][1][1] =
605 4.0 * polyx[i][0] * polyy[j][2];
608 const unsigned int dof_index2(cell_type2_offset +
610 data.shape_grads[dof_index2][q][0][0] =
611 data.shape_grads[dof_index1][q][0][0];
612 data.shape_grads[dof_index2][q][0][1] =
613 data.shape_grads[dof_index1][q][0][1];
614 data.shape_grads[dof_index2][q][1][0] =
615 -1.0 * data.shape_grads[dof_index1][q][1][0];
616 data.shape_grads[dof_index2][q][1][1] =
617 -1.0 * data.shape_grads[dof_index1][q][1][1];
620 const unsigned int dof_index3_1(cell_type3_offset1 +
622 data.shape_grads[dof_index3_1][q][0][0] = 0.0;
623 data.shape_grads[dof_index3_1][q][0][1] =
625 data.shape_grads[dof_index3_1][q][1][0] = 0.0;
626 data.shape_grads[dof_index3_1][q][1][1] = 0.0;
628 const unsigned int dof_index3_2(cell_type3_offset2 +
630 data.shape_grads[dof_index3_2][q][0][0] = 0.0;
631 data.shape_grads[dof_index3_2][q][0][1] = 0.0;
632 data.shape_grads[dof_index3_2][q][1][0] =
634 data.shape_grads[dof_index3_2][q][1][1] = 0.0;
639 for (
unsigned int j = 0; j < degree; ++j)
641 const unsigned int shift_j(j * degree);
642 for (
unsigned int i = 0; i < degree; ++i)
644 const unsigned int shift_ij(i + shift_j);
647 const unsigned int dof_index1(cell_type1_offset +
649 data.shape_hessians[dof_index1][q][0][0][0] =
650 8.0 * polyx[i][3] * polyy[j][0];
651 data.shape_hessians[dof_index1][q][1][0][0] =
652 8.0 * polyx[i][2] * polyy[j][1];
654 data.shape_hessians[dof_index1][q][0][1][0] =
655 data.shape_hessians[dof_index1][q][1][0][0];
656 data.shape_hessians[dof_index1][q][1][1][0] =
657 8.0 * polyx[i][1] * polyy[j][2];
659 data.shape_hessians[dof_index1][q][0][0][1] =
660 data.shape_hessians[dof_index1][q][1][0][0];
661 data.shape_hessians[dof_index1][q][1][0][1] =
662 data.shape_hessians[dof_index1][q][1][1][0];
664 data.shape_hessians[dof_index1][q][0][1][1] =
665 data.shape_hessians[dof_index1][q][1][1][0];
666 data.shape_hessians[dof_index1][q][1][1][1] =
667 8.0 * polyx[i][0] * polyy[j][3];
672 const unsigned int dof_index2(cell_type2_offset +
674 for (
unsigned int d = 0; d < dim; ++d)
676 data.shape_hessians[dof_index2][q][0][0][d] =
677 data.shape_hessians[dof_index1][q][0][0][d];
678 data.shape_hessians[dof_index2][q][0][1][d] =
679 data.shape_hessians[dof_index1][q][0][1][d];
680 data.shape_hessians[dof_index2][q][1][0][d] =
682 data.shape_hessians[dof_index1][q][1][0][d];
683 data.shape_hessians[dof_index2][q][1][1][d] =
685 data.shape_hessians[dof_index1][q][1][1][d];
689 const unsigned int dof_index3_1(cell_type3_offset1 +
691 data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
692 data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
693 data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
694 data.shape_hessians[dof_index3_1][q][0][1][1] =
696 data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
697 data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
698 data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
699 data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
701 const unsigned int dof_index3_2(cell_type3_offset2 +
703 data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
704 data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
705 data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
706 data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
707 data.shape_hessians[dof_index3_2][q][1][0][0] =
709 data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
710 data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
711 data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
722 std::vector<std::vector<double>> sigma(
723 n_q_points, std::vector<double>(lines_per_cell));
724 std::vector<std::vector<double>> lambda(
725 n_q_points, std::vector<double>(lines_per_cell));
726 for (
unsigned int q = 0; q < n_q_points; ++q)
728 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
731 p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
733 (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
734 sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
736 (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
737 sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
738 sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
739 sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
741 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
742 (1.0 - p_list[q][2]);
744 p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
746 (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
747 lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
749 (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
750 lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
751 lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
752 lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
756 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
758 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
760 data.sigma_imj_values[q][i][j] =
761 sigma[q][i] - sigma[q][j];
794 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
795 unsigned int sigma_imj_component[vertices_per_cell]
798 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
800 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
807 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
808 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
812 sigma_imj_component[i][j] = 0;
813 for (
unsigned int d = 0; d < dim; ++d)
816 sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
821 sigma_imj_component[i][j] = d;
828 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
829 2.0 * sigma_imj_sign[i][j];
836 data.edge_sigma_values.resize(lines_per_cell);
837 data.edge_lambda_values.resize(lines_per_cell);
838 data.edge_sigma_grads.resize(lines_per_cell);
839 data.edge_lambda_grads_3d.resize(lines_per_cell);
840 data.edge_lambda_gradgrads_3d.resize(lines_per_cell);
841 for (
unsigned int m = 0; m < lines_per_cell; ++m)
843 data.edge_sigma_values[m].resize(n_q_points);
844 data.edge_lambda_values[m].resize(n_q_points);
847 data.edge_sigma_grads[m].resize(dim);
849 data.edge_lambda_grads_3d[m].resize(n_q_points);
850 for (
unsigned int q = 0; q < n_q_points; ++q)
852 data.edge_lambda_grads_3d[m][q].resize(dim);
856 data.edge_lambda_gradgrads_3d[m].resize(dim);
857 for (
unsigned int d = 0; d < dim; ++d)
859 data.edge_lambda_gradgrads_3d[m][d].resize(dim);
866 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
868 for (
unsigned int m = 0; m < lines_per_cell; ++m)
873 const unsigned int e1(
875 const unsigned int e2(
878 for (
unsigned int q = 0; q < n_q_points; ++q)
880 data.edge_sigma_values[m][q] =
881 data.sigma_imj_values[q][e2][e1];
882 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
885 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
888 for (
unsigned int q = 0; q < n_q_points; ++q)
890 double x(p_list[q][0]);
891 double y(p_list[q][1]);
892 double z(p_list[q][2]);
893 data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
894 data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
895 data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
896 data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
897 data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
898 data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
899 data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
900 data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
901 data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
902 data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
903 data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
904 data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
934 for (
unsigned int m = 0; m < lines_per_cell; ++m)
936 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
937 [edge_lambda_directions[m][1]] =
939 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
940 [edge_lambda_directions[m][0]] =
948 data.face_lambda_values.resize(faces_per_cell);
949 data.face_lambda_grads.resize(faces_per_cell);
951 for (
unsigned int m = 0; m < faces_per_cell; ++m)
953 data.face_lambda_values[m].resize(n_q_points);
954 data.face_lambda_grads[m].resize(3);
957 for (
unsigned int q = 0; q < n_q_points; ++q)
959 double x(p_list[q][0]);
960 double y(p_list[q][1]);
961 double z(p_list[q][2]);
962 data.face_lambda_values[0][q] = 1.0 - x;
963 data.face_lambda_values[1][q] = x;
964 data.face_lambda_values[2][q] = 1.0 - y;
965 data.face_lambda_values[3][q] = y;
966 data.face_lambda_values[4][q] = 1.0 - z;
967 data.face_lambda_values[5][q] = z;
970 data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
971 data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
972 data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
973 data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
974 data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
975 data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
988 const unsigned int cell_type1_offset(n_line_dofs +
1001 const unsigned int cell_type2_offset1(
1002 cell_type1_offset + degree * degree * degree);
1003 const unsigned int cell_type2_offset2(
1004 cell_type2_offset1 + degree * degree * degree);
1014 const unsigned int cell_type3_offset1(
1015 cell_type2_offset2 + degree * degree * degree);
1016 const unsigned int cell_type3_offset2(cell_type3_offset1 +
1018 const unsigned int cell_type3_offset3(cell_type3_offset2 +
1022 std::vector<Point<dim>> cell_points(n_q_points);
1023 for (
unsigned int q = 0; q < n_q_points; ++q)
1025 for (
unsigned int d = 0; d < dim; ++d)
1027 cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
1035 const unsigned int poly_length =
1041 for (
unsigned int q = 0; q < n_q_points; ++q)
1049 std::vector<std::vector<double>> polyx(
1050 degree, std::vector<double>(poly_length));
1051 std::vector<std::vector<double>> polyy(
1052 degree, std::vector<double>(poly_length));
1053 std::vector<std::vector<double>> polyz(
1054 degree, std::vector<double>(poly_length));
1055 for (
unsigned int i = 0; i < degree; ++i)
1058 IntegratedLegendrePolynomials[i + 2].value(
1059 cell_points[q][0], polyx[i]);
1060 IntegratedLegendrePolynomials[i + 2].value(
1061 cell_points[q][1], polyy[i]);
1062 IntegratedLegendrePolynomials[i + 2].value(
1063 cell_points[q][2], polyz[i]);
1068 for (
unsigned int k = 0; k < degree; ++k)
1070 const unsigned int shift_k(k * degree * degree);
1071 const unsigned int shift_j(
1074 for (
unsigned int j = 0; j < degree; ++j)
1076 const unsigned int shift_jk(j * degree +
1078 for (
unsigned int i = 0; i < degree; ++i)
1080 const unsigned int shift_ijk(shift_jk +
1084 const unsigned int dof_index1(
1085 cell_type1_offset + shift_ijk);
1087 data.shape_values[dof_index1][q][0] =
1088 2.0 * polyx[i][1] * polyy[j][0] *
1090 data.shape_values[dof_index1][q][1] =
1091 2.0 * polyx[i][0] * polyy[j][1] *
1093 data.shape_values[dof_index1][q][2] =
1094 2.0 * polyx[i][0] * polyy[j][0] *
1098 const unsigned int dof_index2_1(
1099 cell_type2_offset1 + shift_ijk);
1100 const unsigned int dof_index2_2(
1101 cell_type2_offset2 + shift_ijk);
1103 data.shape_values[dof_index2_1][q][0] =
1104 data.shape_values[dof_index1][q][0];
1105 data.shape_values[dof_index2_1][q][1] =
1107 data.shape_values[dof_index1][q][1];
1108 data.shape_values[dof_index2_1][q][2] =
1109 data.shape_values[dof_index1][q][2];
1111 data.shape_values[dof_index2_2][q][0] =
1112 data.shape_values[dof_index1][q][0];
1113 data.shape_values[dof_index2_2][q][1] =
1115 data.shape_values[dof_index1][q][1];
1116 data.shape_values[dof_index2_2][q][2] =
1118 data.shape_values[dof_index1][q][2];
1122 const unsigned int shift_ij(
1125 const unsigned int dof_index3_1(
1126 cell_type3_offset1 + shift_ij);
1127 const unsigned int dof_index3_2(
1128 cell_type3_offset2 + shift_ij);
1129 const unsigned int dof_index3_3(
1130 cell_type3_offset3 + shift_ij);
1132 data.shape_values[dof_index3_1][q][0] =
1133 polyy[j][0] * polyz[k][0];
1134 data.shape_values[dof_index3_1][q][1] = 0.0;
1135 data.shape_values[dof_index3_1][q][2] = 0.0;
1137 data.shape_values[dof_index3_2][q][0] = 0.0;
1138 data.shape_values[dof_index3_2][q][1] =
1139 polyx[j][0] * polyz[k][0];
1140 data.shape_values[dof_index3_2][q][2] = 0.0;
1142 data.shape_values[dof_index3_3][q][0] = 0.0;
1143 data.shape_values[dof_index3_3][q][1] = 0.0;
1144 data.shape_values[dof_index3_3][q][2] =
1145 polyx[j][0] * polyy[k][0];
1151 for (
unsigned int k = 0; k < degree; ++k)
1153 const unsigned int shift_k(k * degree * degree);
1154 const unsigned int shift_j(
1157 for (
unsigned int j = 0; j < degree; ++j)
1159 const unsigned int shift_jk(j * degree +
1161 for (
unsigned int i = 0; i < degree; ++i)
1163 const unsigned int shift_ijk(shift_jk +
1167 const unsigned int dof_index1(
1168 cell_type1_offset + shift_ijk);
1170 data.shape_grads[dof_index1][q][0][0] =
1171 4.0 * polyx[i][2] * polyy[j][0] *
1173 data.shape_grads[dof_index1][q][0][1] =
1174 4.0 * polyx[i][1] * polyy[j][1] *
1176 data.shape_grads[dof_index1][q][0][2] =
1177 4.0 * polyx[i][1] * polyy[j][0] *
1180 data.shape_grads[dof_index1][q][1][0] =
1181 data.shape_grads[dof_index1][q][0][1];
1182 data.shape_grads[dof_index1][q][1][1] =
1183 4.0 * polyx[i][0] * polyy[j][2] *
1185 data.shape_grads[dof_index1][q][1][2] =
1186 4.0 * polyx[i][0] * polyy[j][1] *
1189 data.shape_grads[dof_index1][q][2][0] =
1190 data.shape_grads[dof_index1][q][0][2];
1191 data.shape_grads[dof_index1][q][2][1] =
1192 data.shape_grads[dof_index1][q][1][2];
1193 data.shape_grads[dof_index1][q][2][2] =
1194 4.0 * polyx[i][0] * polyy[j][0] *
1198 const unsigned int dof_index2_1(
1199 cell_type2_offset1 + shift_ijk);
1200 const unsigned int dof_index2_2(
1201 cell_type2_offset2 + shift_ijk);
1203 for (
unsigned int d = 0; d < dim; ++d)
1205 data.shape_grads[dof_index2_1][q][0]
1208 .shape_grads[dof_index1][q][0][d];
1209 data.shape_grads[dof_index2_1][q][1]
1213 .shape_grads[dof_index1][q][1][d];
1214 data.shape_grads[dof_index2_1][q][2]
1217 .shape_grads[dof_index1][q][2][d];
1219 data.shape_grads[dof_index2_2][q][0]
1222 .shape_grads[dof_index1][q][0][d];
1223 data.shape_grads[dof_index2_2][q][1]
1227 .shape_grads[dof_index1][q][1][d];
1228 data.shape_grads[dof_index2_2][q][2]
1232 .shape_grads[dof_index1][q][2][d];
1237 const unsigned int shift_ij(
1240 const unsigned int dof_index3_1(
1241 cell_type3_offset1 + shift_ij);
1242 const unsigned int dof_index3_2(
1243 cell_type3_offset2 + shift_ij);
1244 const unsigned int dof_index3_3(
1245 cell_type3_offset3 + shift_ij);
1246 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1248 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1250 data.shape_grads[dof_index3_1][q][d1]
1252 data.shape_grads[dof_index3_2][q][d1]
1254 data.shape_grads[dof_index3_3][q][d1]
1258 data.shape_grads[dof_index3_1][q][0][1] =
1259 2.0 * polyy[j][1] * polyz[k][0];
1260 data.shape_grads[dof_index3_1][q][0][2] =
1261 2.0 * polyy[j][0] * polyz[k][1];
1263 data.shape_grads[dof_index3_2][q][1][0] =
1264 2.0 * polyx[j][1] * polyz[k][0];
1265 data.shape_grads[dof_index3_2][q][1][2] =
1266 2.0 * polyx[j][0] * polyz[k][1];
1268 data.shape_grads[dof_index3_3][q][2][0] =
1269 2.0 * polyx[j][1] * polyy[k][0];
1270 data.shape_grads[dof_index3_3][q][2][1] =
1271 2.0 * polyx[j][0] * polyy[k][1];
1277 for (
unsigned int k = 0; k < degree; ++k)
1279 const unsigned int shift_k(k * degree * degree);
1280 const unsigned int shift_j(
1284 for (
unsigned int j = 0; j < degree; ++j)
1286 const unsigned int shift_jk(j * degree +
1288 for (
unsigned int i = 0; i < degree; ++i)
1290 const unsigned int shift_ijk(shift_jk +
1294 const unsigned int dof_index1(
1295 cell_type1_offset + shift_ijk);
1297 data.shape_hessians[dof_index1][q][0][0]
1299 8.0 * polyx[i][3] * polyy[j][0] *
1301 data.shape_hessians[dof_index1][q][1][0]
1303 8.0 * polyx[i][2] * polyy[j][1] *
1305 data.shape_hessians[dof_index1][q][2][0]
1307 8.0 * polyx[i][2] * polyy[j][0] *
1310 data.shape_hessians[dof_index1][q][0][1]
1312 data.shape_hessians[dof_index1][q][1][0]
1314 data.shape_hessians[dof_index1][q][1][1]
1316 8.0 * polyx[i][1] * polyy[j][2] *
1318 data.shape_hessians[dof_index1][q][2][1]
1320 8.0 * polyx[i][1] * polyy[j][1] *
1323 data.shape_hessians[dof_index1][q][0][2]
1325 data.shape_hessians[dof_index1][q][2][0]
1327 data.shape_hessians[dof_index1][q][1][2]
1329 data.shape_hessians[dof_index1][q][2][1]
1331 data.shape_hessians[dof_index1][q][2][2]
1333 8.0 * polyx[i][1] * polyy[j][0] *
1337 data.shape_hessians[dof_index1][q][0][0]
1339 data.shape_hessians[dof_index1][q][1][0]
1341 data.shape_hessians[dof_index1][q][1][0]
1343 data.shape_hessians[dof_index1][q][1][1]
1345 data.shape_hessians[dof_index1][q][2][0]
1347 data.shape_hessians[dof_index1][q][2][1]
1350 data.shape_hessians[dof_index1][q][0][1]
1352 data.shape_hessians[dof_index1][q][1][1]
1354 data.shape_hessians[dof_index1][q][1][1]
1356 8.0 * polyx[i][0] * polyy[j][3] *
1358 data.shape_hessians[dof_index1][q][2][1]
1360 8.0 * polyx[i][0] * polyy[j][2] *
1363 data.shape_hessians[dof_index1][q][0][2]
1365 data.shape_hessians[dof_index1][q][2][1]
1367 data.shape_hessians[dof_index1][q][1][2]
1369 data.shape_hessians[dof_index1][q][2][1]
1371 data.shape_hessians[dof_index1][q][2][2]
1373 8.0 * polyx[i][0] * polyy[j][1] *
1377 data.shape_hessians[dof_index1][q][0][0]
1379 data.shape_hessians[dof_index1][q][2][0]
1381 data.shape_hessians[dof_index1][q][1][0]
1383 data.shape_hessians[dof_index1][q][2][1]
1385 data.shape_hessians[dof_index1][q][2][0]
1387 data.shape_hessians[dof_index1][q][2][2]
1390 data.shape_hessians[dof_index1][q][0][1]
1392 data.shape_hessians[dof_index1][q][2][1]
1394 data.shape_hessians[dof_index1][q][1][1]
1396 data.shape_hessians[dof_index1][q][2][1]
1398 data.shape_hessians[dof_index1][q][2][1]
1400 data.shape_hessians[dof_index1][q][2][2]
1403 data.shape_hessians[dof_index1][q][0][2]
1405 data.shape_hessians[dof_index1][q][2][2]
1407 data.shape_hessians[dof_index1][q][1][2]
1409 data.shape_hessians[dof_index1][q][2][2]
1411 data.shape_hessians[dof_index1][q][2][2]
1413 8.0 * polyx[i][0] * polyy[j][0] *
1418 const unsigned int dof_index2_1(
1419 cell_type2_offset1 + shift_ijk);
1420 const unsigned int dof_index2_2(
1421 cell_type2_offset2 + shift_ijk);
1423 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1425 for (
unsigned int d2 = 0; d2 < dim;
1429 .shape_hessians[dof_index2_1][q]
1432 .shape_hessians[dof_index1][q]
1435 .shape_hessians[dof_index2_1][q]
1439 .shape_hessians[dof_index1][q]
1442 .shape_hessians[dof_index2_1][q]
1445 .shape_hessians[dof_index1][q]
1449 .shape_hessians[dof_index2_2][q]
1452 .shape_hessians[dof_index1][q]
1455 .shape_hessians[dof_index2_2][q]
1459 .shape_hessians[dof_index1][q]
1462 .shape_hessians[dof_index2_2][q]
1466 .shape_hessians[dof_index1][q]
1473 const unsigned int shift_ij(
1476 const unsigned int dof_index3_1(
1477 cell_type3_offset1 + shift_ij);
1478 const unsigned int dof_index3_2(
1479 cell_type3_offset2 + shift_ij);
1480 const unsigned int dof_index3_3(
1481 cell_type3_offset3 + shift_ij);
1482 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1484 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1486 for (
unsigned int d3 = 0; d3 < dim;
1490 .shape_hessians[dof_index3_1][q]
1494 .shape_hessians[dof_index3_2][q]
1498 .shape_hessians[dof_index3_3][q]
1505 .shape_hessians[dof_index3_1][q][0][1][1] =
1506 4.0 * polyy[j][2] * polyz[k][0];
1508 .shape_hessians[dof_index3_1][q][0][1][2] =
1509 4.0 * polyy[j][1] * polyz[k][1];
1512 .shape_hessians[dof_index3_1][q][0][2][1] =
1514 .shape_hessians[dof_index3_1][q][0][1][2];
1516 .shape_hessians[dof_index3_1][q][0][2][2] =
1517 4.0 * polyy[j][0] * polyz[k][2];
1521 .shape_hessians[dof_index3_2][q][1][0][0] =
1522 4.0 * polyx[j][2] * polyz[k][0];
1524 .shape_hessians[dof_index3_2][q][1][0][2] =
1525 4.0 * polyx[j][1] * polyz[k][1];
1528 .shape_hessians[dof_index3_2][q][1][2][0] =
1530 .shape_hessians[dof_index3_2][q][1][0][2];
1532 .shape_hessians[dof_index3_2][q][1][2][2] =
1533 4.0 * polyx[j][0] * polyz[k][2];
1537 .shape_hessians[dof_index3_3][q][2][0][0] =
1538 4.0 * polyx[j][2] * polyy[k][0];
1540 .shape_hessians[dof_index3_3][q][2][0][1] =
1541 4.0 * polyx[j][1] * polyy[k][1];
1544 .shape_hessians[dof_index3_3][q][2][1][0] =
1546 .shape_hessians[dof_index3_3][q][2][0][1];
1548 .shape_hessians[dof_index3_3][q][2][1][1] =
1549 4.0 * polyx[j][0] * polyy[k][2];
1589 const unsigned int n_q_points = quadrature.
size();
1592 fe_data.
shape_values.size() == this->n_dofs_per_cell(),
1594 this->n_dofs_per_cell()));
1600 const unsigned int degree(
1605 const unsigned int vertices_per_line(2);
1609 std::vector<std::vector<unsigned int>> edge_order(
1610 lines_per_cell, std::vector<unsigned int>(vertices_per_line));
1622 std::vector<int> edge_sign(lines_per_cell);
1623 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1625 unsigned int v0_loc =
1627 unsigned int v1_loc =
1629 unsigned int v0_glob = cell->vertex_index(v0_loc);
1630 unsigned int v1_glob = cell->vertex_index(v1_loc);
1635 if (cell->face(m)->at_boundary() ==
false)
1636 if (cell->neighbor_is_coarser(m))
1638 v0_glob = cell->parent()->vertex_index(v0_loc);
1639 v1_glob = cell->parent()->vertex_index(v1_loc);
1642 if (v0_glob > v1_glob)
1645 edge_sign[m] = -1.0;
1689 std::vector<std::vector<double>> edge_sigma_values(
1691 std::vector<std::vector<double>> edge_sigma_grads(
1694 std::vector<std::vector<double>> edge_lambda_values(
1696 std::vector<std::vector<double>> edge_lambda_grads(
1700 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1702 std::transform(edge_sigma_values[m].begin(),
1703 edge_sigma_values[m].end(),
1704 edge_sigma_values[m].begin(),
1705 [&](
const double edge_sigma_value) {
1706 return edge_sign[m] * edge_sigma_value;
1709 std::transform(edge_sigma_grads[m].begin(),
1710 edge_sigma_grads[m].end(),
1711 edge_sigma_grads[m].begin(),
1712 [&](
const double edge_sigma_grad) {
1713 return edge_sign[m] * edge_sigma_grad;
1720 const unsigned int poly_length =
1726 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1728 const unsigned int shift_m(m * this->n_dofs_per_line());
1729 for (
unsigned int q = 0; q < n_q_points; ++q)
1732 std::vector<std::vector<double>> poly(
1733 degree, std::vector<double>(poly_length));
1734 for (
unsigned int i = 1; i < degree + 1; ++i)
1739 IntegratedLegendrePolynomials[i + 1].value(
1740 edge_sigma_values[m][q], poly[i - 1]);
1745 for (
unsigned int d = 0; d < dim; ++d)
1748 0.5 * edge_sigma_grads[m][d] *
1749 edge_lambda_values[m][q];
1752 for (
unsigned int i = 1; i < degree + 1; ++i)
1754 const unsigned int poly_index = i - 1;
1755 const unsigned int dof_index(i + shift_m);
1756 for (
unsigned int d = 0; d < dim; ++d)
1759 edge_sigma_grads[m][d] *
1760 poly[poly_index][1] *
1761 edge_lambda_values[m][q] +
1762 poly[poly_index][0] *
1763 edge_lambda_grads[m][d];
1770 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1772 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1777 0.5 * edge_sigma_grads[m][d1] *
1778 edge_lambda_grads[m][d2];
1782 for (
unsigned int i = 1; i < degree + 1; ++i)
1784 const unsigned int poly_index = i - 1;
1785 const unsigned int dof_index(i + shift_m);
1788 edge_sigma_grads[m][0] *
1789 edge_sigma_grads[m][0] *
1790 edge_lambda_values[m][q] * poly[poly_index][2];
1793 (edge_sigma_grads[m][0] *
1794 edge_lambda_grads[m][1] +
1795 edge_sigma_grads[m][1] *
1796 edge_lambda_grads[m][0]) *
1797 poly[poly_index][1];
1803 edge_sigma_grads[m][1] *
1804 edge_sigma_grads[m][1] *
1805 edge_lambda_values[m][q] * poly[poly_index][2];
1811 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1813 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1815 for (
unsigned int d3 = 0; d3 < dim; ++d3)
1824 for (
unsigned int i = 0; i < degree; ++i)
1826 const unsigned int dof_index(i + 1 + shift_m);
1828 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1830 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1832 for (
unsigned int d3 = 0; d3 < dim; ++d3)
1836 edge_sigma_grads[m][d1] *
1837 edge_sigma_grads[m][d2] *
1838 edge_sigma_grads[m][d3] *
1840 edge_lambda_values[m][q] +
1842 (edge_sigma_grads[m][d1] *
1843 edge_sigma_grads[m][d2] *
1844 edge_lambda_grads[m][d3] +
1845 edge_sigma_grads[m][d3] *
1846 edge_sigma_grads[m][d1] *
1847 edge_lambda_grads[m][d2] +
1848 edge_sigma_grads[m][d3] *
1849 edge_sigma_grads[m][d2] *
1850 edge_lambda_grads[m][d1]);
1868 std::vector<int> edge_sign(lines_per_cell);
1869 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1871 unsigned int v0_loc =
1873 unsigned int v1_loc =
1875 unsigned int v0_glob = cell->vertex_index(v0_loc);
1876 unsigned int v1_glob = cell->vertex_index(v1_loc);
1878 if (v0_glob > v1_glob)
1881 edge_sign[m] = -1.0;
1892 for (
unsigned int f = 0; f < 6; ++f)
1893 if (cell->face(f)->at_boundary() ==
false)
1894 if (cell->neighbor_is_coarser(f))
1895 for (
unsigned int m = 0; m < 4; ++m)
1897 unsigned int parent_m =
1900 unsigned int v0_loc =
1902 unsigned int v1_loc =
1905 unsigned int v0_glob =
1906 cell->parent()->vertex_index(v0_loc);
1907 unsigned int v1_glob =
1908 cell->parent()->vertex_index(v1_loc);
1910 if (v0_glob > v1_glob)
1914 edge_sign[parent_m] = -1.0;
1920 edge_sign[parent_m] = 1.0;
1950 std::vector<unsigned int> adjacent_faces = {2, 2, 4, 4, 0, 0};
1951 for (
unsigned int f = 0; f < 6; ++f)
1953 if (!cell->face(f)->at_boundary() &&
1954 !cell->face(adjacent_faces[f])->at_boundary())
1955 if (!cell->neighbor_is_coarser(f) &&
1956 !cell->neighbor_is_coarser(adjacent_faces[f]))
1957 if (!cell->neighbor(f)
1958 ->face(adjacent_faces[f])
1960 if (cell->neighbor(f)->neighbor_is_coarser(
1963 unsigned int parent_m =
1965 unsigned int v0_loc =
1968 unsigned int v1_loc =
1971 unsigned int v0_glob =
1972 cell->parent()->vertex_index(v0_loc);
1973 unsigned int v1_glob =
1974 cell->parent()->vertex_index(v1_loc);
1975 if (v0_glob > v1_glob)
1979 edge_sign[parent_m] = -1.0;
1985 edge_sign[parent_m] = 1.0;
1989 if (!cell->face(f)->at_boundary() &&
1990 !cell->face(adjacent_faces[f] + 1)->at_boundary())
1991 if (!cell->neighbor_is_coarser(f) &&
1992 !cell->neighbor_is_coarser(adjacent_faces[f] + 1))
1993 if (!cell->neighbor(f)
1994 ->face(adjacent_faces[f] + 1)
1996 if (cell->neighbor(f)->neighbor_is_coarser(
1997 adjacent_faces[f] + 1))
1999 unsigned int parent_m =
2001 unsigned int v0_loc =
2004 unsigned int v1_loc =
2007 unsigned int v0_glob =
2008 cell->parent()->vertex_index(v0_loc);
2009 unsigned int v1_glob =
2010 cell->parent()->vertex_index(v1_loc);
2011 if (v0_glob > v1_glob)
2015 edge_sign[parent_m] = -1.0;
2021 edge_sign[parent_m] = 1.0;
2063 std::vector<std::vector<double>> edge_sigma_values(
2065 std::vector<std::vector<double>> edge_lambda_values(
2067 std::vector<std::vector<double>> edge_sigma_grads(
2069 std::vector<std::vector<std::vector<double>>> edge_lambda_grads(
2071 std::vector<std::vector<std::vector<double>>>
2075 for (
unsigned int m = 0; m < lines_per_cell; ++m)
2077 std::transform(edge_sigma_values[m].begin(),
2078 edge_sigma_values[m].end(),
2079 edge_sigma_values[m].begin(),
2080 [&](
const double edge_sigma_value) {
2081 return edge_sign[m] * edge_sigma_value;
2083 std::transform(edge_sigma_grads[m].begin(),
2084 edge_sigma_grads[m].end(),
2085 edge_sigma_grads[m].begin(),
2086 [&](
const double edge_sigma_grad) {
2087 return edge_sign[m] * edge_sigma_grad;
2095 const unsigned int poly_length =
2100 std::vector<std::vector<double>> poly(
2101 degree, std::vector<double>(poly_length));
2102 for (
unsigned int m = 0; m < lines_per_cell; ++m)
2104 const unsigned int shift_m(m * this->n_dofs_per_line());
2105 for (
unsigned int q = 0; q < n_q_points; ++q)
2111 for (
unsigned int i = 0; i < degree; ++i)
2113 IntegratedLegendrePolynomials[i + 2].value(
2114 edge_sigma_values[m][q], poly[i]);
2120 for (
unsigned int d = 0; d < dim; ++d)
2123 0.5 * edge_sigma_grads[m][d] *
2124 edge_lambda_values[m][q];
2129 for (
unsigned int i = 0; i < degree; ++i)
2131 const unsigned int dof_index(i + 1 + shift_m);
2132 for (
unsigned int d = 0; d < dim; ++d)
2135 edge_sigma_grads[m][d] * poly[i][1] *
2136 edge_lambda_values[m][q] +
2137 poly[i][0] * edge_lambda_grads[m][q][d];
2145 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2147 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2150 0.5 * edge_sigma_grads[m][d1] *
2151 edge_lambda_grads[m][q][d2];
2157 for (
unsigned int i = 0; i < degree; ++i)
2159 const unsigned int dof_index(i + 1 + shift_m);
2161 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2163 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2167 edge_sigma_grads[m][d1] *
2168 edge_sigma_grads[m][d2] *
2170 edge_lambda_values[m][q] +
2171 (edge_sigma_grads[m][d1] *
2172 edge_lambda_grads[m][q][d2] +
2173 edge_sigma_grads[m][d2] *
2174 edge_lambda_grads[m][q][d1]) *
2176 edge_lambda_gradgrads_3d[m][d1]
2187 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2189 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2191 for (
unsigned int d3 = 0; d3 < dim; ++d3)
2195 0.5 * edge_sigma_grads[m][d1] *
2196 edge_lambda_gradgrads_3d[m][d3][d2];
2204 for (
unsigned int i = 0; i < degree; ++i)
2206 const unsigned int dof_index(i + 1 + shift_m);
2208 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2210 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2212 for (
unsigned int d3 = 0; d3 < dim;
2218 edge_sigma_grads[m][d1] *
2219 edge_sigma_grads[m][d2] *
2220 edge_sigma_grads[m][d3] *
2222 edge_lambda_values[m][q] +
2224 (edge_sigma_grads[m][d1] *
2225 edge_sigma_grads[m][d2] *
2226 edge_lambda_grads[m][q]
2228 edge_sigma_grads[m][d3] *
2229 edge_sigma_grads[m][d1] *
2230 edge_lambda_grads[m][q]
2232 edge_sigma_grads[m][d3] *
2233 edge_sigma_grads[m][d2] *
2234 edge_lambda_grads[m][q]
2237 (edge_sigma_grads[m][d1] *
2238 edge_lambda_gradgrads_3d
2240 edge_sigma_grads[m][d2] *
2241 edge_lambda_gradgrads_3d
2243 edge_sigma_grads[m][d3] *
2244 edge_lambda_gradgrads_3d
2286 const unsigned int degree(
2297 const unsigned int n_q_points = quadrature.
size();
2300 fe_data.
shape_values.size() == this->n_dofs_per_cell(),
2302 this->n_dofs_per_cell()));
2309 const unsigned int vertices_per_face(
2314 const unsigned int n_line_dofs =
2318 std::vector<std::vector<unsigned int>> face_orientation(
2319 faces_per_cell, std::vector<unsigned int>(vertices_per_face));
2329 {1, 2}, {0, 3}, {3, 0}, {2, 1}};
2331 for (
unsigned int m = 0; m < faces_per_cell; ++m)
2334 bool cell_has_coarser_neighbor =
false;
2335 if (cell->face(m)->at_boundary() ==
false)
2336 if (cell->neighbor_is_coarser(m))
2337 cell_has_coarser_neighbor =
true;
2341 unsigned int current_max = 0;
2346 if (cell_has_coarser_neighbor)
2348 unsigned int current_glob = cell->parent()->vertex_index(
2350 for (
unsigned int v = 1; v < vertices_per_face; ++v)
2353 cell->parent()->vertex_index(
2357 current_glob = cell->parent()->vertex_index(
2366 unsigned int current_glob = cell->vertex_index(
2368 for (
unsigned int v = 1; v < vertices_per_face; ++v)
2375 current_glob = cell->vertex_index(
2381 face_orientation[m][0] =
2386 m, vertex_opposite_on_face[current_max]);
2394 if (cell_has_coarser_neighbor)
2396 if (cell->parent()->vertex_index(
2398 m, vertices_adjacent_on_face[current_max][0])) >
2399 cell->parent()->vertex_index(
2401 m, vertices_adjacent_on_face[current_max][1])))
2403 face_orientation[m][1] =
2405 m, vertices_adjacent_on_face[current_max][0]);
2406 face_orientation[m][3] =
2408 m, vertices_adjacent_on_face[current_max][1]);
2412 face_orientation[m][1] =
2414 m, vertices_adjacent_on_face[current_max][1]);
2415 face_orientation[m][3] =
2417 m, vertices_adjacent_on_face[current_max][0]);
2422 if (cell->vertex_index(
2424 m, vertices_adjacent_on_face[current_max][0])) >
2427 m, vertices_adjacent_on_face[current_max][1])))
2429 face_orientation[m][1] =
2431 m, vertices_adjacent_on_face[current_max][0]);
2432 face_orientation[m][3] =
2434 m, vertices_adjacent_on_face[current_max][1]);
2438 face_orientation[m][1] =
2440 m, vertices_adjacent_on_face[current_max][1]);
2441 face_orientation[m][3] =
2443 m, vertices_adjacent_on_face[current_max][0]);
2449 std::vector<std::vector<double>> face_xi_values(
2450 faces_per_cell, std::vector<double>(n_q_points));
2451 std::vector<std::vector<double>> face_xi_grads(
2452 faces_per_cell, std::vector<double>(dim));
2453 std::vector<std::vector<double>> face_eta_values(
2454 faces_per_cell, std::vector<double>(n_q_points));
2455 std::vector<std::vector<double>> face_eta_grads(
2456 faces_per_cell, std::vector<double>(dim));
2458 std::vector<std::vector<double>> face_lambda_values(
2459 faces_per_cell, std::vector<double>(n_q_points));
2460 std::vector<std::vector<double>> face_lambda_grads(
2461 faces_per_cell, std::vector<double>(dim));
2462 for (
unsigned int m = 0; m < faces_per_cell; ++m)
2464 for (
unsigned int q = 0; q < n_q_points; ++q)
2466 face_xi_values[m][q] =
2468 [face_orientation[m][1]];
2469 face_eta_values[m][q] =
2471 [face_orientation[m][3]];
2474 for (
unsigned int d = 0; d < dim; ++d)
2476 face_xi_grads[m][d] =
2478 [face_orientation[m][1]][d];
2479 face_eta_grads[m][d] =
2481 [face_orientation[m][3]][d];
2487 const unsigned int poly_length =
2492 std::vector<std::vector<double>> polyxi(
2493 degree, std::vector<double>(poly_length));
2494 std::vector<std::vector<double>> polyeta(
2495 degree, std::vector<double>(poly_length));
2498 for (
unsigned int m = 0; m < faces_per_cell; ++m)
2501 const unsigned int shift_m(m * this->n_dofs_per_quad(0));
2509 const unsigned int face_type1_offset(n_line_dofs + shift_m);
2519 const unsigned int face_type2_offset(face_type1_offset +
2532 const unsigned int face_type3_offset1(face_type2_offset +
2534 const unsigned int face_type3_offset2(face_type3_offset1 +
2538 for (
unsigned int q = 0; q < n_q_points; ++q)
2547 for (
unsigned int i = 0; i < degree; ++i)
2550 IntegratedLegendrePolynomials[i + 2].value(
2551 face_xi_values[m][q], polyxi[i]);
2552 IntegratedLegendrePolynomials[i + 2].value(
2553 face_eta_values[m][q], polyeta[i]);
2558 for (
unsigned int j = 0; j < degree; ++j)
2560 const unsigned int shift_j(j * degree);
2561 for (
unsigned int i = 0; i < degree; ++i)
2563 const unsigned int shift_ij(shift_j + i);
2565 const unsigned int dof_index1(face_type1_offset +
2567 for (
unsigned int d = 0; d < dim; ++d)
2570 (face_xi_grads[m][d] * polyxi[i][1] *
2572 face_eta_grads[m][d] * polyxi[i][0] *
2574 face_lambda_values[m][q] +
2575 face_lambda_grads[m][d] * polyxi[i][0] *
2579 const unsigned int dof_index2(face_type2_offset +
2581 for (
unsigned int d = 0; d < dim; ++d)
2584 (face_xi_grads[m][d] * polyxi[i][1] *
2586 face_eta_grads[m][d] * polyxi[i][0] *
2588 face_lambda_values[m][q];
2592 const unsigned int dof_index3_1(face_type3_offset1 +
2594 const unsigned int dof_index3_2(face_type3_offset2 +
2596 for (
unsigned int d = 0; d < dim; ++d)
2599 face_xi_grads[m][d] * polyeta[j][0] *
2600 face_lambda_values[m][q];
2603 face_eta_grads[m][d] * polyxi[j][0] *
2604 face_lambda_values[m][q];
2610 for (
unsigned int j = 0; j < degree; ++j)
2612 const unsigned int shift_j(j * degree);
2613 for (
unsigned int i = 0; i < degree; ++i)
2615 const unsigned int shift_ij(shift_j + i);
2617 const unsigned int dof_index1(face_type1_offset +
2619 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2621 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2625 (face_xi_grads[m][d1] *
2626 face_xi_grads[m][d2] * polyxi[i][2] *
2628 (face_xi_grads[m][d1] *
2629 face_eta_grads[m][d2] +
2630 face_xi_grads[m][d2] *
2631 face_eta_grads[m][d1]) *
2632 polyxi[i][1] * polyeta[j][1] +
2633 face_eta_grads[m][d1] *
2634 face_eta_grads[m][d2] *
2635 polyxi[i][0] * polyeta[j][2]) *
2636 face_lambda_values[m][q] +
2637 (face_xi_grads[m][d2] * polyxi[i][1] *
2639 face_eta_grads[m][d2] * polyxi[i][0] *
2641 face_lambda_grads[m][d1] +
2642 (face_xi_grads[m][d1] * polyxi[i][1] *
2644 face_eta_grads[m][d1] * polyxi[i][0] *
2646 face_lambda_grads[m][d2];
2650 const unsigned int dof_index2(face_type2_offset +
2652 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2654 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2658 (face_xi_grads[m][d1] *
2659 face_xi_grads[m][d2] * polyxi[i][2] *
2661 (face_xi_grads[m][d1] *
2662 face_eta_grads[m][d2] -
2663 face_xi_grads[m][d2] *
2664 face_eta_grads[m][d1]) *
2665 polyxi[i][1] * polyeta[j][1] -
2666 face_eta_grads[m][d1] *
2667 face_eta_grads[m][d2] *
2668 polyxi[i][0] * polyeta[j][2]) *
2669 face_lambda_values[m][q] +
2670 (face_xi_grads[m][d1] * polyxi[i][1] *
2672 face_eta_grads[m][d1] * polyxi[i][0] *
2674 face_lambda_grads[m][d2];
2679 const unsigned int dof_index3_1(face_type3_offset1 +
2681 const unsigned int dof_index3_2(face_type3_offset2 +
2683 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2685 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2688 face_xi_grads[m][d1] *
2689 (face_eta_grads[m][d2] * polyeta[j][1] *
2690 face_lambda_values[m][q] +
2691 face_lambda_grads[m][d2] * polyeta[j][0]);
2694 face_eta_grads[m][d1] *
2695 (face_xi_grads[m][d2] * polyxi[j][1] *
2696 face_lambda_values[m][q] +
2697 face_lambda_grads[m][d2] * polyxi[j][0]);
2704 for (
unsigned int j = 0; j < degree; ++j)
2706 const unsigned int shift_j(j * degree);
2707 for (
unsigned int i = 0; i < degree; ++i)
2709 const unsigned int shift_ij(shift_j + i);
2712 const unsigned int dof_index1(face_type1_offset +
2714 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2716 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2718 for (
unsigned int d3 = 0; d3 < dim; ++d3)
2723 face_xi_grads[m][d3] *
2724 (face_eta_grads[m][d1] *
2726 face_eta_grads[m][d2] *
2727 face_lambda_values[m][q] +
2729 face_lambda_grads[m][d2]) +
2731 face_eta_grads[m][d2] *
2732 face_lambda_grads[m][d1]) +
2735 face_eta_grads[m][d1] *
2736 face_eta_grads[m][d2] *
2737 face_eta_grads[m][d3] *
2738 face_lambda_values[m][q] +
2740 (face_eta_grads[m][d1] *
2741 face_eta_grads[m][d2] *
2742 face_lambda_grads[m][d3] +
2743 face_eta_grads[m][d3] *
2744 (face_eta_grads[m][d1] *
2745 face_lambda_grads[m]
2747 face_eta_grads[m][d2] *
2750 (polyxi[i][1] * polyeta[j][1] *
2751 face_eta_grads[m][d3] +
2752 polyxi[i][2] * polyeta[j][0] *
2753 face_xi_grads[m][d3]) *
2754 (face_xi_grads[m][d1] *
2755 face_lambda_grads[m][d2] +
2756 face_xi_grads[m][d2] *
2757 face_lambda_grads[m][d1]) +
2758 face_lambda_grads[m][d3] *
2759 (polyxi[i][2] * polyeta[j][0] *
2760 face_xi_grads[m][d1] *
2761 face_xi_grads[m][d2] +
2762 polyxi[i][1] * polyeta[j][1] *
2763 (face_xi_grads[m][d1] *
2764 face_eta_grads[m][d2] +
2765 face_xi_grads[m][d2] *
2766 face_eta_grads[m][d1])) +
2767 face_lambda_values[m][q] *
2768 (polyxi[i][3] * polyeta[j][0] *
2769 face_xi_grads[m][d1] *
2770 face_xi_grads[m][d2] *
2771 face_xi_grads[m][d3] +
2772 polyxi[i][1] * polyeta[j][2] *
2773 face_eta_grads[m][d3] *
2774 (face_xi_grads[m][d1] *
2775 face_eta_grads[m][d2] +
2776 face_xi_grads[m][d2] *
2777 face_eta_grads[m][d1]) +
2778 polyxi[i][2] * polyeta[j][1] *
2779 (face_xi_grads[m][d3] *
2780 face_xi_grads[m][d2] *
2781 face_eta_grads[m][d1] +
2782 face_xi_grads[m][d1] *
2783 (face_xi_grads[m][d2] *
2784 face_eta_grads[m][d3] +
2785 face_xi_grads[m][d3] *
2786 face_eta_grads[m][d2])));
2792 const unsigned int dof_index2(face_type2_offset +
2794 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2796 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2798 for (
unsigned int d3 = 0; d3 < dim; ++d3)
2802 face_xi_grads[m][d1] *
2803 (polyxi[i][1] * polyeta[j][1] *
2804 (face_eta_grads[m][d2] *
2805 face_lambda_grads[m][d3] +
2806 face_eta_grads[m][d3] *
2807 face_lambda_grads[m][d2]) +
2808 polyxi[i][2] * polyeta[j][0] *
2809 (face_xi_grads[m][d2] *
2810 face_lambda_grads[m][d3] +
2811 face_xi_grads[m][d3] *
2812 face_lambda_grads[m][d2]) +
2813 face_lambda_values[m][q] *
2814 (face_eta_grads[m][d2] *
2817 face_eta_grads[m][d3] +
2820 face_xi_grads[m][d3]) +
2821 face_xi_grads[m][d2] *
2824 face_eta_grads[m][d3] +
2827 face_xi_grads[m][d3]))) -
2829 face_eta_grads[m][d1] *
2830 (face_eta_grads[m][d2] *
2832 face_eta_grads[m][d3] *
2833 face_lambda_values[m][q] +
2835 face_lambda_grads[m][d3]) +
2837 face_eta_grads[m][d3] *
2838 face_lambda_grads[m][d2]) -
2839 face_eta_grads[m][d1] *
2841 face_xi_grads[m][d3] *
2843 face_eta_grads[m][d2] *
2844 face_lambda_values[m][q] +
2846 face_lambda_grads[m][d2]) +
2847 face_xi_grads[m][d2] *
2850 face_eta_grads[m][d3] *
2851 face_lambda_values[m]
2854 face_lambda_grads[m]
2856 polyxi[i][2] * polyeta[j][1] *
2857 face_xi_grads[m][d3] *
2858 face_lambda_values[m][q]));
2864 const unsigned int dof_index3_1(face_type3_offset1 +
2866 const unsigned int dof_index3_2(face_type3_offset2 +
2868 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2870 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2872 for (
unsigned int d3 = 0; d3 < dim; ++d3)
2876 face_xi_grads[m][d1] *
2877 (face_eta_grads[m][d2] *
2879 face_eta_grads[m][d3] *
2880 face_lambda_values[m][q] +
2882 face_lambda_grads[m][d3]) +
2883 face_lambda_grads[m][d2] *
2885 face_eta_grads[m][d3]);
2889 face_eta_grads[m][d1] *
2890 (face_xi_grads[m][d2] *
2892 face_xi_grads[m][d3] *
2893 face_lambda_values[m][q] +
2895 face_lambda_grads[m][d3]) +
2896 face_lambda_grads[m][d2] *
2897 polyxi[j][1] * face_xi_grads[m][d3]);