23 template <
int dim,
int spacedim>
30 std::vector<
bool>(compute_num_dofs(order), true),
32 std::vector<
bool>(dim, true)))
40 for (
unsigned int comp = 0; comp < this->
n_components(); ++comp)
52 template <
int dim,
int spacedim>
63 template <
int dim,
int spacedim>
68 const unsigned int )
const
77 template <
int dim,
int spacedim>
88 template <
int dim,
int spacedim>
93 const unsigned int )
const
101 template <
int dim,
int spacedim>
112 template <
int dim,
int spacedim>
117 const unsigned int )
const
125 template <
int dim,
int spacedim>
126 std::unique_ptr<typename ::FiniteElement<dim, spacedim>::InternalDataBase>
136 typename ::FiniteElement<dim, spacedim>::InternalDataBase>
137 data_ptr = std_cxx14::make_unique<InternalData>();
139 data.
update_each = requires_update_flags(update_flags);
142 const unsigned int degree(this->degree - 1);
148 const unsigned int n_line_dofs = this->dofs_per_line * lines_per_cell;
149 const unsigned int n_face_dofs = this->dofs_per_quad * faces_per_cell;
152 const unsigned int n_q_points = quadrature.
size();
155 data.sigma_imj_values.resize(
157 std::vector<std::vector<double>>(vertices_per_cell,
158 std::vector<double>(vertices_per_cell)));
160 data.sigma_imj_grads.resize(vertices_per_cell,
161 std::vector<std::vector<double>>(
162 vertices_per_cell, std::vector<double>(dim)));
167 data.shape_values.resize(this->dofs_per_cell,
173 data.shape_grads.resize(this->dofs_per_cell,
183 std::vector<Point<dim>> p_list(n_q_points);
193 std::vector<std::vector<double>> sigma(
194 n_q_points, std::vector<double>(lines_per_cell));
195 std::vector<std::vector<double>> lambda(
196 n_q_points, std::vector<double>(lines_per_cell));
198 for (
unsigned int q = 0; q < n_q_points; ++q)
200 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
201 sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
202 sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
203 sigma[q][3] = p_list[q][0] + p_list[q][1];
205 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
206 lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
207 lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
208 lambda[q][3] = p_list[q][0] * p_list[q][1];
209 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
211 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
213 data.sigma_imj_values[q][i][j] =
214 sigma[q][i] - sigma[q][j];
227 {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
228 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
229 unsigned int sigma_imj_component[vertices_per_cell]
232 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
234 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
241 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
242 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
246 sigma_imj_component[i][j] = 0;
247 for (
unsigned int d = 0;
d < dim; ++
d)
250 sigma_comp_signs[i][
d] - sigma_comp_signs[j][
d];
255 sigma_imj_component[i][j] =
d;
262 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
263 2.0 * sigma_imj_sign[i][j];
271 data.edge_sigma_values.resize(lines_per_cell);
272 data.edge_sigma_grads.resize(lines_per_cell);
273 for (
unsigned int m = 0; m < lines_per_cell; ++m)
275 data.edge_sigma_values[m].resize(n_q_points);
278 data.edge_sigma_grads[m].resize(dim);
288 data.edge_lambda_values.resize(lines_per_cell,
289 std::vector<double>(n_q_points));
290 data.edge_lambda_grads_2d.resize(lines_per_cell,
291 std::vector<double>(dim));
292 for (
unsigned int m = 0; m < lines_per_cell; ++m)
297 const unsigned int e1(
299 const unsigned int e2(
301 for (
unsigned int q = 0; q < n_q_points; ++q)
303 data.edge_sigma_values[m][q] =
304 data.sigma_imj_values[q][e2][e1];
305 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
308 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
310 data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
311 data.edge_lambda_grads_2d[1] = {1.0, 0.0};
312 data.edge_lambda_grads_2d[2] = {0.0, -1.0};
313 data.edge_lambda_grads_2d[3] = {0.0, 1.0};
381 const unsigned int cell_type1_offset = n_line_dofs;
388 const unsigned int cell_type2_offset =
389 cell_type1_offset + degree * degree;
397 const unsigned int cell_type3_offset1 =
398 cell_type2_offset + degree * degree;
399 const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
404 std::vector<Point<dim>> cell_points(n_q_points);
405 for (
unsigned int q = 0; q < n_q_points; ++q)
407 for (
unsigned int d = 0;
d < dim; ++
d)
409 cell_points[q][
d] = 2.0 * p_list[q][
d] - 1.0;
414 for (
unsigned int q = 0; q < n_q_points; ++q)
428 const unsigned int poly_length(
431 std::vector<std::vector<double>> polyx(
432 degree, std::vector<double>(poly_length));
433 std::vector<std::vector<double>> polyy(
434 degree, std::vector<double>(poly_length));
435 for (
unsigned int i = 0; i < degree; ++i)
440 IntegratedLegendrePolynomials[i + 2].value(
441 cell_points[q][0], polyx[i]);
442 IntegratedLegendrePolynomials[i + 2].value(
443 cell_points[q][1], polyy[i]);
448 for (
unsigned int j = 0; j < degree; ++j)
450 const unsigned int shift_j(j * degree);
451 for (
unsigned int i = 0; i < degree; ++i)
453 const unsigned int shift_ij(i + shift_j);
456 const unsigned int dof_index1(cell_type1_offset +
458 data.shape_values[dof_index1][q][0] =
459 2.0 * polyx[i][1] * polyy[j][0];
460 data.shape_values[dof_index1][q][1] =
461 2.0 * polyx[i][0] * polyy[j][1];
464 const unsigned int dof_index2(cell_type2_offset +
466 data.shape_values[dof_index2][q][0] =
467 data.shape_values[dof_index1][q][0];
468 data.shape_values[dof_index2][q][1] =
469 -1.0 * data.shape_values[dof_index1][q][1];
472 const unsigned int dof_index3_1(cell_type3_offset1 +
474 data.shape_values[dof_index3_1][q][0] = polyy[j][0];
475 data.shape_values[dof_index3_1][q][1] = 0.0;
477 const unsigned int dof_index3_2(cell_type3_offset2 +
479 data.shape_values[dof_index3_2][q][0] = 0.0;
480 data.shape_values[dof_index3_2][q][1] = polyx[j][0];
485 for (
unsigned int j = 0; j < degree; ++j)
487 const unsigned int shift_j(j * degree);
488 for (
unsigned int i = 0; i < degree; ++i)
490 const unsigned int shift_ij(i + shift_j);
493 const unsigned int dof_index1(cell_type1_offset +
495 data.shape_grads[dof_index1][q][0][0] =
496 4.0 * polyx[i][2] * polyy[j][0];
497 data.shape_grads[dof_index1][q][0][1] =
498 4.0 * polyx[i][1] * polyy[j][1];
499 data.shape_grads[dof_index1][q][1][0] =
500 data.shape_grads[dof_index1][q][0][1];
501 data.shape_grads[dof_index1][q][1][1] =
502 4.0 * polyx[i][0] * polyy[j][2];
505 const unsigned int dof_index2(cell_type2_offset +
507 data.shape_grads[dof_index2][q][0][0] =
508 data.shape_grads[dof_index1][q][0][0];
509 data.shape_grads[dof_index2][q][0][1] =
510 data.shape_grads[dof_index1][q][0][1];
511 data.shape_grads[dof_index2][q][1][0] =
512 -1.0 * data.shape_grads[dof_index1][q][1][0];
513 data.shape_grads[dof_index2][q][1][1] =
514 -1.0 * data.shape_grads[dof_index1][q][1][1];
517 const unsigned int dof_index3_1(cell_type3_offset1 +
519 data.shape_grads[dof_index3_1][q][0][0] = 0.0;
520 data.shape_grads[dof_index3_1][q][0][1] =
522 data.shape_grads[dof_index3_1][q][1][0] = 0.0;
523 data.shape_grads[dof_index3_1][q][1][1] = 0.0;
525 const unsigned int dof_index3_2(cell_type3_offset2 +
527 data.shape_grads[dof_index3_2][q][0][0] = 0.0;
528 data.shape_grads[dof_index3_2][q][0][1] = 0.0;
529 data.shape_grads[dof_index3_2][q][1][0] =
531 data.shape_grads[dof_index3_2][q][1][1] = 0.0;
542 std::vector<std::vector<double>> sigma(
543 n_q_points, std::vector<double>(lines_per_cell));
544 std::vector<std::vector<double>> lambda(
545 n_q_points, std::vector<double>(lines_per_cell));
546 for (
unsigned int q = 0; q < n_q_points; ++q)
548 sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
551 p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
553 (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
554 sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
556 (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
557 sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
558 sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
559 sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
561 lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
562 (1.0 - p_list[q][2]);
564 p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
566 (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
567 lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
569 (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
570 lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
571 lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
572 lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
576 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
578 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
580 data.sigma_imj_values[q][i][j] =
581 sigma[q][i] - sigma[q][j];
614 int sigma_imj_sign[vertices_per_cell][vertices_per_cell];
615 unsigned int sigma_imj_component[vertices_per_cell]
618 for (
unsigned int i = 0; i < vertices_per_cell; ++i)
620 for (
unsigned int j = 0; j < vertices_per_cell; ++j)
627 sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
628 sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
632 sigma_imj_component[i][j] = 0;
633 for (
unsigned int d = 0;
d < dim; ++
d)
636 sigma_comp_signs[i][
d] - sigma_comp_signs[j][
d];
641 sigma_imj_component[i][j] =
d;
648 data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
649 2.0 * sigma_imj_sign[i][j];
656 data.edge_sigma_values.resize(lines_per_cell);
657 data.edge_lambda_values.resize(lines_per_cell);
658 data.edge_sigma_grads.resize(lines_per_cell);
659 data.edge_lambda_grads_3d.resize(lines_per_cell);
660 data.edge_lambda_gradgrads_3d.resize(lines_per_cell);
661 for (
unsigned int m = 0; m < lines_per_cell; ++m)
663 data.edge_sigma_values[m].resize(n_q_points);
664 data.edge_lambda_values[m].resize(n_q_points);
667 data.edge_sigma_grads[m].resize(dim);
669 data.edge_lambda_grads_3d[m].resize(n_q_points);
670 for (
unsigned int q = 0; q < n_q_points; ++q)
672 data.edge_lambda_grads_3d[m][q].resize(dim);
676 data.edge_lambda_gradgrads_3d[m].resize(dim);
677 for (
unsigned int d = 0;
d < dim; ++
d)
679 data.edge_lambda_gradgrads_3d[m][
d].resize(dim);
686 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
688 for (
unsigned int m = 0; m < lines_per_cell; ++m)
693 const unsigned int e1(
695 const unsigned int e2(
698 for (
unsigned int q = 0; q < n_q_points; ++q)
700 data.edge_sigma_values[m][q] =
701 data.sigma_imj_values[q][e2][e1];
702 data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
705 data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
708 for (
unsigned int q = 0; q < n_q_points; ++q)
710 double x(p_list[q][0]);
711 double y(p_list[q][1]);
712 double z(p_list[q][2]);
713 data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0};
714 data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x};
715 data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0};
716 data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y};
717 data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x};
718 data.edge_lambda_grads_3d[5][q] = {z, 0.0, x};
719 data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y};
720 data.edge_lambda_grads_3d[7][q] = {0.0, z, y};
721 data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0};
722 data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0};
723 data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
724 data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
754 for (
unsigned int m = 0; m < lines_per_cell; ++m)
756 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
757 [edge_lambda_directions[m][1]] =
759 data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
760 [edge_lambda_directions[m][0]] =
768 data.face_lambda_values.resize(faces_per_cell);
769 data.face_lambda_grads.resize(faces_per_cell);
771 for (
unsigned int m = 0; m < faces_per_cell; ++m)
773 data.face_lambda_values[m].resize(n_q_points);
774 data.face_lambda_grads[m].resize(3);
777 for (
unsigned int q = 0; q < n_q_points; ++q)
779 double x(p_list[q][0]);
780 double y(p_list[q][1]);
781 double z(p_list[q][2]);
782 data.face_lambda_values[0][q] = 1.0 - x;
783 data.face_lambda_values[1][q] = x;
784 data.face_lambda_values[2][q] = 1.0 - y;
785 data.face_lambda_values[3][q] = y;
786 data.face_lambda_values[4][q] = 1.0 - z;
787 data.face_lambda_values[5][q] = z;
790 data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
791 data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
792 data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
793 data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
794 data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
795 data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
808 const unsigned int cell_type1_offset(n_line_dofs +
821 const unsigned int cell_type2_offset1(
822 cell_type1_offset + degree * degree * degree);
823 const unsigned int cell_type2_offset2(
824 cell_type2_offset1 + degree * degree * degree);
834 const unsigned int cell_type3_offset1(
835 cell_type2_offset2 + degree * degree * degree);
836 const unsigned int cell_type3_offset2(cell_type3_offset1 +
838 const unsigned int cell_type3_offset3(cell_type3_offset2 +
842 std::vector<Point<dim>> cell_points(n_q_points);
843 for (
unsigned int q = 0; q < n_q_points; ++q)
845 for (
unsigned int d = 0;
d < dim; ++
d)
847 cell_points[q][
d] = 2.0 * p_list[q][
d] - 1.0;
854 const unsigned int poly_length(
857 for (
unsigned int q = 0; q < n_q_points; ++q)
865 std::vector<std::vector<double>> polyx(
866 degree, std::vector<double>(poly_length));
867 std::vector<std::vector<double>> polyy(
868 degree, std::vector<double>(poly_length));
869 std::vector<std::vector<double>> polyz(
870 degree, std::vector<double>(poly_length));
871 for (
unsigned int i = 0; i < degree; ++i)
874 IntegratedLegendrePolynomials[i + 2].value(
875 cell_points[q][0], polyx[i]);
876 IntegratedLegendrePolynomials[i + 2].value(
877 cell_points[q][1], polyy[i]);
878 IntegratedLegendrePolynomials[i + 2].value(
879 cell_points[q][2], polyz[i]);
884 for (
unsigned int k = 0; k < degree; ++k)
886 const unsigned int shift_k(k * degree * degree);
887 const unsigned int shift_j(
890 for (
unsigned int j = 0; j < degree; ++j)
892 const unsigned int shift_jk(j * degree +
894 for (
unsigned int i = 0; i < degree; ++i)
896 const unsigned int shift_ijk(shift_jk +
900 const unsigned int dof_index1(
901 cell_type1_offset + shift_ijk);
903 data.shape_values[dof_index1][q][0] =
904 2.0 * polyx[i][1] * polyy[j][0] *
906 data.shape_values[dof_index1][q][1] =
907 2.0 * polyx[i][0] * polyy[j][1] *
909 data.shape_values[dof_index1][q][2] =
910 2.0 * polyx[i][0] * polyy[j][0] *
914 const unsigned int dof_index2_1(
915 cell_type2_offset1 + shift_ijk);
916 const unsigned int dof_index2_2(
917 cell_type2_offset2 + shift_ijk);
919 data.shape_values[dof_index2_1][q][0] =
920 data.shape_values[dof_index1][q][0];
921 data.shape_values[dof_index2_1][q][1] =
923 data.shape_values[dof_index1][q][1];
924 data.shape_values[dof_index2_1][q][2] =
925 data.shape_values[dof_index1][q][2];
927 data.shape_values[dof_index2_2][q][0] =
928 data.shape_values[dof_index1][q][0];
929 data.shape_values[dof_index2_2][q][1] =
931 data.shape_values[dof_index1][q][1];
932 data.shape_values[dof_index2_2][q][2] =
934 data.shape_values[dof_index1][q][2];
938 const unsigned int shift_ij(
941 const unsigned int dof_index3_1(
942 cell_type3_offset1 + shift_ij);
943 const unsigned int dof_index3_2(
944 cell_type3_offset2 + shift_ij);
945 const unsigned int dof_index3_3(
946 cell_type3_offset3 + shift_ij);
948 data.shape_values[dof_index3_1][q][0] =
949 polyy[j][0] * polyz[k][0];
950 data.shape_values[dof_index3_1][q][1] = 0.0;
951 data.shape_values[dof_index3_1][q][2] = 0.0;
953 data.shape_values[dof_index3_2][q][0] = 0.0;
954 data.shape_values[dof_index3_2][q][1] =
955 polyx[j][0] * polyz[k][0];
956 data.shape_values[dof_index3_2][q][2] = 0.0;
958 data.shape_values[dof_index3_3][q][0] = 0.0;
959 data.shape_values[dof_index3_3][q][1] = 0.0;
960 data.shape_values[dof_index3_3][q][2] =
961 polyx[j][0] * polyy[k][0];
967 for (
unsigned int k = 0; k < degree; ++k)
969 const unsigned int shift_k(k * degree * degree);
970 const unsigned int shift_j(
973 for (
unsigned int j = 0; j < degree; ++j)
975 const unsigned int shift_jk(j * degree +
977 for (
unsigned int i = 0; i < degree; ++i)
979 const unsigned int shift_ijk(shift_jk +
983 const unsigned int dof_index1(
984 cell_type1_offset + shift_ijk);
986 data.shape_grads[dof_index1][q][0][0] =
987 4.0 * polyx[i][2] * polyy[j][0] *
989 data.shape_grads[dof_index1][q][0][1] =
990 4.0 * polyx[i][1] * polyy[j][1] *
992 data.shape_grads[dof_index1][q][0][2] =
993 4.0 * polyx[i][1] * polyy[j][0] *
996 data.shape_grads[dof_index1][q][1][0] =
997 data.shape_grads[dof_index1][q][0][1];
998 data.shape_grads[dof_index1][q][1][1] =
999 4.0 * polyx[i][0] * polyy[j][2] *
1001 data.shape_grads[dof_index1][q][1][2] =
1002 4.0 * polyx[i][0] * polyy[j][1] *
1005 data.shape_grads[dof_index1][q][2][0] =
1006 data.shape_grads[dof_index1][q][0][2];
1007 data.shape_grads[dof_index1][q][2][1] =
1008 data.shape_grads[dof_index1][q][1][2];
1009 data.shape_grads[dof_index1][q][2][2] =
1010 4.0 * polyx[i][0] * polyy[j][0] *
1014 const unsigned int dof_index2_1(
1015 cell_type2_offset1 + shift_ijk);
1016 const unsigned int dof_index2_2(
1017 cell_type2_offset2 + shift_ijk);
1019 for (
unsigned int d = 0;
d < dim; ++
d)
1021 data.shape_grads[dof_index2_1][q][0]
1024 .shape_grads[dof_index1][q][0][
d];
1025 data.shape_grads[dof_index2_1][q][1]
1029 .shape_grads[dof_index1][q][1][
d];
1030 data.shape_grads[dof_index2_1][q][2]
1033 .shape_grads[dof_index1][q][2][
d];
1035 data.shape_grads[dof_index2_2][q][0]
1038 .shape_grads[dof_index1][q][0][
d];
1039 data.shape_grads[dof_index2_2][q][1]
1043 .shape_grads[dof_index1][q][1][
d];
1044 data.shape_grads[dof_index2_2][q][2]
1048 .shape_grads[dof_index1][q][2][
d];
1053 const unsigned int shift_ij(
1056 const unsigned int dof_index3_1(
1057 cell_type3_offset1 + shift_ij);
1058 const unsigned int dof_index3_2(
1059 cell_type3_offset2 + shift_ij);
1060 const unsigned int dof_index3_3(
1061 cell_type3_offset3 + shift_ij);
1062 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1064 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1066 data.shape_grads[dof_index3_1][q][d1]
1068 data.shape_grads[dof_index3_2][q][d1]
1070 data.shape_grads[dof_index3_3][q][d1]
1074 data.shape_grads[dof_index3_1][q][0][1] =
1075 2.0 * polyy[j][1] * polyz[k][0];
1076 data.shape_grads[dof_index3_1][q][0][2] =
1077 2.0 * polyy[j][0] * polyz[k][1];
1079 data.shape_grads[dof_index3_2][q][1][0] =
1080 2.0 * polyx[j][1] * polyz[k][0];
1081 data.shape_grads[dof_index3_2][q][1][2] =
1082 2.0 * polyx[j][0] * polyz[k][1];
1084 data.shape_grads[dof_index3_3][q][2][0] =
1085 2.0 * polyx[j][1] * polyy[k][0];
1086 data.shape_grads[dof_index3_3][q][2][1] =
1087 2.0 * polyx[j][0] * polyy[k][1];
1106 template <
int dim,
int spacedim>
1127 const unsigned int n_q_points = quadrature.
size();
1132 this->dofs_per_cell));
1138 const unsigned int degree(
1143 const unsigned int vertices_per_line(2);
1147 std::vector<std::vector<unsigned int>> edge_order(
1148 lines_per_cell, std::vector<unsigned int>(vertices_per_line));
1160 std::vector<int> edge_sign(lines_per_cell);
1161 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1163 unsigned int v0_loc =
1165 unsigned int v1_loc =
1167 unsigned int v0_glob = cell->vertex_index(v0_loc);
1168 unsigned int v1_glob = cell->vertex_index(v1_loc);
1170 if (v0_glob > v1_glob)
1173 edge_sign[m] = -1.0;
1217 std::vector<std::vector<double>> edge_sigma_values(
1219 std::vector<std::vector<double>> edge_sigma_grads(
1222 std::vector<std::vector<double>> edge_lambda_values(
1224 std::vector<std::vector<double>> edge_lambda_grads(
1228 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1231 edge_sigma_values[m].
end(),
1232 edge_sigma_values[m].
begin(),
1233 [&](
const double edge_sigma_value) {
1234 return edge_sign[m] * edge_sigma_value;
1238 edge_sigma_grads[m].
end(),
1239 edge_sigma_grads[m].
begin(),
1240 [&](
const double edge_sigma_grad) {
1241 return edge_sign[m] * edge_sigma_grad;
1251 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1253 const unsigned int shift_m(m * this->dofs_per_line);
1254 for (
unsigned int q = 0; q < n_q_points; ++q)
1257 std::vector<std::vector<double>> poly(
1258 degree, std::vector<double>(poly_length));
1259 for (
unsigned int i = 1; i < degree + 1; ++i)
1264 IntegratedLegendrePolynomials[i + 1].value(
1265 edge_sigma_values[m][q], poly[i - 1]);
1270 for (
unsigned int d = 0;
d < dim; ++
d)
1273 0.5 * edge_sigma_grads[m][
d] *
1274 edge_lambda_values[m][q];
1277 for (
unsigned int i = 1; i < degree + 1; ++i)
1279 const unsigned int poly_index = i - 1;
1280 const unsigned int dof_index(i + shift_m);
1281 for (
unsigned int d = 0;
d < dim; ++
d)
1284 edge_sigma_grads[m][
d] *
1285 poly[poly_index][1] *
1286 edge_lambda_values[m][q] +
1287 poly[poly_index][0] *
1288 edge_lambda_grads[m][
d];
1295 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1297 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1302 0.5 * edge_sigma_grads[m][d1] *
1303 edge_lambda_grads[m][d2];
1307 for (
unsigned int i = 1; i < degree + 1; ++i)
1309 const unsigned int poly_index = i - 1;
1310 const unsigned int dof_index(i + shift_m);
1313 edge_sigma_grads[m][0] *
1314 edge_sigma_grads[m][0] *
1315 edge_lambda_values[m][q] * poly[poly_index][2];
1318 (edge_sigma_grads[m][0] *
1319 edge_lambda_grads[m][1] +
1320 edge_sigma_grads[m][1] *
1321 edge_lambda_grads[m][0]) *
1322 poly[poly_index][1];
1328 edge_sigma_grads[m][1] *
1329 edge_sigma_grads[m][1] *
1330 edge_lambda_values[m][q] * poly[poly_index][2];
1345 std::vector<int> edge_sign(lines_per_cell);
1346 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1348 unsigned int v0_loc =
1350 unsigned int v1_loc =
1352 unsigned int v0_glob = cell->vertex_index(v0_loc);
1353 unsigned int v1_glob = cell->vertex_index(v1_loc);
1355 if (v0_glob > v1_glob)
1358 edge_sign[m] = -1.0;
1404 std::vector<std::vector<double>> edge_sigma_values(
1406 std::vector<std::vector<double>> edge_lambda_values(
1408 std::vector<std::vector<double>> edge_sigma_grads(
1410 std::vector<std::vector<std::vector<double>>> edge_lambda_grads(
1412 std::vector<std::vector<std::vector<double>>>
1416 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1419 edge_sigma_values[m].
end(),
1420 edge_sigma_values[m].
begin(),
1421 [&](
const double edge_sigma_value) {
1422 return edge_sign[m] * edge_sigma_value;
1425 edge_sigma_grads[m].
end(),
1426 edge_sigma_grads[m].
begin(),
1427 [&](
const double edge_sigma_grad) {
1428 return edge_sign[m] * edge_sigma_grad;
1438 std::vector<std::vector<double>> poly(
1439 degree, std::vector<double>(poly_length));
1440 for (
unsigned int m = 0; m < lines_per_cell; ++m)
1442 const unsigned int shift_m(m * this->dofs_per_line);
1443 for (
unsigned int q = 0; q < n_q_points; ++q)
1448 for (
unsigned int i = 0; i < degree; ++i)
1450 IntegratedLegendrePolynomials[i + 2].value(
1451 edge_sigma_values[m][q], poly[i]);
1457 for (
unsigned int d = 0;
d < dim; ++
d)
1460 0.5 * edge_sigma_grads[m][
d] *
1461 edge_lambda_values[m][q];
1465 for (
unsigned int i = 0; i < degree; ++i)
1467 const unsigned int dof_index(i + 1 + shift_m);
1468 for (
unsigned int d = 0;
d < dim; ++
d)
1471 edge_sigma_grads[m][
d] * poly[i][1] *
1472 edge_lambda_values[m][q] +
1473 poly[i][0] * edge_lambda_grads[m][q][
d];
1481 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1483 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1486 0.5 * edge_sigma_grads[m][d1] *
1487 edge_lambda_grads[m][q][d2];
1492 for (
unsigned int i = 0; i < degree; ++i)
1494 const unsigned int dof_index(i + 1 + shift_m);
1496 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1498 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1502 edge_sigma_grads[m][d1] *
1503 edge_sigma_grads[m][d2] *
1505 edge_lambda_values[m][q] +
1506 (edge_sigma_grads[m][d1] *
1507 edge_lambda_grads[m][q][d2] +
1508 edge_sigma_grads[m][d2] *
1509 edge_lambda_grads[m][q][d1]) *
1511 edge_lambda_gradgrads_3d[m][d1]
1533 template <
int dim,
int spacedim>
1553 const unsigned int degree(
1564 const unsigned int n_q_points = quadrature.
size();
1569 this->dofs_per_cell));
1576 const unsigned int vertices_per_face(
1581 const unsigned int n_line_dofs =
1585 std::vector<std::vector<unsigned int>> face_orientation(
1586 faces_per_cell, std::vector<unsigned int>(vertices_per_face));
1596 {1, 2}, {0, 3}, {0, 3}, {1, 2}};
1598 for (
unsigned int m = 0; m < faces_per_cell; ++m)
1602 unsigned int current_max = 0;
1603 unsigned int current_glob = cell->vertex_index(
1605 for (
unsigned int v = 1; v < vertices_per_face; ++v)
1612 current_glob = cell->vertex_index(
1616 face_orientation[m][0] =
1621 m, vertex_opposite_on_face[current_max]);
1626 m, vertices_adjacent_on_face[current_max][0])) >
1628 m, vertices_adjacent_on_face[current_max][1])))
1630 face_orientation[m][1] =
1632 m, vertices_adjacent_on_face[current_max][0]);
1633 face_orientation[m][3] =
1635 m, vertices_adjacent_on_face[current_max][1]);
1639 face_orientation[m][1] =
1641 m, vertices_adjacent_on_face[current_max][1]);
1642 face_orientation[m][3] =
1644 m, vertices_adjacent_on_face[current_max][0]);
1650 std::vector<std::vector<double>> face_xi_values(
1651 faces_per_cell, std::vector<double>(n_q_points));
1652 std::vector<std::vector<double>> face_xi_grads(
1653 faces_per_cell, std::vector<double>(dim));
1654 std::vector<std::vector<double>> face_eta_values(
1655 faces_per_cell, std::vector<double>(n_q_points));
1656 std::vector<std::vector<double>> face_eta_grads(
1657 faces_per_cell, std::vector<double>(dim));
1659 std::vector<std::vector<double>> face_lambda_values(
1660 faces_per_cell, std::vector<double>(n_q_points));
1661 std::vector<std::vector<double>> face_lambda_grads(
1662 faces_per_cell, std::vector<double>(dim));
1663 for (
unsigned int m = 0; m < faces_per_cell; ++m)
1665 for (
unsigned int q = 0; q < n_q_points; ++q)
1667 face_xi_values[m][q] =
1669 [face_orientation[m][1]];
1670 face_eta_values[m][q] =
1672 [face_orientation[m][3]];
1675 for (
unsigned int d = 0;
d < dim; ++
d)
1677 face_xi_grads[m][
d] =
1679 [face_orientation[m][1]][
d];
1680 face_eta_grads[m][
d] =
1682 [face_orientation[m][3]][
d];
1689 std::vector<std::vector<double>> polyxi(
1690 degree, std::vector<double>(poly_length));
1691 std::vector<std::vector<double>> polyeta(
1692 degree, std::vector<double>(poly_length));
1695 for (
unsigned int m = 0; m < faces_per_cell; ++m)
1697 const unsigned int shift_m(m * this->dofs_per_quad);
1705 const unsigned int face_type1_offset(n_line_dofs + shift_m);
1715 const unsigned int face_type2_offset(face_type1_offset +
1727 const unsigned int face_type3_offset1(face_type2_offset +
1729 const unsigned int face_type3_offset2(face_type3_offset1 +
1733 for (
unsigned int q = 0; q < n_q_points; ++q)
1742 for (
unsigned int i = 0; i < degree; ++i)
1745 IntegratedLegendrePolynomials[i + 2].value(
1746 face_xi_values[m][q], polyxi[i]);
1747 IntegratedLegendrePolynomials[i + 2].value(
1748 face_eta_values[m][q], polyeta[i]);
1753 for (
unsigned int j = 0; j < degree; ++j)
1755 const unsigned int shift_j(j * degree);
1756 for (
unsigned int i = 0; i < degree; ++i)
1758 const unsigned int shift_ij(shift_j + i);
1760 const unsigned int dof_index1(face_type1_offset +
1762 for (
unsigned int d = 0;
d < dim; ++
d)
1765 (face_xi_grads[m][
d] * polyxi[i][1] *
1767 face_eta_grads[m][
d] * polyxi[i][0] *
1769 face_lambda_values[m][q] +
1770 face_lambda_grads[m][
d] * polyxi[i][0] *
1774 const unsigned int dof_index2(face_type2_offset +
1776 for (
unsigned int d = 0;
d < dim; ++
d)
1779 (face_xi_grads[m][
d] * polyxi[i][1] *
1781 face_eta_grads[m][
d] * polyxi[i][0] *
1783 face_lambda_values[m][q];
1787 const unsigned int dof_index3_1(face_type3_offset1 +
1789 const unsigned int dof_index3_2(face_type3_offset2 +
1791 for (
unsigned int d = 0;
d < dim; ++
d)
1794 face_xi_grads[m][
d] * polyeta[j][0] *
1795 face_lambda_values[m][q];
1798 face_eta_grads[m][
d] * polyxi[j][0] *
1799 face_lambda_values[m][q];
1805 for (
unsigned int j = 0; j < degree; ++j)
1807 const unsigned int shift_j(j * degree);
1808 for (
unsigned int i = 0; i < degree; ++i)
1810 const unsigned int shift_ij(shift_j + i);
1812 const unsigned int dof_index1(face_type1_offset +
1814 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1816 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1820 (face_xi_grads[m][d1] *
1821 face_xi_grads[m][d2] * polyxi[i][2] *
1823 (face_xi_grads[m][d1] *
1824 face_eta_grads[m][d2] +
1825 face_xi_grads[m][d2] *
1826 face_eta_grads[m][d1]) *
1827 polyxi[i][1] * polyeta[j][1] +
1828 face_eta_grads[m][d1] *
1829 face_eta_grads[m][d2] *
1830 polyxi[i][0] * polyeta[j][2]) *
1831 face_lambda_values[m][q] +
1832 (face_xi_grads[m][d2] * polyxi[i][1] *
1834 face_eta_grads[m][d2] * polyxi[i][0] *
1836 face_lambda_grads[m][d1] +
1837 (face_xi_grads[m][d1] * polyxi[i][1] *
1839 face_eta_grads[m][d1] * polyxi[i][0] *
1841 face_lambda_grads[m][d2];
1845 const unsigned int dof_index2(face_type2_offset +
1847 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1849 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1853 (face_xi_grads[m][d1] *
1854 face_xi_grads[m][d2] * polyxi[i][2] *
1856 (face_xi_grads[m][d1] *
1857 face_eta_grads[m][d2] -
1858 face_xi_grads[m][d2] *
1859 face_eta_grads[m][d1]) *
1860 polyxi[i][1] * polyeta[j][1] -
1861 face_eta_grads[m][d1] *
1862 face_eta_grads[m][d2] *
1863 polyxi[i][0] * polyeta[j][2]) *
1864 face_lambda_values[m][q] +
1865 (face_xi_grads[m][d1] * polyxi[i][1] *
1867 face_eta_grads[m][d1] * polyxi[i][0] *
1869 face_lambda_grads[m][d2];
1874 const unsigned int dof_index3_1(face_type3_offset1 +
1876 const unsigned int dof_index3_2(face_type3_offset2 +
1878 for (
unsigned int d1 = 0; d1 < dim; ++d1)
1880 for (
unsigned int d2 = 0; d2 < dim; ++d2)
1883 face_xi_grads[m][d1] *
1884 (face_eta_grads[m][d2] * polyeta[j][1] *
1885 face_lambda_values[m][q] +
1886 face_lambda_grads[m][d2] * polyeta[j][0]);
1889 face_eta_grads[m][d1] *
1890 (face_xi_grads[m][d2] * polyxi[j][1] *
1891 face_lambda_values[m][q] +
1892 face_lambda_grads[m][d2] * polyxi[j][0]);
1909 template <
int dim,
int spacedim>
1917 const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
1931 fill_edge_values(cell, quadrature, fe_data);
1932 if (dim == 3 && this->degree > 1)
1934 fill_face_values(cell, quadrature, fe_data);
1938 const unsigned int n_q_points = quadrature.
size();
1943 this->dofs_per_cell));
1952 std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
1953 for (
unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
1955 const unsigned int first =
1957 this->get_nonzero_components(dof)
1958 .first_selected_component()];
1964 for (
unsigned int q = 0; q < n_q_points; ++q)
1966 for (
unsigned int d = 0;
d < dim; ++
d)
1969 transformed_shape_values[q][
d];
1978 std::vector<Tensor<2, dim>> input(n_q_points);
1979 std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
1980 for (
unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
1982 for (
unsigned int q = 0; q < n_q_points; ++q)
1991 const unsigned int first =
1993 this->get_nonzero_components(dof)
1994 .first_selected_component()];
1996 for (
unsigned int q = 0; q < n_q_points; ++q)
1998 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2000 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2002 transformed_shape_grads[q][d1] -=
2004 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
2009 for (
unsigned int q = 0; q < n_q_points; ++q)
2011 for (
unsigned int d = 0;
d < dim; ++
d)
2014 transformed_shape_grads[q][
d];
2023 template <
int dim,
int spacedim>
2027 const unsigned int face_no,
2031 const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2057 fill_edge_values(cell,
2060 if (dim == 3 && this->degree > 1)
2062 fill_face_values(cell,
2068 const unsigned int n_q_points = quadrature.
size();
2071 cell->face_orientation(face_no),
2072 cell->face_flip(face_no),
2073 cell->face_rotation(face_no),
2080 std::vector<Tensor<1, dim>> transformed_shape_values(n_q_points);
2081 for (
unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
2090 const unsigned int first =
2092 this->get_nonzero_components(dof)
2093 .first_selected_component()];
2095 for (
unsigned int q = 0; q < n_q_points; ++q)
2097 for (
unsigned int d = 0;
d < dim; ++
d)
2100 transformed_shape_values[q][
d];
2109 std::vector<Tensor<2, dim>> input(n_q_points);
2110 std::vector<Tensor<2, dim>> transformed_shape_grads(n_q_points);
2111 for (
unsigned int dof = 0; dof < this->dofs_per_cell; ++dof)
2113 for (
unsigned int q = 0; q < n_q_points; ++q)
2122 const unsigned int first =
2124 this->get_nonzero_components(dof)
2125 .first_selected_component()];
2127 for (
unsigned int q = 0; q < n_q_points; ++q)
2129 for (
unsigned int d1 = 0; d1 < dim; ++d1)
2131 for (
unsigned int d2 = 0; d2 < dim; ++d2)
2133 transformed_shape_grads[q][d1] -=
2135 mapping_data.jacobian_pushed_forward_grads[q][d2][d1];
2140 for (
unsigned int q = 0; q < n_q_points; ++q)
2142 for (
unsigned int d = 0;
d < dim; ++
d)
2145 transformed_shape_grads[q][
d];
2154 template <
int dim,
int spacedim>
2158 const unsigned int ,
2159 const unsigned int ,
2163 const ::internal::FEValuesImplementation::MappingRelatedData<dim, dim>
2174 template <
int dim,
int spacedim>
2201 template <
int dim,
int spacedim>
2212 std::ostringstream namebuf;
2213 namebuf <<
"FE_NedelecSZ<" << dim <<
">(" << this->degree - 1 <<
")";
2215 return namebuf.str();
2220 template <
int dim,
int spacedim>
2221 std::unique_ptr<FiniteElement<dim, dim>>
2224 return std_cxx14::make_unique<FE_NedelecSZ<dim, spacedim>>(*this);
2229 template <
int dim,
int spacedim>
2230 std::vector<unsigned int>
2239 std::vector<unsigned int> dpo(dim + 1);
2241 dpo[1] = degree + 1;
2242 dpo[2] = 2 * degree * (degree + 1);
2245 dpo[3] = 3 * degree * degree * (degree + 1);
2252 template <
int dim,
int spacedim>
2261 return 2 * (degree + 1) * (degree + 2);
2264 return 3 * (degree + 1) * (degree + 2) * (degree + 2);
2276 template <
int dim,
int spacedim>
2281 IntegratedLegendrePolynomials =
2288 #include "fe_nedelec_sz.inst"