From: Sebastian Kinnewig Date: Sun, 14 Apr 2024 13:43:25 +0000 (+0200) Subject: Move the computation of the InternalData into a separate function. X-Git-Tag: v9.6.0-rc1~368^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c58e0a39b3e963b73028b3853330e17c32f9fed;p=dealii.git Move the computation of the InternalData into a separate function. --- diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h index ee2852cece..d0508b2776 100644 --- a/include/deal.II/fe/fe_nedelec_sz.h +++ b/include/deal.II/fe/fe_nedelec_sz.h @@ -163,6 +163,17 @@ protected: */ MappingKind mapping_kind; + /** + * Compute the value and the derivatives of the Nedelec functions at + * the points given in p_list. + */ + void + evaluate(const std::vector> &p_list, + const UpdateFlags update_flags, + std::unique_ptr< + typename dealii::FiniteElement::InternalDataBase> + &data_ptr) const; + virtual std::unique_ptr< typename dealii::FiniteElement::InternalDataBase> get_data( diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc index d891c9631f..6f05014fa2 100644 --- a/source/fe/fe_nedelec_sz.cc +++ b/source/fe/fe_nedelec_sz.cc @@ -222,35 +222,43 @@ FE_NedelecSZ::shape_grad_grad_component( template -std::unique_ptr::InternalDataBase> -FE_NedelecSZ::get_data( - const UpdateFlags update_flags, - const Mapping & /*mapping*/, - const Quadrature &quadrature, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - & /*output_data*/) const -{ +void +FE_NedelecSZ::evaluate( + const std::vector> /*p_list*/, + const UpdateFlags /*update_flags*/, std::unique_ptr< typename dealii::FiniteElement::InternalDataBase> - data_ptr = std::make_unique(); + & /*data_ptr*/) const +{ + DEAL_II_NOT_IMPLEMENTED(); +} + + + +template <> +void +FE_NedelecSZ<2, 2>::evaluate( + const std::vector> &p_list, + const UpdateFlags update_flags, + std::unique_ptr::InternalDataBase> + &data_ptr) const +{ auto &data = dynamic_cast(*data_ptr); data.update_each = requires_update_flags(update_flags); // Useful quantities: + const unsigned int dim = 2; const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1 - const unsigned int vertices_per_cell = GeometryInfo::vertices_per_cell; - const unsigned int lines_per_cell = GeometryInfo::lines_per_cell; - const unsigned int faces_per_cell = GeometryInfo::faces_per_cell; + const unsigned int vertices_per_cell = + GeometryInfo<2 /*dim*/>::vertices_per_cell; + const unsigned int lines_per_cell = GeometryInfo<2 /*dim*/>::lines_per_cell; const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell; - // we assume that all quads have the same number of dofs - const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell; + const unsigned int n_q_points = p_list.size(); - const UpdateFlags flags(data.update_each); - const unsigned int n_q_points = quadrature.size(); + const UpdateFlags flags(data.update_each); // Resize the internal data storage: data.sigma_imj_values.resize( @@ -264,1302 +272,1151 @@ FE_NedelecSZ::get_data( // Resize shape function arrays according to update flags: if (flags & update_values) - { - data.shape_values.resize(this->n_dofs_per_cell(), - std::vector>(n_q_points)); - } + data.shape_values.resize(this->n_dofs_per_cell(), + std::vector>(n_q_points)); if (flags & update_gradients) + data.shape_grads.resize(this->n_dofs_per_cell(), + std::vector>( + n_q_points)); + + if (flags & update_hessians) + data.shape_hessians.resize(this->n_dofs_per_cell(), + std::vector>( + n_q_points)); + + // Compute values of sigma & lambda and the sigma differences and + // lambda additions. + std::vector> sigma(n_q_points, + std::vector(lines_per_cell)); + std::vector> lambda(n_q_points, + std::vector(lines_per_cell)); + + for (unsigned int q = 0; q < n_q_points; ++q) { - data.shape_grads.resize(this->n_dofs_per_cell(), - std::vector>( - n_q_points)); + sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]); + sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]); + sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1]; + sigma[q][3] = p_list[q][0] + p_list[q][1]; + + lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]); + lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]); + lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1]; + lambda[q][3] = p_list[q][0] * p_list[q][1]; + for (unsigned int i = 0; i < vertices_per_cell; ++i) + for (unsigned int j = 0; j < vertices_per_cell; ++j) + data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j]; } - if (flags & update_hessians) + // Calculate the gradient of sigma_imj_values[q][i][j] = + // sigma[q][i]-sigma[q][j] + // - this depends on the component and the direction of the + // corresponding edge. + // - the direction of the edge is determined by + // sigma_imj_sign[i][j]. + // Helper arrays: + const int sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1}, + {1, -1}, + {-1, 1}, + {1, 1}}; + int sigma_imj_sign[vertices_per_cell][vertices_per_cell]; + unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell]; + + for (unsigned int i = 0; i < vertices_per_cell; ++i) + for (unsigned int j = 0; j < vertices_per_cell; ++j) + { + // sigma_imj_sign is the sign (+/-) of the coefficient of + // x/y/z in sigma_imj_values Due to the numbering of vertices + // on the reference element it is easy to find edges in the + // positive direction are from smaller to higher local vertex + // numbering. + sigma_imj_sign[i][j] = (i < j) ? -1 : 1; + sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j]; + + // Now store the component which the sigma_i - sigma_j + // corresponds to: + sigma_imj_component[i][j] = 0; + for (unsigned int d = 0; d < dim; ++d) + { + int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d]; + // Only interested in the first non-zero + // as if there is a second, it can not be a valid edge. + if (temp_imj != 0) + { + sigma_imj_component[i][j] = d; + break; + } + } + // Can now calculate the gradient, only non-zero in the + // component given: Note some i,j combinations will be + // incorrect, but only on invalid edges. + data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] = + 2.0 * sigma_imj_sign[i][j]; + } + + // Now compute the edge parameterisations for a single element + // with global numbering matching that of the reference element: + + // Resize the edge parameterisations + data.edge_sigma_values.resize(lines_per_cell, + std::vector(n_q_points)); + data.edge_sigma_grads.resize(lines_per_cell, std::vector(dim)); + + // Fill the values for edge lambda and edge sigma: + const unsigned int edge_sigma_direction[lines_per_cell] = {1, 1, 0, 0}; + + data.edge_lambda_values.resize(lines_per_cell, + std::vector(n_q_points)); + + data.edge_lambda_grads_2d.resize(lines_per_cell, std::vector(dim)); + + for (unsigned int m = 0; m < lines_per_cell; ++m) { - data.shape_hessians.resize(this->n_dofs_per_cell(), - std::vector>( - n_q_points)); + // e1=max(reference vertex numbering on this edge) + // e2=min(reference vertex numbering on this edge) + // Which is guaranteed to be: + const unsigned int e1(GeometryInfo::line_to_cell_vertices(m, 1)); + const unsigned int e2(GeometryInfo::line_to_cell_vertices(m, 0)); + for (unsigned int q = 0; q < n_q_points; ++q) + { + data.edge_sigma_values[m][q] = data.sigma_imj_values[q][e2][e1]; + data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2]; + } + + data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0; } - std::vector> p_list(n_q_points); - p_list = quadrature.get_points(); + data.edge_lambda_grads_2d[0] = {-1.0, 0.0}; + data.edge_lambda_grads_2d[1] = {1.0, 0.0}; + data.edge_lambda_grads_2d[2] = {0.0, -1.0}; + data.edge_lambda_grads_2d[3] = {0.0, 1.0}; + // If the polynomial order is 0, then no more work to do: + if (degree < 1) + return; - switch (dim) + // Otherwise, we can compute the non-cell dependent shape functions. + // + // Note: the local dof numberings follow the usual order of lines -> + // faces -> cells + // (we have no vertex-based DoFs in this element). + // For a given cell we have: + // n_line_dofs = dofs_per_line*lines_per_cell. + // n_face_dofs = dofs_per_face*faces_per_cell. + // n_cell_dofs = dofs_per_quad (2d) + // = dofs_per_hex (3d) + // + // i.e. For the local dof numbering: + // the first line dof is 0, + // the first face dof is n_line_dofs, + // the first cell dof is n_line_dofs + n_face_dofs. + // + // On a line, DoFs are ordered first by line_dof and then line_index: + // i.e. line_dof_index = line_dof + line_index*(dofs_per_line) + // + // and similarly for faces: + // i.e. face_dof_index = face_dof + face_index*(dofs_per_face). + // + // HOWEVER, we have different types of DoFs on a line/face/cell. + // On a line we have two types, lowest order and higher order + // gradients. + // - The numbering is such the lowest order is first, then higher + // order. + // This is simple enough as there is only 1 lowest order and + // degree higher orders DoFs per line. + // + // On a 2d cell, we have 3 types: Type 1/2/3: + // - The ordering done by type: + // - Type 1: 0 <= i1,j1 < degree. degree^2 in total. + // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index + // = ij1. + // - Type 2: 0 <= i2,j2 < degree. degree^2 in total. + // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index + // = degree^2 + ij2 + // - Type 3: 0 <= i3 < 2*degree. 2*degree in total. + // Numbered: ij3 = i3. i.e. cell_dof_index + // = 2*(degree^2) + ij3. + // + // These then fit into the local dof numbering described above: + // - local dof numberings are: + // line_dofs: local_dof = line_dof_index. 0 <= local_dof < + // dofs_per_line*lines_per_cell face_dofs: local_dof = + // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof + // = n_lines_dof + n_face_dofs + cell_dof_index. + // + // The cell-based shape functions are: + // + // Type 1 (gradients): + // \phi^{C_{1}}_{ij} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ), + // + // 0 <= i,j < degree. + // + // NOTE: The derivative produced by IntegratedLegendrePolynomials does + // not account for the + // (2*x-1) or (2*y-1) so we must take this into account when + // taking derivatives. + const unsigned int cell_type1_offset = n_line_dofs; + + // Type 2: + // \phi^{C_{2}}_{ij} = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x} + // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y}, + // + // 0 <= i,j < degree. + const unsigned int cell_type2_offset = cell_type1_offset + degree * degree; + + // Type 3 (two subtypes): + // \phi^{C_{3}}_{j} = L_{j+2}(2y-1) \mathbf{e}_{x} + // + // \phi^{C_{3}}_{j+degree} = L_{j+2}(2x-1) \mathbf{e}_{y} + // + // 0 <= j < degree + const unsigned int cell_type3_offset1 = cell_type2_offset + degree * degree; + const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree; + + if (flags & (update_values | update_gradients | update_hessians)) { - case 2: + // compute all points we must evaluate the 1d polynomials at: + std::vector> cell_points(n_q_points); + for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int d = 0; d < dim; ++d) + cell_points[q][d] = 2.0 * p_list[q][d] - 1.0; + + // Loop through quad points: + for (unsigned int q = 0; q < n_q_points; ++q) { - // Compute values of sigma & lambda and the sigma differences and - // lambda additions. - std::vector> sigma( - n_q_points, std::vector(lines_per_cell)); - std::vector> lambda( - n_q_points, std::vector(lines_per_cell)); + // pre-compute values & required derivatives at this quad + // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1), + // + // for each polyc[d], c=x,y, contains the d-th derivative with + // respect to the coordinate c. - for (unsigned int q = 0; q < n_q_points; ++q) + // We only need poly values and 1st derivative for + // update_values, but need the 2nd derivative too for + // update_gradients. For update_hessians we also need the 3rd + // derivatives. + const unsigned int poly_length = + (flags & update_hessians) ? 4 : + ((flags & update_gradients) ? 3 : 2); + + std::vector> polyx( + degree, std::vector(poly_length)); + std::vector> polyy( + degree, std::vector(poly_length)); + for (unsigned int i = 0; i < degree; ++i) + { + // Compute all required 1d polynomials and their + // derivatives, starting at degree 2. e.g. to access + // L'_{3}(2x-1) use polyx[1][1]. + IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0], + polyx[i]); + IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1], + polyy[i]); + } + // Now use these to compute the shape functions: + if (flags & update_values) { - sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]); - sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]); - sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1]; - sigma[q][3] = p_list[q][0] + p_list[q][1]; - - lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]); - lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]); - lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1]; - lambda[q][3] = p_list[q][0] * p_list[q][1]; - for (unsigned int i = 0; i < vertices_per_cell; ++i) + for (unsigned int j = 0; j < degree; ++j) { - for (unsigned int j = 0; j < vertices_per_cell; ++j) + const unsigned int shift_j(j * degree); + for (unsigned int i = 0; i < degree; ++i) { - data.sigma_imj_values[q][i][j] = - sigma[q][i] - sigma[q][j]; + const unsigned int shift_ij(i + shift_j); + + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ij); + data.shape_values[dof_index1][q][0] = + 2.0 * polyx[i][1] * polyy[j][0]; + data.shape_values[dof_index1][q][1] = + 2.0 * polyx[i][0] * polyy[j][1]; + + // Type 2: + const unsigned int dof_index2(cell_type2_offset + + shift_ij); + data.shape_values[dof_index2][q][0] = + data.shape_values[dof_index1][q][0]; + data.shape_values[dof_index2][q][1] = + -1.0 * data.shape_values[dof_index1][q][1]; } + // Type 3: + const unsigned int dof_index3_1(cell_type3_offset1 + j); + data.shape_values[dof_index3_1][q][0] = polyy[j][0]; + data.shape_values[dof_index3_1][q][1] = 0.0; + + const unsigned int dof_index3_2(cell_type3_offset2 + j); + data.shape_values[dof_index3_2][q][0] = 0.0; + data.shape_values[dof_index3_2][q][1] = polyx[j][0]; } } - - // Calculate the gradient of sigma_imj_values[q][i][j] = - // sigma[q][i]-sigma[q][j] - // - this depends on the component and the direction of the - // corresponding edge. - // - the direction of the edge is determined by - // sigma_imj_sign[i][j]. - // Helper arrays: - const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = { - {-1, -1}, {1, -1}, {-1, 1}, {1, 1}}; - int sigma_imj_sign[vertices_per_cell][vertices_per_cell]; - unsigned int sigma_imj_component[vertices_per_cell] - [vertices_per_cell]; - - for (unsigned int i = 0; i < vertices_per_cell; ++i) + if (flags & update_gradients) { - for (unsigned int j = 0; j < vertices_per_cell; ++j) + for (unsigned int j = 0; j < degree; ++j) { - // sigma_imj_sign is the sign (+/-) of the coefficient of - // x/y/z in sigma_imj_values Due to the numbering of vertices - // on the reference element it is easy to find edges in the - // positive direction are from smaller to higher local vertex - // numbering. - sigma_imj_sign[i][j] = (i < j) ? -1 : 1; - sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j]; - - // Now store the component which the sigma_i - sigma_j - // corresponds to: - sigma_imj_component[i][j] = 0; - for (unsigned int d = 0; d < dim; ++d) + const unsigned int shift_j(j * degree); + for (unsigned int i = 0; i < degree; ++i) { - int temp_imj = - sigma_comp_signs[i][d] - sigma_comp_signs[j][d]; - // Only interested in the first non-zero - // as if there is a second, it can not be a valid edge. - if (temp_imj != 0) - { - sigma_imj_component[i][j] = d; - break; - } + const unsigned int shift_ij(i + shift_j); + + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ij); + data.shape_grads[dof_index1][q][0][0] = + 4.0 * polyx[i][2] * polyy[j][0]; + data.shape_grads[dof_index1][q][0][1] = + 4.0 * polyx[i][1] * polyy[j][1]; + data.shape_grads[dof_index1][q][1][0] = + data.shape_grads[dof_index1][q][0][1]; + data.shape_grads[dof_index1][q][1][1] = + 4.0 * polyx[i][0] * polyy[j][2]; + + // Type 2: + const unsigned int dof_index2(cell_type2_offset + + shift_ij); + data.shape_grads[dof_index2][q][0][0] = + data.shape_grads[dof_index1][q][0][0]; + data.shape_grads[dof_index2][q][0][1] = + data.shape_grads[dof_index1][q][0][1]; + data.shape_grads[dof_index2][q][1][0] = + -1.0 * data.shape_grads[dof_index1][q][1][0]; + data.shape_grads[dof_index2][q][1][1] = + -1.0 * data.shape_grads[dof_index1][q][1][1]; } - // Can now calculate the gradient, only non-zero in the - // component given: Note some i,j combinations will be - // incorrect, but only on invalid edges. - data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] = - 2.0 * sigma_imj_sign[i][j]; + // Type 3: + const unsigned int dof_index3_1(cell_type3_offset1 + j); + data.shape_grads[dof_index3_1][q][0][0] = 0.0; + data.shape_grads[dof_index3_1][q][0][1] = 2.0 * polyy[j][1]; + data.shape_grads[dof_index3_1][q][1][0] = 0.0; + data.shape_grads[dof_index3_1][q][1][1] = 0.0; + + const unsigned int dof_index3_2(cell_type3_offset2 + j); + data.shape_grads[dof_index3_2][q][0][0] = 0.0; + data.shape_grads[dof_index3_2][q][0][1] = 0.0; + data.shape_grads[dof_index3_2][q][1][0] = 2.0 * polyx[j][1]; + data.shape_grads[dof_index3_2][q][1][1] = 0.0; } } - - // Now compute the edge parameterisations for a single element - // with global numbering matching that of the reference element: - - // Resize the edge parameterisations - data.edge_sigma_values.resize(lines_per_cell); - data.edge_sigma_grads.resize(lines_per_cell); - for (unsigned int m = 0; m < lines_per_cell; ++m) + if (flags & update_hessians) { - data.edge_sigma_values[m].resize(n_q_points); + for (unsigned int j = 0; j < degree; ++j) + { + const unsigned int shift_j(j * degree); + for (unsigned int i = 0; i < degree; ++i) + { + const unsigned int shift_ij(i + shift_j); - // sigma grads are constant in a cell (no need for quad points) - data.edge_sigma_grads[m].resize(dim); - } + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ij); + data.shape_hessians[dof_index1][q][0][0][0] = + 8.0 * polyx[i][3] * polyy[j][0]; + data.shape_hessians[dof_index1][q][1][0][0] = + 8.0 * polyx[i][2] * polyy[j][1]; - // Fill the values for edge lambda and edge sigma: - const unsigned int - edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1, - 1, - 0, - 0}; - - data.edge_lambda_values.resize(lines_per_cell, - std::vector(n_q_points)); - data.edge_lambda_grads_2d.resize(lines_per_cell, - std::vector(dim)); - for (unsigned int m = 0; m < lines_per_cell; ++m) - { - // e1=max(reference vertex numbering on this edge) - // e2=min(reference vertex numbering on this edge) - // Which is guaranteed to be: - const unsigned int e1( - GeometryInfo::line_to_cell_vertices(m, 1)); - const unsigned int e2( - GeometryInfo::line_to_cell_vertices(m, 0)); - for (unsigned int q = 0; q < n_q_points; ++q) - { - data.edge_sigma_values[m][q] = - data.sigma_imj_values[q][e2][e1]; - data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2]; - } + data.shape_hessians[dof_index1][q][0][1][0] = + data.shape_hessians[dof_index1][q][1][0][0]; + data.shape_hessians[dof_index1][q][1][1][0] = + 8.0 * polyx[i][1] * polyy[j][2]; - data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0; - } - data.edge_lambda_grads_2d[0] = {-1.0, 0.0}; - data.edge_lambda_grads_2d[1] = {1.0, 0.0}; - data.edge_lambda_grads_2d[2] = {0.0, -1.0}; - data.edge_lambda_grads_2d[3] = {0.0, 1.0}; + data.shape_hessians[dof_index1][q][0][0][1] = + data.shape_hessians[dof_index1][q][1][0][0]; + data.shape_hessians[dof_index1][q][1][0][1] = + data.shape_hessians[dof_index1][q][1][1][0]; - // If the polynomial order is 0, then no more work to do: - if (degree < 1) - { - break; - } + data.shape_hessians[dof_index1][q][0][1][1] = + data.shape_hessians[dof_index1][q][1][1][0]; + data.shape_hessians[dof_index1][q][1][1][1] = + 8.0 * polyx[i][0] * polyy[j][3]; - // Otherwise, we can compute the non-cell dependent shape functions. - // - // Note: the local dof numberings follow the usual order of lines -> - // faces -> cells - // (we have no vertex-based DoFs in this element). - // For a given cell we have: - // n_line_dofs = dofs_per_line*lines_per_cell. - // n_face_dofs = dofs_per_face*faces_per_cell. - // n_cell_dofs = dofs_per_quad (2d) - // = dofs_per_hex (3d) - // - // i.e. For the local dof numbering: - // the first line dof is 0, - // the first face dof is n_line_dofs, - // the first cell dof is n_line_dofs + n_face_dofs. - // - // On a line, DoFs are ordered first by line_dof and then line_index: - // i.e. line_dof_index = line_dof + line_index*(dofs_per_line) - // - // and similarly for faces: - // i.e. face_dof_index = face_dof + face_index*(dofs_per_face). - // - // HOWEVER, we have different types of DoFs on a line/face/cell. - // On a line we have two types, lowest order and higher order - // gradients. - // - The numbering is such the lowest order is first, then higher - // order. - // This is simple enough as there is only 1 lowest order and - // degree higher orders DoFs per line. - // - // On a 2d cell, we have 3 types: Type 1/2/3: - // - The ordering done by type: - // - Type 1: 0 <= i1,j1 < degree. degree^2 in total. - // Numbered: ij1 = i1 + j1*(degree). i.e. cell_dof_index - // = ij1. - // - Type 2: 0 <= i2,j2 < degree. degree^2 in total. - // Numbered: ij2 = i2 + j2*(degree). i.e. cell_dof_index - // = degree^2 + ij2 - // - Type 3: 0 <= i3 < 2*degree. 2*degree in total. - // Numbered: ij3 = i3. i.e. cell_dof_index - // = 2*(degree^2) + ij3. - // - // These then fit into the local dof numbering described above: - // - local dof numberings are: - // line_dofs: local_dof = line_dof_index. 0 <= local_dof < - // dofs_per_line*lines_per_cell face_dofs: local_dof = - // n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof - // = n_lines_dof + n_face_dofs + cell_dof_index. - // - // The cell-based shape functions are: - // - // Type 1 (gradients): - // \phi^{C_{1}}_{ij) = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ), - // - // 0 <= i,j < degree. - // - // NOTE: The derivative produced by IntegratedLegendrePolynomials does - // not account for the - // (2*x-1) or (2*y-1) so we must take this into account when - // taking derivatives. - const unsigned int cell_type1_offset = n_line_dofs; - - // Type 2: - // \phi^{C_{2}}_{ij) = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x} - // - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y}, - // - // 0 <= i,j < degree. - const unsigned int cell_type2_offset = - cell_type1_offset + degree * degree; - // Type 3 (two subtypes): - // \phi^{C_{3}}_{j) = L_{j+2}(2y-1) \mathbf{e}_{x} - // - // \phi^{C_{3}}_{j+degree) = L_{j+2}(2x-1) \mathbf{e}_{y} - // - // 0 <= j < degree - const unsigned int cell_type3_offset1 = - cell_type2_offset + degree * degree; - const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree; - if (flags & (update_values | update_gradients | update_hessians)) - { - // compute all points we must evaluate the 1d polynomials at: - std::vector> cell_points(n_q_points); - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int d = 0; d < dim; ++d) - { - cell_points[q][d] = 2.0 * p_list[q][d] - 1.0; + // Type 2: + const unsigned int dof_index2(cell_type2_offset + + shift_ij); + for (unsigned int d = 0; d < dim; ++d) + { + data.shape_hessians[dof_index2][q][0][0][d] = + data.shape_hessians[dof_index1][q][0][0][d]; + data.shape_hessians[dof_index2][q][0][1][d] = + data.shape_hessians[dof_index1][q][0][1][d]; + data.shape_hessians[dof_index2][q][1][0][d] = + -1.0 * data.shape_hessians[dof_index1][q][1][0][d]; + data.shape_hessians[dof_index2][q][1][1][d] = + -1.0 * data.shape_hessians[dof_index1][q][1][1][d]; + } } + // Type 3: + const unsigned int dof_index3_1(cell_type3_offset1 + j); + data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0; + data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0; + data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0; + data.shape_hessians[dof_index3_1][q][0][1][1] = + 4.0 * polyy[j][2]; + data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0; + data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0; + data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0; + data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0; + + const unsigned int dof_index3_2(cell_type3_offset2 + j); + data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0; + data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0; + data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0; + data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0; + data.shape_hessians[dof_index3_2][q][1][0][0] = + 4.0 * polyx[j][2]; + data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0; + data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0; + data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0; } + } + } + } +} - // Loop through quad points: - for (unsigned int q = 0; q < n_q_points; ++q) - { - // pre-compute values & required derivatives at this quad - // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1), - // - // for each polyc[d], c=x,y, contains the d-th derivative with - // respect to the coordinate c. - - // We only need poly values and 1st derivative for - // update_values, but need the 2nd derivative too for - // update_gradients. For update_hessians we also need the 3rd - // derivatives. - const unsigned int poly_length = - (flags & update_hessians) ? - 4 : - ((flags & update_gradients) ? 3 : 2); - - std::vector> polyx( - degree, std::vector(poly_length)); - std::vector> polyy( - degree, std::vector(poly_length)); - for (unsigned int i = 0; i < degree; ++i) - { - // Compute all required 1d polynomials and their - // derivatives, starting at degree 2. e.g. to access - // L'_{3}(2x-1) use polyx[1][1]. - IntegratedLegendrePolynomials[i + 2].value( - cell_points[q][0], polyx[i]); - IntegratedLegendrePolynomials[i + 2].value( - cell_points[q][1], polyy[i]); - } - // Now use these to compute the shape functions: - if (flags & update_values) - { - for (unsigned int j = 0; j < degree; ++j) - { - const unsigned int shift_j(j * degree); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ij(i + shift_j); - // Type 1: - const unsigned int dof_index1(cell_type1_offset + - shift_ij); - data.shape_values[dof_index1][q][0] = - 2.0 * polyx[i][1] * polyy[j][0]; - data.shape_values[dof_index1][q][1] = - 2.0 * polyx[i][0] * polyy[j][1]; - // Type 2: - const unsigned int dof_index2(cell_type2_offset + - shift_ij); - data.shape_values[dof_index2][q][0] = - data.shape_values[dof_index1][q][0]; - data.shape_values[dof_index2][q][1] = - -1.0 * data.shape_values[dof_index1][q][1]; - } - // Type 3: - const unsigned int dof_index3_1(cell_type3_offset1 + - j); - data.shape_values[dof_index3_1][q][0] = polyy[j][0]; - data.shape_values[dof_index3_1][q][1] = 0.0; +template <> +void +FE_NedelecSZ<3, 3>::evaluate( + const std::vector> p_list, + const UpdateFlags update_flags, + std::unique_ptr::InternalDataBase> + &data_ptr) const +{ + auto &data = dynamic_cast(*data_ptr); + data.update_each = requires_update_flags(update_flags); - const unsigned int dof_index3_2(cell_type3_offset2 + - j); - data.shape_values[dof_index3_2][q][0] = 0.0; - data.shape_values[dof_index3_2][q][1] = polyx[j][0]; - } - } - if (flags & update_gradients) - { - for (unsigned int j = 0; j < degree; ++j) - { - const unsigned int shift_j(j * degree); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ij(i + shift_j); + // Useful quantities: + const unsigned int dim = 3; + const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1 - // Type 1: - const unsigned int dof_index1(cell_type1_offset + - shift_ij); - data.shape_grads[dof_index1][q][0][0] = - 4.0 * polyx[i][2] * polyy[j][0]; - data.shape_grads[dof_index1][q][0][1] = - 4.0 * polyx[i][1] * polyy[j][1]; - data.shape_grads[dof_index1][q][1][0] = - data.shape_grads[dof_index1][q][0][1]; - data.shape_grads[dof_index1][q][1][1] = - 4.0 * polyx[i][0] * polyy[j][2]; + const unsigned int vertices_per_cell = + GeometryInfo<3 /*dim*/>::vertices_per_cell; + const unsigned int lines_per_cell = GeometryInfo<3 /*dim*/>::lines_per_cell; + const unsigned int faces_per_cell = GeometryInfo<3 /*dim*/>::faces_per_cell; - // Type 2: - const unsigned int dof_index2(cell_type2_offset + - shift_ij); - data.shape_grads[dof_index2][q][0][0] = - data.shape_grads[dof_index1][q][0][0]; - data.shape_grads[dof_index2][q][0][1] = - data.shape_grads[dof_index1][q][0][1]; - data.shape_grads[dof_index2][q][1][0] = - -1.0 * data.shape_grads[dof_index1][q][1][0]; - data.shape_grads[dof_index2][q][1][1] = - -1.0 * data.shape_grads[dof_index1][q][1][1]; - } - // Type 3: - const unsigned int dof_index3_1(cell_type3_offset1 + - j); - data.shape_grads[dof_index3_1][q][0][0] = 0.0; - data.shape_grads[dof_index3_1][q][0][1] = - 2.0 * polyy[j][1]; - data.shape_grads[dof_index3_1][q][1][0] = 0.0; - data.shape_grads[dof_index3_1][q][1][1] = 0.0; + const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell; - const unsigned int dof_index3_2(cell_type3_offset2 + - j); - data.shape_grads[dof_index3_2][q][0][0] = 0.0; - data.shape_grads[dof_index3_2][q][0][1] = 0.0; - data.shape_grads[dof_index3_2][q][1][0] = - 2.0 * polyx[j][1]; - data.shape_grads[dof_index3_2][q][1][1] = 0.0; - } - } - if (flags & update_hessians) - { - for (unsigned int j = 0; j < degree; ++j) - { - const unsigned int shift_j(j * degree); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ij(i + shift_j); + // we assume that all quads have the same number of dofs + const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell; + const unsigned int n_q_points = p_list.size(); - // Type 1: - const unsigned int dof_index1(cell_type1_offset + - shift_ij); - data.shape_hessians[dof_index1][q][0][0][0] = - 8.0 * polyx[i][3] * polyy[j][0]; - data.shape_hessians[dof_index1][q][1][0][0] = - 8.0 * polyx[i][2] * polyy[j][1]; + const UpdateFlags flags(data.update_each); - data.shape_hessians[dof_index1][q][0][1][0] = - data.shape_hessians[dof_index1][q][1][0][0]; - data.shape_hessians[dof_index1][q][1][1][0] = - 8.0 * polyx[i][1] * polyy[j][2]; + // Resize the internal data storage: + data.sigma_imj_values.resize( + n_q_points, + std::vector>(vertices_per_cell, + std::vector(vertices_per_cell))); - data.shape_hessians[dof_index1][q][0][0][1] = - data.shape_hessians[dof_index1][q][1][0][0]; - data.shape_hessians[dof_index1][q][1][0][1] = - data.shape_hessians[dof_index1][q][1][1][0]; + data.sigma_imj_grads.resize(vertices_per_cell, + std::vector>( + vertices_per_cell, std::vector(dim))); - data.shape_hessians[dof_index1][q][0][1][1] = - data.shape_hessians[dof_index1][q][1][1][0]; - data.shape_hessians[dof_index1][q][1][1][1] = - 8.0 * polyx[i][0] * polyy[j][3]; + // Resize shape function arrays according to update flags: + if (flags & update_values) + data.shape_values.resize(this->n_dofs_per_cell(), + std::vector>(n_q_points)); + if (flags & update_gradients) + data.shape_grads.resize(this->n_dofs_per_cell(), + std::vector>( + n_q_points)); + if (flags & update_hessians) + data.shape_hessians.resize(this->n_dofs_per_cell(), + std::vector>( + n_q_points)); + + // Compute values of sigma & lambda and the sigma differences and + // lambda additions. + std::vector> sigma(n_q_points, + std::vector(lines_per_cell)); + std::vector> lambda(n_q_points, + std::vector(lines_per_cell)); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + sigma[q][0] = + (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + (1 - p_list[q][2]); + sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]); + sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]); + sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]); + sigma[q][4] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2]; + sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2]; + sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2]; + sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2]; + + lambda[q][0] = + (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]); + lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]); + lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]); + lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]); + lambda[q][4] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2]; + lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2]; + lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2]; + lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2]; + + // Compute values of sigma_imj = \sigma_{i} - \sigma_{j} + // and lambda_ipj = \lambda_{i} + \lambda_{j}. + for (unsigned int i = 0; i < vertices_per_cell; ++i) + for (unsigned int j = 0; j < vertices_per_cell; ++j) + data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j]; + } - // Type 2: - const unsigned int dof_index2(cell_type2_offset + - shift_ij); - for (unsigned int d = 0; d < dim; ++d) - { - data.shape_hessians[dof_index2][q][0][0][d] = - data.shape_hessians[dof_index1][q][0][0][d]; - data.shape_hessians[dof_index2][q][0][1][d] = - data.shape_hessians[dof_index1][q][0][1][d]; - data.shape_hessians[dof_index2][q][1][0][d] = - -1.0 * - data.shape_hessians[dof_index1][q][1][0][d]; - data.shape_hessians[dof_index2][q][1][1][d] = - -1.0 * - data.shape_hessians[dof_index1][q][1][1][d]; - } - } - // Type 3: - const unsigned int dof_index3_1(cell_type3_offset1 + - j); - data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0; - data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0; - data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0; - data.shape_hessians[dof_index3_1][q][0][1][1] = - 4.0 * polyy[j][2]; - data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0; - data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0; - data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0; - data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0; - - const unsigned int dof_index3_2(cell_type3_offset2 + - j); - data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0; - data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0; - data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0; - data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0; - data.shape_hessians[dof_index3_2][q][1][0][0] = - 4.0 * polyx[j][2]; - data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0; - data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0; - data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0; - } - } - } - } - break; + // We now want some additional information about + // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to + // calculate values & derivatives of the shape functions we need to + // know: + // - The component the sigma_imj value corresponds to - this varies + // with i & j. + // - The gradient of the sigma_imj value + // - this depends on the component and the direction of the + // corresponding edge. + // - the direction of the edge is determined by + // sigma_imj_sign[i][j]. + // + // Note that not every i,j combination is a valid edge (there are only + // 12 valid edges in 3d), but we compute them all as it simplifies + // things. + + // store the sign of each component x, y, z in the sigma list. + // can use this to fill in the sigma_imj_component data. + const int sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1, -1}, + {1, -1, -1}, + {-1, 1, -1}, + {1, 1, -1}, + {-1, -1, 1}, + {1, -1, 1}, + {-1, 1, 1}, + {1, 1, 1}}; + + int sigma_imj_sign[vertices_per_cell][vertices_per_cell]; + unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell]; + + for (unsigned int i = 0; i < vertices_per_cell; ++i) + for (unsigned int j = 0; j < vertices_per_cell; ++j) + { + // sigma_imj_sign is the sign (+/-) of the coefficient of + // x/y/z in sigma_imj. Due to the numbering of vertices on the + // reference element this is easy to work out because edges in + // the positive direction go from smaller to higher local + // vertex numbering. + sigma_imj_sign[i][j] = (i < j) ? -1 : 1; + sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j]; + + // Now store the component which the sigma_i - sigma_j + // corresponds to: + sigma_imj_component[i][j] = 0; + for (unsigned int d = 0; d < dim; ++d) + { + int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d]; + // Only interested in the first non-zero + // as if there is a second, it will not be a valid edge. + if (temp_imj != 0) + { + sigma_imj_component[i][j] = d; + break; + } + } + // Can now calculate the gradient, only non-zero in the + // component given: Note some i,j combinations will be + // incorrect, but only on invalid edges. + data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] = + 2.0 * sigma_imj_sign[i][j]; + } + + // Now compute the edge parameterisations for a single element + // with global numbering matching that of the reference element: + + // resize the edge parameterisations + data.edge_sigma_values.resize(lines_per_cell, + std::vector(n_q_points)); + data.edge_lambda_values.resize(lines_per_cell, + std::vector(n_q_points)); + data.edge_sigma_grads.resize(lines_per_cell, std::vector(dim)); + data.edge_lambda_grads_3d.resize( + lines_per_cell, + std::vector>(n_q_points, std::vector(dim))); + data.edge_lambda_gradgrads_3d.resize( + lines_per_cell, + std::vector>(dim, std::vector(dim))); + + // Fill the values: + const unsigned int edge_sigma_direction[lines_per_cell] = { + 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2}; + + for (unsigned int m = 0; m < lines_per_cell; ++m) + { + // e1=max(reference vertex numbering on this edge) + // e2=min(reference vertex numbering on this edge) + // Which is guaranteed to be: + const unsigned int e1(GeometryInfo::line_to_cell_vertices(m, 1)); + const unsigned int e2(GeometryInfo::line_to_cell_vertices(m, 0)); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + data.edge_sigma_values[m][q] = data.sigma_imj_values[q][e2][e1]; + data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2]; } - case 3: + + data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0; + } + + // edge_lambda_grads + for (unsigned int q = 0; q < n_q_points; ++q) + { + double x(p_list[q][0]); + double y(p_list[q][1]); + double z(p_list[q][2]); + data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0}; + data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x}; + data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0}; + data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y}; + data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x}; + data.edge_lambda_grads_3d[5][q] = {z, 0.0, x}; + data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y}; + data.edge_lambda_grads_3d[7][q] = {0.0, z, y}; + data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0}; + data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0}; + data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0}; + data.edge_lambda_grads_3d[11][q] = {y, x, 0.0}; + } + + // edge_lambda gradgrads: + const int edge_lambda_sign[lines_per_cell] = { + 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1}; // sign of the 2nd derivative for + // each edge. + + const unsigned int edge_lambda_directions[lines_per_cell][2] = { + {0, 2}, + {0, 2}, + {1, 2}, + {1, 2}, + {0, 2}, + {0, 2}, + {1, 2}, + {1, 2}, + {0, 1}, + {0, 1}, + {0, 1}, + {0, 1}}; // component which edge_lambda[m] depends on. + + for (unsigned int m = 0; m < lines_per_cell; ++m) + { + data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]] + [edge_lambda_directions[m][1]] = + edge_lambda_sign[m]; + data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]] + [edge_lambda_directions[m][0]] = + edge_lambda_sign[m]; + } + + // If the polynomial order is 0, then no more work to do: + if (degree < 1) + return; + + // resize required data: + data.face_lambda_values.resize(faces_per_cell, + std::vector(n_q_points)); + data.face_lambda_grads.resize(faces_per_cell, std::vector(dim)); + + // Fill in the values (these don't change between cells). + for (unsigned int q = 0; q < n_q_points; ++q) + { + double x(p_list[q][0]); + double y(p_list[q][1]); + double z(p_list[q][2]); + data.face_lambda_values[0][q] = 1.0 - x; + data.face_lambda_values[1][q] = x; + data.face_lambda_values[2][q] = 1.0 - y; + data.face_lambda_values[3][q] = y; + data.face_lambda_values[4][q] = 1.0 - z; + data.face_lambda_values[5][q] = z; + } + + // gradients are constant: + data.face_lambda_grads[0] = {-1.0, 0.0, 0.0}; + data.face_lambda_grads[1] = {1.0, 0.0, 0.0}; + data.face_lambda_grads[2] = {0.0, -1.0, 0.0}; + data.face_lambda_grads[3] = {0.0, 1.0, 0.0}; + data.face_lambda_grads[4] = {0.0, 0.0, -1.0}; + data.face_lambda_grads[5] = {0.0, 0.0, 1.0}; + + // for cell-based shape functions: + // these don't depend on the cell, so can precompute all here: + if (flags & (update_values | update_gradients | update_hessians)) + { + // Cell-based shape functions: + // + // Type-1 (gradients): + // \phi^{C_{1}}_{ijk} = grad( + // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ), + // + // 0 <= i,j,k < degree. (in a group of degree*degree*degree) + const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs); + // Type-2: + // + // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk} + // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1, + // -1)\phi^{C_{1}}_{ijk} + // + // 0 <= i,j,k < degree. (subtypes in groups of + // degree*degree*degree) + // + // here we order so that all of subtype 1 comes first, then + // subtype 2. + const unsigned int cell_type2_offset1(cell_type1_offset + + degree * degree * degree); + const unsigned int cell_type2_offset2(cell_type2_offset1 + + degree * degree * degree); + // Type-3 + // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x} + // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y} + // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z} + // + // 0 <= i,j,k < degree. (subtypes in groups of degree*degree) + // + // again we order so we compute all of subtype 1 first, then + // subtype 2, etc. + const unsigned int cell_type3_offset1(cell_type2_offset2 + + degree * degree * degree); + const unsigned int cell_type3_offset2(cell_type3_offset1 + + degree * degree); + const unsigned int cell_type3_offset3(cell_type3_offset2 + + degree * degree); + + // compute all points we must evaluate the 1d polynomials at: + std::vector> cell_points(n_q_points); + for (unsigned int q = 0; q < n_q_points; ++q) { - // Compute values of sigma & lambda and the sigma differences and - // lambda additions. - std::vector> sigma( - n_q_points, std::vector(lines_per_cell)); - std::vector> lambda( - n_q_points, std::vector(lines_per_cell)); - for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int d = 0; d < dim; ++d) { - sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + - (1 - p_list[q][2]); - sigma[q][1] = - p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]); - sigma[q][2] = - (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]); - sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]); - sigma[q][4] = - (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2]; - sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2]; - sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2]; - sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2]; - - lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * - (1.0 - p_list[q][2]); - lambda[q][1] = - p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]); - lambda[q][2] = - (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]); - lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]); - lambda[q][4] = - (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2]; - lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2]; - lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2]; - lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2]; - - // Compute values of sigma_imj = \sigma_{i} - \sigma_{j} - // and lambda_ipj = \lambda_{i} + \lambda_{j}. - for (unsigned int i = 0; i < vertices_per_cell; ++i) - { - for (unsigned int j = 0; j < vertices_per_cell; ++j) - { - data.sigma_imj_values[q][i][j] = - sigma[q][i] - sigma[q][j]; - } - } + cell_points[q][d] = 2.0 * p_list[q][d] - 1.0; } + } + + // We only need poly values and 1st derivative for + // update_values, but need the 2nd derivative too for + // update_gradients. For update_hessians we also need 3rd + // derivative. + const unsigned int poly_length = + (flags & update_hessians) ? 4 : ((flags & update_gradients) ? 3 : 2); - // We now want some additional information about - // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to - // calculate values & derivatives of the shape functions we need to - // know: - // - The component the sigma_imj value corresponds to - this varies - // with i & j. - // - The gradient of the sigma_imj value - // - this depends on the component and the direction of the - // corresponding edge. - // - the direction of the edge is determined by - // sigma_imj_sign[i][j]. + // Loop through quad points: + for (unsigned int q = 0; q < n_q_points; ++q) + { + // pre-compute values & required derivatives at this quad + // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy = + // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1). // - // Note that not every i,j combination is a valid edge (there are only - // 12 valid edges in 3d), but we compute them all as it simplifies - // things. - - // store the sign of each component x, y, z in the sigma list. - // can use this to fill in the sigma_imj_component data. - const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = { - {-1, -1, -1}, - {1, -1, -1}, - {-1, 1, -1}, - {1, 1, -1}, - {-1, -1, 1}, - {1, -1, 1}, - {-1, 1, 1}, - {1, 1, 1}}; - - int sigma_imj_sign[vertices_per_cell][vertices_per_cell]; - unsigned int sigma_imj_component[vertices_per_cell] - [vertices_per_cell]; - - for (unsigned int i = 0; i < vertices_per_cell; ++i) + // for each polyc[d], c=x,y,z, contains the d-th + // derivative with respect to the coordinate c. + std::vector> polyx( + degree, std::vector(poly_length)); + std::vector> polyy( + degree, std::vector(poly_length)); + std::vector> polyz( + degree, std::vector(poly_length)); + for (unsigned int i = 0; i < degree; ++i) + { + // compute all required 1d polynomials for i + IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0], + polyx[i]); + IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1], + polyy[i]); + IntegratedLegendrePolynomials[i + 2].value(cell_points[q][2], + polyz[i]); + } + // Now use these to compute the shape functions: + if (flags & update_values) { - for (unsigned int j = 0; j < vertices_per_cell; ++j) + for (unsigned int k = 0; k < degree; ++k) { - // sigma_imj_sign is the sign (+/-) of the coefficient of - // x/y/z in sigma_imj. Due to the numbering of vertices on the - // reference element this is easy to work out because edges in - // the positive direction go from smaller to higher local - // vertex numbering. - sigma_imj_sign[i][j] = (i < j) ? -1 : 1; - sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j]; - - // Now store the component which the sigma_i - sigma_j - // corresponds to: - sigma_imj_component[i][j] = 0; - for (unsigned int d = 0; d < dim; ++d) + const unsigned int shift_k(k * degree * degree); + const unsigned int shift_j(k * + degree); // Used below when subbing + // k for j (type 3) + for (unsigned int j = 0; j < degree; ++j) { - int temp_imj = - sigma_comp_signs[i][d] - sigma_comp_signs[j][d]; - // Only interested in the first non-zero - // as if there is a second, it will not be a valid edge. - if (temp_imj != 0) + const unsigned int shift_jk(j * degree + shift_k); + for (unsigned int i = 0; i < degree; ++i) { - sigma_imj_component[i][j] = d; - break; + const unsigned int shift_ijk(shift_jk + i); + + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ijk); + + data.shape_values[dof_index1][q][0] = + 2.0 * polyx[i][1] * polyy[j][0] * polyz[k][0]; + data.shape_values[dof_index1][q][1] = + 2.0 * polyx[i][0] * polyy[j][1] * polyz[k][0]; + data.shape_values[dof_index1][q][2] = + 2.0 * polyx[i][0] * polyy[j][0] * polyz[k][1]; + + // Type 2: + const unsigned int dof_index2_1(cell_type2_offset1 + + shift_ijk); + const unsigned int dof_index2_2(cell_type2_offset2 + + shift_ijk); + + data.shape_values[dof_index2_1][q][0] = + data.shape_values[dof_index1][q][0]; + data.shape_values[dof_index2_1][q][1] = + -1.0 * data.shape_values[dof_index1][q][1]; + data.shape_values[dof_index2_1][q][2] = + data.shape_values[dof_index1][q][2]; + + data.shape_values[dof_index2_2][q][0] = + data.shape_values[dof_index1][q][0]; + data.shape_values[dof_index2_2][q][1] = + -1.0 * data.shape_values[dof_index1][q][1]; + data.shape_values[dof_index2_2][q][2] = + -1.0 * data.shape_values[dof_index1][q][2]; } + // Type 3: (note we re-use k and j for + // convenience): + const unsigned int shift_ij(j + + shift_j); // here we've subbed + // j for i, k for j. + const unsigned int dof_index3_1(cell_type3_offset1 + + shift_ij); + const unsigned int dof_index3_2(cell_type3_offset2 + + shift_ij); + const unsigned int dof_index3_3(cell_type3_offset3 + + shift_ij); + + data.shape_values[dof_index3_1][q][0] = + polyy[j][0] * polyz[k][0]; + data.shape_values[dof_index3_1][q][1] = 0.0; + data.shape_values[dof_index3_1][q][2] = 0.0; + + data.shape_values[dof_index3_2][q][0] = 0.0; + data.shape_values[dof_index3_2][q][1] = + polyx[j][0] * polyz[k][0]; + data.shape_values[dof_index3_2][q][2] = 0.0; + + data.shape_values[dof_index3_3][q][0] = 0.0; + data.shape_values[dof_index3_3][q][1] = 0.0; + data.shape_values[dof_index3_3][q][2] = + polyx[j][0] * polyy[k][0]; } - // Can now calculate the gradient, only non-zero in the - // component given: Note some i,j combinations will be - // incorrect, but only on invalid edges. - data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] = - 2.0 * sigma_imj_sign[i][j]; } } - // Now compute the edge parameterisations for a single element - // with global numbering matching that of the reference element: - - // resize the edge parameterisations - data.edge_sigma_values.resize(lines_per_cell); - data.edge_lambda_values.resize(lines_per_cell); - data.edge_sigma_grads.resize(lines_per_cell); - data.edge_lambda_grads_3d.resize(lines_per_cell); - data.edge_lambda_gradgrads_3d.resize(lines_per_cell); - for (unsigned int m = 0; m < lines_per_cell; ++m) + if (flags & update_gradients) { - data.edge_sigma_values[m].resize(n_q_points); - data.edge_lambda_values[m].resize(n_q_points); - - // sigma grads are constant in a cell (no need for quad points) - data.edge_sigma_grads[m].resize(dim); - - data.edge_lambda_grads_3d[m].resize(n_q_points); - for (unsigned int q = 0; q < n_q_points; ++q) - { - data.edge_lambda_grads_3d[m][q].resize(dim); - } - // lambda_gradgrads are constant in a cell (no need for quad - // points) - data.edge_lambda_gradgrads_3d[m].resize(dim); - for (unsigned int d = 0; d < dim; ++d) + for (unsigned int k = 0; k < degree; ++k) { - data.edge_lambda_gradgrads_3d[m][d].resize(dim); - } - } - - // Fill the values: - const unsigned int - edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = { - 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2}; - - for (unsigned int m = 0; m < lines_per_cell; ++m) - { - // e1=max(reference vertex numbering on this edge) - // e2=min(reference vertex numbering on this edge) - // Which is guaranteed to be: - const unsigned int e1( - GeometryInfo::line_to_cell_vertices(m, 1)); - const unsigned int e2( - GeometryInfo::line_to_cell_vertices(m, 0)); + const unsigned int shift_k(k * degree * degree); + const unsigned int shift_j(k * + degree); // Used below when subbing + // k for j (type 3) + for (unsigned int j = 0; j < degree; ++j) + { + const unsigned int shift_jk(j * degree + shift_k); + for (unsigned int i = 0; i < degree; ++i) + { + const unsigned int shift_ijk(shift_jk + i); + + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ijk); + + data.shape_grads[dof_index1][q][0][0] = + 4.0 * polyx[i][2] * polyy[j][0] * polyz[k][0]; + data.shape_grads[dof_index1][q][0][1] = + 4.0 * polyx[i][1] * polyy[j][1] * polyz[k][0]; + data.shape_grads[dof_index1][q][0][2] = + 4.0 * polyx[i][1] * polyy[j][0] * polyz[k][1]; + + data.shape_grads[dof_index1][q][1][0] = + data.shape_grads[dof_index1][q][0][1]; + data.shape_grads[dof_index1][q][1][1] = + 4.0 * polyx[i][0] * polyy[j][2] * polyz[k][0]; + data.shape_grads[dof_index1][q][1][2] = + 4.0 * polyx[i][0] * polyy[j][1] * polyz[k][1]; + + data.shape_grads[dof_index1][q][2][0] = + data.shape_grads[dof_index1][q][0][2]; + data.shape_grads[dof_index1][q][2][1] = + data.shape_grads[dof_index1][q][1][2]; + data.shape_grads[dof_index1][q][2][2] = + 4.0 * polyx[i][0] * polyy[j][0] * polyz[k][2]; + + // Type 2: + const unsigned int dof_index2_1(cell_type2_offset1 + + shift_ijk); + const unsigned int dof_index2_2(cell_type2_offset2 + + shift_ijk); - for (unsigned int q = 0; q < n_q_points; ++q) - { - data.edge_sigma_values[m][q] = - data.sigma_imj_values[q][e2][e1]; - data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2]; + for (unsigned int d = 0; d < dim; ++d) + { + data.shape_grads[dof_index2_1][q][0][d] = + data.shape_grads[dof_index1][q][0][d]; + data.shape_grads[dof_index2_1][q][1][d] = + -1.0 * data.shape_grads[dof_index1][q][1][d]; + data.shape_grads[dof_index2_1][q][2][d] = + data.shape_grads[dof_index1][q][2][d]; + + data.shape_grads[dof_index2_2][q][0][d] = + data.shape_grads[dof_index1][q][0][d]; + data.shape_grads[dof_index2_2][q][1][d] = + -1.0 * data.shape_grads[dof_index1][q][1][d]; + data.shape_grads[dof_index2_2][q][2][d] = + -1.0 * data.shape_grads[dof_index1][q][2][d]; + } + } + // Type 3: (note we re-use k and j for + // convenience): + const unsigned int shift_ij(j + + shift_j); // here we've subbed + // j for i, k for j. + const unsigned int dof_index3_1(cell_type3_offset1 + + shift_ij); + const unsigned int dof_index3_2(cell_type3_offset2 + + shift_ij); + const unsigned int dof_index3_3(cell_type3_offset3 + + shift_ij); + for (unsigned int d1 = 0; d1 < dim; ++d1) + { + for (unsigned int d2 = 0; d2 < dim; ++d2) + { + data.shape_grads[dof_index3_1][q][d1][d2] = 0.0; + data.shape_grads[dof_index3_2][q][d1][d2] = 0.0; + data.shape_grads[dof_index3_3][q][d1][d2] = 0.0; + } + } + data.shape_grads[dof_index3_1][q][0][1] = + 2.0 * polyy[j][1] * polyz[k][0]; + data.shape_grads[dof_index3_1][q][0][2] = + 2.0 * polyy[j][0] * polyz[k][1]; + + data.shape_grads[dof_index3_2][q][1][0] = + 2.0 * polyx[j][1] * polyz[k][0]; + data.shape_grads[dof_index3_2][q][1][2] = + 2.0 * polyx[j][0] * polyz[k][1]; + + data.shape_grads[dof_index3_3][q][2][0] = + 2.0 * polyx[j][1] * polyy[k][0]; + data.shape_grads[dof_index3_3][q][2][1] = + 2.0 * polyx[j][0] * polyy[k][1]; + } } - - data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0; } - // edge_lambda_grads - for (unsigned int q = 0; q < n_q_points; ++q) + if (flags & update_hessians) { - double x(p_list[q][0]); - double y(p_list[q][1]); - double z(p_list[q][2]); - data.edge_lambda_grads_3d[0][q] = {z - 1.0, 0.0, x - 1.0}; - data.edge_lambda_grads_3d[1][q] = {1.0 - z, 0.0, -x}; - data.edge_lambda_grads_3d[2][q] = {0.0, z - 1.0, y - 1.0}; - data.edge_lambda_grads_3d[3][q] = {0.0, 1.0 - z, -y}; - data.edge_lambda_grads_3d[4][q] = {-z, 0.0, 1.0 - x}; - data.edge_lambda_grads_3d[5][q] = {z, 0.0, x}; - data.edge_lambda_grads_3d[6][q] = {0.0, -z, 1.0 - y}; - data.edge_lambda_grads_3d[7][q] = {0.0, z, y}; - data.edge_lambda_grads_3d[8][q] = {y - 1.0, x - 1.0, 0.0}; - data.edge_lambda_grads_3d[9][q] = {1.0 - y, -x, 0.0}; - data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0}; - data.edge_lambda_grads_3d[11][q] = {y, x, 0.0}; - } - // edge_lambda gradgrads: - const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = { - 1, - -1, - 1, - -1, - -1, - 1, - -1, - 1, - 1, - -1, - -1, - 1}; // sign of the 2nd derivative for each edge. - const unsigned int - edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = { - {0, 2}, - {0, 2}, - {1, 2}, - {1, 2}, - {0, 2}, - {0, 2}, - {1, 2}, - {1, 2}, - {0, 1}, - {0, 1}, - {0, 1}, - {0, 1}}; // component which edge_lambda[m] depends on. - for (unsigned int m = 0; m < lines_per_cell; ++m) - { - data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]] - [edge_lambda_directions[m][1]] = - edge_lambda_sign[m]; - data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]] - [edge_lambda_directions[m][0]] = - edge_lambda_sign[m]; - } - // Precomputation for higher order shape functions, - // and the face parameterisation. - if (degree > 0) - { - // resize required data: - data.face_lambda_values.resize(faces_per_cell); - data.face_lambda_grads.resize(faces_per_cell); - // for face-based shape functions: - for (unsigned int m = 0; m < faces_per_cell; ++m) - { - data.face_lambda_values[m].resize(n_q_points); - data.face_lambda_grads[m].resize(3); - } - // Fill in the values (these don't change between cells). - for (unsigned int q = 0; q < n_q_points; ++q) - { - double x(p_list[q][0]); - double y(p_list[q][1]); - double z(p_list[q][2]); - data.face_lambda_values[0][q] = 1.0 - x; - data.face_lambda_values[1][q] = x; - data.face_lambda_values[2][q] = 1.0 - y; - data.face_lambda_values[3][q] = y; - data.face_lambda_values[4][q] = 1.0 - z; - data.face_lambda_values[5][q] = z; - } - // gradients are constant: - data.face_lambda_grads[0] = {-1.0, 0.0, 0.0}; - data.face_lambda_grads[1] = {1.0, 0.0, 0.0}; - data.face_lambda_grads[2] = {0.0, -1.0, 0.0}; - data.face_lambda_grads[3] = {0.0, 1.0, 0.0}; - data.face_lambda_grads[4] = {0.0, 0.0, -1.0}; - data.face_lambda_grads[5] = {0.0, 0.0, 1.0}; - - // for cell-based shape functions: - // these don't depend on the cell, so can precompute all here: - if (flags & (update_values | update_gradients | update_hessians)) + for (unsigned int k = 0; k < degree; ++k) { - // Cell-based shape functions: - // - // Type-1 (gradients): - // \phi^{C_{1}}_{ijk} = grad( - // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ), - // - // 0 <= i,j,k < degree. (in a group of degree*degree*degree) - const unsigned int cell_type1_offset(n_line_dofs + - n_face_dofs); - // Type-2: - // - // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk} - // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1, - // -1)\phi^{C_{1}}_{ijk} - // - // 0 <= i,j,k < degree. (subtypes in groups of - // degree*degree*degree) - // - // here we order so that all of subtype 1 comes first, then - // subtype 2. - const unsigned int cell_type2_offset1( - cell_type1_offset + degree * degree * degree); - const unsigned int cell_type2_offset2( - cell_type2_offset1 + degree * degree * degree); - // Type-3 - // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x} - // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y} - // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z} - // - // 0 <= i,j,k < degree. (subtypes in groups of degree*degree) - // - // again we order so we compute all of subtype 1 first, then - // subtype 2, etc. - const unsigned int cell_type3_offset1( - cell_type2_offset2 + degree * degree * degree); - const unsigned int cell_type3_offset2(cell_type3_offset1 + - degree * degree); - const unsigned int cell_type3_offset3(cell_type3_offset2 + - degree * degree); - - // compute all points we must evaluate the 1d polynomials at: - std::vector> cell_points(n_q_points); - for (unsigned int q = 0; q < n_q_points; ++q) - { - for (unsigned int d = 0; d < dim; ++d) - { - cell_points[q][d] = 2.0 * p_list[q][d] - 1.0; - } - } + const unsigned int shift_k(k * degree * degree); + const unsigned int shift_j(k * + degree); // Used below when subbing + // k for j type 3 - // We only need poly values and 1st derivative for - // update_values, but need the 2nd derivative too for - // update_gradients. For update_hessians we also need 3rd - // derivative. - const unsigned int poly_length = - (flags & update_hessians) ? - 4 : - ((flags & update_gradients) ? 3 : 2); - - // Loop through quad points: - for (unsigned int q = 0; q < n_q_points; ++q) + for (unsigned int j = 0; j < degree; ++j) { - // pre-compute values & required derivatives at this quad - // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy = - // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1). - // - // for each polyc[d], c=x,y,z, contains the d-th - // derivative with respect to the coordinate c. - std::vector> polyx( - degree, std::vector(poly_length)); - std::vector> polyy( - degree, std::vector(poly_length)); - std::vector> polyz( - degree, std::vector(poly_length)); + const unsigned int shift_jk(j * degree + shift_k); for (unsigned int i = 0; i < degree; ++i) { - // compute all required 1d polynomials for i - IntegratedLegendrePolynomials[i + 2].value( - cell_points[q][0], polyx[i]); - IntegratedLegendrePolynomials[i + 2].value( - cell_points[q][1], polyy[i]); - IntegratedLegendrePolynomials[i + 2].value( - cell_points[q][2], polyz[i]); - } - // Now use these to compute the shape functions: - if (flags & update_values) - { - for (unsigned int k = 0; k < degree; ++k) + const unsigned int shift_ijk(shift_jk + i); + + // Type 1: + const unsigned int dof_index1(cell_type1_offset + + shift_ijk); + + data.shape_hessians[dof_index1][q][0][0][0] = + 8.0 * polyx[i][3] * polyy[j][0] * polyz[k][0]; + data.shape_hessians[dof_index1][q][1][0][0] = + 8.0 * polyx[i][2] * polyy[j][1] * polyz[k][0]; + data.shape_hessians[dof_index1][q][2][0][0] = + 8.0 * polyx[i][2] * polyy[j][0] * polyz[k][1]; + + data.shape_hessians[dof_index1][q][0][1][0] = + data.shape_hessians[dof_index1][q][1][0][0]; + data.shape_hessians[dof_index1][q][1][1][0] = + 8.0 * polyx[i][1] * polyy[j][2] * polyz[k][0]; + data.shape_hessians[dof_index1][q][2][1][0] = + 8.0 * polyx[i][1] * polyy[j][1] * polyz[k][1]; + + data.shape_hessians[dof_index1][q][0][2][0] = + data.shape_hessians[dof_index1][q][2][0][0]; + data.shape_hessians[dof_index1][q][1][2][0] = + data.shape_hessians[dof_index1][q][2][1][0]; + data.shape_hessians[dof_index1][q][2][2][0] = + 8.0 * polyx[i][1] * polyy[j][0] * polyz[k][2]; + + + data.shape_hessians[dof_index1][q][0][0][1] = + data.shape_hessians[dof_index1][q][1][0][0]; + data.shape_hessians[dof_index1][q][1][0][1] = + data.shape_hessians[dof_index1][q][1][1][0]; + data.shape_hessians[dof_index1][q][2][0][1] = + data.shape_hessians[dof_index1][q][2][1][0]; + + data.shape_hessians[dof_index1][q][0][1][1] = + data.shape_hessians[dof_index1][q][1][1][0]; + data.shape_hessians[dof_index1][q][1][1][1] = + 8.0 * polyx[i][0] * polyy[j][3] * polyz[k][0]; + data.shape_hessians[dof_index1][q][2][1][1] = + 8.0 * polyx[i][0] * polyy[j][2] * polyz[k][1]; + + data.shape_hessians[dof_index1][q][0][2][1] = + data.shape_hessians[dof_index1][q][2][1][0]; + data.shape_hessians[dof_index1][q][1][2][1] = + data.shape_hessians[dof_index1][q][2][1][1]; + data.shape_hessians[dof_index1][q][2][2][1] = + 8.0 * polyx[i][0] * polyy[j][1] * polyz[k][2]; + + + data.shape_hessians[dof_index1][q][0][0][2] = + data.shape_hessians[dof_index1][q][2][0][0]; + data.shape_hessians[dof_index1][q][1][0][2] = + data.shape_hessians[dof_index1][q][2][1][0]; + data.shape_hessians[dof_index1][q][2][0][2] = + data.shape_hessians[dof_index1][q][2][2][0]; + + data.shape_hessians[dof_index1][q][0][1][2] = + data.shape_hessians[dof_index1][q][2][1][0]; + data.shape_hessians[dof_index1][q][1][1][2] = + data.shape_hessians[dof_index1][q][2][1][1]; + data.shape_hessians[dof_index1][q][2][1][2] = + data.shape_hessians[dof_index1][q][2][2][1]; + + data.shape_hessians[dof_index1][q][0][2][2] = + data.shape_hessians[dof_index1][q][2][2][0]; + data.shape_hessians[dof_index1][q][1][2][2] = + data.shape_hessians[dof_index1][q][2][2][1]; + data.shape_hessians[dof_index1][q][2][2][2] = + 8.0 * polyx[i][0] * polyy[j][0] * polyz[k][3]; + + + // Type 2: + const unsigned int dof_index2_1(cell_type2_offset1 + + shift_ijk); + const unsigned int dof_index2_2(cell_type2_offset2 + + shift_ijk); + + for (unsigned int d1 = 0; d1 < dim; ++d1) { - const unsigned int shift_k(k * degree * degree); - const unsigned int shift_j( - k * degree); // Used below when subbing k for j - // (type 3) - for (unsigned int j = 0; j < degree; ++j) + for (unsigned int d2 = 0; d2 < dim; ++d2) { - const unsigned int shift_jk(j * degree + - shift_k); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ijk(shift_jk + - i); - - // Type 1: - const unsigned int dof_index1( - cell_type1_offset + shift_ijk); - - data.shape_values[dof_index1][q][0] = - 2.0 * polyx[i][1] * polyy[j][0] * - polyz[k][0]; - data.shape_values[dof_index1][q][1] = - 2.0 * polyx[i][0] * polyy[j][1] * - polyz[k][0]; - data.shape_values[dof_index1][q][2] = - 2.0 * polyx[i][0] * polyy[j][0] * - polyz[k][1]; - - // Type 2: - const unsigned int dof_index2_1( - cell_type2_offset1 + shift_ijk); - const unsigned int dof_index2_2( - cell_type2_offset2 + shift_ijk); - - data.shape_values[dof_index2_1][q][0] = - data.shape_values[dof_index1][q][0]; - data.shape_values[dof_index2_1][q][1] = - -1.0 * - data.shape_values[dof_index1][q][1]; - data.shape_values[dof_index2_1][q][2] = - data.shape_values[dof_index1][q][2]; - - data.shape_values[dof_index2_2][q][0] = - data.shape_values[dof_index1][q][0]; - data.shape_values[dof_index2_2][q][1] = - -1.0 * - data.shape_values[dof_index1][q][1]; - data.shape_values[dof_index2_2][q][2] = - -1.0 * - data.shape_values[dof_index1][q][2]; - } - // Type 3: (note we re-use k and j for - // convenience): - const unsigned int shift_ij( - j + shift_j); // here we've subbed j for i, - // k for j. - const unsigned int dof_index3_1( - cell_type3_offset1 + shift_ij); - const unsigned int dof_index3_2( - cell_type3_offset2 + shift_ij); - const unsigned int dof_index3_3( - cell_type3_offset3 + shift_ij); - - data.shape_values[dof_index3_1][q][0] = - polyy[j][0] * polyz[k][0]; - data.shape_values[dof_index3_1][q][1] = 0.0; - data.shape_values[dof_index3_1][q][2] = 0.0; - - data.shape_values[dof_index3_2][q][0] = 0.0; - data.shape_values[dof_index3_2][q][1] = - polyx[j][0] * polyz[k][0]; - data.shape_values[dof_index3_2][q][2] = 0.0; - - data.shape_values[dof_index3_3][q][0] = 0.0; - data.shape_values[dof_index3_3][q][1] = 0.0; - data.shape_values[dof_index3_3][q][2] = - polyx[j][0] * polyy[k][0]; + data.shape_hessians[dof_index2_1][q][0][d1] + [d2] = + data + .shape_hessians[dof_index1][q][0][d1][d2]; + data.shape_hessians[dof_index2_1][q][1][d1] + [d2] = + -1.0 * + data + .shape_hessians[dof_index1][q][1][d1][d2]; + data.shape_hessians[dof_index2_1][q][2][d1] + [d2] = + data + .shape_hessians[dof_index1][q][2][d1][d2]; + + data.shape_hessians[dof_index2_2][q][0][d1] + [d2] = + data + .shape_hessians[dof_index1][q][0][d1][d2]; + data.shape_hessians[dof_index2_2][q][1][d1] + [d2] = + -1.0 * + data + .shape_hessians[dof_index1][q][1][d1][d2]; + data.shape_hessians[dof_index2_2][q][2][d1] + [d2] = + -1.0 * + data + .shape_hessians[dof_index1][q][2][d1][d2]; } } } - if (flags & update_gradients) + // Type 3: (note we re-use k and j for + // convenience): + const unsigned int shift_ij(j + + shift_j); // here we've subbed + // j for i, k for j. + const unsigned int dof_index3_1(cell_type3_offset1 + + shift_ij); + const unsigned int dof_index3_2(cell_type3_offset2 + + shift_ij); + const unsigned int dof_index3_3(cell_type3_offset3 + + shift_ij); + for (unsigned int d1 = 0; d1 < dim; ++d1) { - for (unsigned int k = 0; k < degree; ++k) + for (unsigned int d2 = 0; d2 < dim; ++d2) { - const unsigned int shift_k(k * degree * degree); - const unsigned int shift_j( - k * degree); // Used below when subbing k for j - // (type 3) - for (unsigned int j = 0; j < degree; ++j) + for (unsigned int d3 = 0; d3 < dim; ++d3) { - const unsigned int shift_jk(j * degree + - shift_k); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ijk(shift_jk + - i); - - // Type 1: - const unsigned int dof_index1( - cell_type1_offset + shift_ijk); - - data.shape_grads[dof_index1][q][0][0] = - 4.0 * polyx[i][2] * polyy[j][0] * - polyz[k][0]; - data.shape_grads[dof_index1][q][0][1] = - 4.0 * polyx[i][1] * polyy[j][1] * - polyz[k][0]; - data.shape_grads[dof_index1][q][0][2] = - 4.0 * polyx[i][1] * polyy[j][0] * - polyz[k][1]; - - data.shape_grads[dof_index1][q][1][0] = - data.shape_grads[dof_index1][q][0][1]; - data.shape_grads[dof_index1][q][1][1] = - 4.0 * polyx[i][0] * polyy[j][2] * - polyz[k][0]; - data.shape_grads[dof_index1][q][1][2] = - 4.0 * polyx[i][0] * polyy[j][1] * - polyz[k][1]; - - data.shape_grads[dof_index1][q][2][0] = - data.shape_grads[dof_index1][q][0][2]; - data.shape_grads[dof_index1][q][2][1] = - data.shape_grads[dof_index1][q][1][2]; - data.shape_grads[dof_index1][q][2][2] = - 4.0 * polyx[i][0] * polyy[j][0] * - polyz[k][2]; - - // Type 2: - const unsigned int dof_index2_1( - cell_type2_offset1 + shift_ijk); - const unsigned int dof_index2_2( - cell_type2_offset2 + shift_ijk); - - for (unsigned int d = 0; d < dim; ++d) - { - data.shape_grads[dof_index2_1][q][0] - [d] = - data - .shape_grads[dof_index1][q][0][d]; - data.shape_grads[dof_index2_1][q][1] - [d] = - -1.0 * - data - .shape_grads[dof_index1][q][1][d]; - data.shape_grads[dof_index2_1][q][2] - [d] = - data - .shape_grads[dof_index1][q][2][d]; - - data.shape_grads[dof_index2_2][q][0] - [d] = - data - .shape_grads[dof_index1][q][0][d]; - data.shape_grads[dof_index2_2][q][1] - [d] = - -1.0 * - data - .shape_grads[dof_index1][q][1][d]; - data.shape_grads[dof_index2_2][q][2] - [d] = - -1.0 * - data - .shape_grads[dof_index1][q][2][d]; - } - } - // Type 3: (note we re-use k and j for - // convenience): - const unsigned int shift_ij( - j + shift_j); // here we've subbed j for i, - // k for j. - const unsigned int dof_index3_1( - cell_type3_offset1 + shift_ij); - const unsigned int dof_index3_2( - cell_type3_offset2 + shift_ij); - const unsigned int dof_index3_3( - cell_type3_offset3 + shift_ij); - for (unsigned int d1 = 0; d1 < dim; ++d1) - { - for (unsigned int d2 = 0; d2 < dim; ++d2) - { - data.shape_grads[dof_index3_1][q][d1] - [d2] = 0.0; - data.shape_grads[dof_index3_2][q][d1] - [d2] = 0.0; - data.shape_grads[dof_index3_3][q][d1] - [d2] = 0.0; - } - } - data.shape_grads[dof_index3_1][q][0][1] = - 2.0 * polyy[j][1] * polyz[k][0]; - data.shape_grads[dof_index3_1][q][0][2] = - 2.0 * polyy[j][0] * polyz[k][1]; - - data.shape_grads[dof_index3_2][q][1][0] = - 2.0 * polyx[j][1] * polyz[k][0]; - data.shape_grads[dof_index3_2][q][1][2] = - 2.0 * polyx[j][0] * polyz[k][1]; - - data.shape_grads[dof_index3_3][q][2][0] = - 2.0 * polyx[j][1] * polyy[k][0]; - data.shape_grads[dof_index3_3][q][2][1] = - 2.0 * polyx[j][0] * polyy[k][1]; + data.shape_hessians[dof_index3_1][q][d1][d2] + [d3] = 0.0; + data.shape_hessians[dof_index3_2][q][d1][d2] + [d3] = 0.0; + data.shape_hessians[dof_index3_3][q][d1][d2] + [d3] = 0.0; } } } - if (flags & update_hessians) - { - for (unsigned int k = 0; k < degree; ++k) - { - const unsigned int shift_k(k * degree * degree); - const unsigned int shift_j( - k * degree); // Used below when subbing k for j - // type 3 + data.shape_hessians[dof_index3_1][q][0][1][1] = + 4.0 * polyy[j][2] * polyz[k][0]; + data.shape_hessians[dof_index3_1][q][0][1][2] = + 4.0 * polyy[j][1] * polyz[k][1]; - for (unsigned int j = 0; j < degree; ++j) - { - const unsigned int shift_jk(j * degree + - shift_k); - for (unsigned int i = 0; i < degree; ++i) - { - const unsigned int shift_ijk(shift_jk + - i); - - // Type 1: - const unsigned int dof_index1( - cell_type1_offset + shift_ijk); - - data.shape_hessians[dof_index1][q][0][0] - [0] = - 8.0 * polyx[i][3] * polyy[j][0] * - polyz[k][0]; - data.shape_hessians[dof_index1][q][1][0] - [0] = - 8.0 * polyx[i][2] * polyy[j][1] * - polyz[k][0]; - data.shape_hessians[dof_index1][q][2][0] - [0] = - 8.0 * polyx[i][2] * polyy[j][0] * - polyz[k][1]; - - data.shape_hessians[dof_index1][q][0][1] - [0] = - data.shape_hessians[dof_index1][q][1][0] - [0]; - data.shape_hessians[dof_index1][q][1][1] - [0] = - 8.0 * polyx[i][1] * polyy[j][2] * - polyz[k][0]; - data.shape_hessians[dof_index1][q][2][1] - [0] = - 8.0 * polyx[i][1] * polyy[j][1] * - polyz[k][1]; - - data.shape_hessians[dof_index1][q][0][2] - [0] = - data.shape_hessians[dof_index1][q][2][0] - [0]; - data.shape_hessians[dof_index1][q][1][2] - [0] = - data.shape_hessians[dof_index1][q][2][1] - [0]; - data.shape_hessians[dof_index1][q][2][2] - [0] = - 8.0 * polyx[i][1] * polyy[j][0] * - polyz[k][2]; - - - data.shape_hessians[dof_index1][q][0][0] - [1] = - data.shape_hessians[dof_index1][q][1][0] - [0]; - data.shape_hessians[dof_index1][q][1][0] - [1] = - data.shape_hessians[dof_index1][q][1][1] - [0]; - data.shape_hessians[dof_index1][q][2][0] - [1] = - data.shape_hessians[dof_index1][q][2][1] - [0]; - - data.shape_hessians[dof_index1][q][0][1] - [1] = - data.shape_hessians[dof_index1][q][1][1] - [0]; - data.shape_hessians[dof_index1][q][1][1] - [1] = - 8.0 * polyx[i][0] * polyy[j][3] * - polyz[k][0]; - data.shape_hessians[dof_index1][q][2][1] - [1] = - 8.0 * polyx[i][0] * polyy[j][2] * - polyz[k][1]; - - data.shape_hessians[dof_index1][q][0][2] - [1] = - data.shape_hessians[dof_index1][q][2][1] - [0]; - data.shape_hessians[dof_index1][q][1][2] - [1] = - data.shape_hessians[dof_index1][q][2][1] - [1]; - data.shape_hessians[dof_index1][q][2][2] - [1] = - 8.0 * polyx[i][0] * polyy[j][1] * - polyz[k][2]; - - - data.shape_hessians[dof_index1][q][0][0] - [2] = - data.shape_hessians[dof_index1][q][2][0] - [0]; - data.shape_hessians[dof_index1][q][1][0] - [2] = - data.shape_hessians[dof_index1][q][2][1] - [0]; - data.shape_hessians[dof_index1][q][2][0] - [2] = - data.shape_hessians[dof_index1][q][2][2] - [0]; - - data.shape_hessians[dof_index1][q][0][1] - [2] = - data.shape_hessians[dof_index1][q][2][1] - [0]; - data.shape_hessians[dof_index1][q][1][1] - [2] = - data.shape_hessians[dof_index1][q][2][1] - [1]; - data.shape_hessians[dof_index1][q][2][1] - [2] = - data.shape_hessians[dof_index1][q][2][2] - [1]; - - data.shape_hessians[dof_index1][q][0][2] - [2] = - data.shape_hessians[dof_index1][q][2][2] - [0]; - data.shape_hessians[dof_index1][q][1][2] - [2] = - data.shape_hessians[dof_index1][q][2][2] - [1]; - data.shape_hessians[dof_index1][q][2][2] - [2] = - 8.0 * polyx[i][0] * polyy[j][0] * - polyz[k][3]; - - - // Type 2: - const unsigned int dof_index2_1( - cell_type2_offset1 + shift_ijk); - const unsigned int dof_index2_2( - cell_type2_offset2 + shift_ijk); - - for (unsigned int d1 = 0; d1 < dim; ++d1) - { - for (unsigned int d2 = 0; d2 < dim; - ++d2) - { - data - .shape_hessians[dof_index2_1][q] - [0][d1][d2] = - data - .shape_hessians[dof_index1][q] - [0][d1][d2]; - data - .shape_hessians[dof_index2_1][q] - [1][d1][d2] = - -1.0 * - data - .shape_hessians[dof_index1][q] - [1][d1][d2]; - data - .shape_hessians[dof_index2_1][q] - [2][d1][d2] = - data - .shape_hessians[dof_index1][q] - [2][d1][d2]; - - data - .shape_hessians[dof_index2_2][q] - [0][d1][d2] = - data - .shape_hessians[dof_index1][q] - [0][d1][d2]; - data - .shape_hessians[dof_index2_2][q] - [1][d1][d2] = - -1.0 * - data - .shape_hessians[dof_index1][q] - [1][d1][d2]; - data - .shape_hessians[dof_index2_2][q] - [2][d1][d2] = - -1.0 * - data - .shape_hessians[dof_index1][q] - [2][d1][d2]; - } - } - } - // Type 3: (note we re-use k and j for - // convenience): - const unsigned int shift_ij( - j + shift_j); // here we've subbed j for i, - // k for j. - const unsigned int dof_index3_1( - cell_type3_offset1 + shift_ij); - const unsigned int dof_index3_2( - cell_type3_offset2 + shift_ij); - const unsigned int dof_index3_3( - cell_type3_offset3 + shift_ij); - for (unsigned int d1 = 0; d1 < dim; ++d1) - { - for (unsigned int d2 = 0; d2 < dim; ++d2) - { - for (unsigned int d3 = 0; d3 < dim; - ++d3) - { - data - .shape_hessians[dof_index3_1][q] - [d1][d2][d3] = - 0.0; - data - .shape_hessians[dof_index3_2][q] - [d1][d2][d3] = - 0.0; - data - .shape_hessians[dof_index3_3][q] - [d1][d2][d3] = - 0.0; - } - } - } - data - .shape_hessians[dof_index3_1][q][0][1][1] = - 4.0 * polyy[j][2] * polyz[k][0]; - data - .shape_hessians[dof_index3_1][q][0][1][2] = - 4.0 * polyy[j][1] * polyz[k][1]; - - data - .shape_hessians[dof_index3_1][q][0][2][1] = - data - .shape_hessians[dof_index3_1][q][0][1][2]; - data - .shape_hessians[dof_index3_1][q][0][2][2] = - 4.0 * polyy[j][0] * polyz[k][2]; + data.shape_hessians[dof_index3_1][q][0][2][1] = + data.shape_hessians[dof_index3_1][q][0][1][2]; + data.shape_hessians[dof_index3_1][q][0][2][2] = + 4.0 * polyy[j][0] * polyz[k][2]; - data - .shape_hessians[dof_index3_2][q][1][0][0] = - 4.0 * polyx[j][2] * polyz[k][0]; - data - .shape_hessians[dof_index3_2][q][1][0][2] = - 4.0 * polyx[j][1] * polyz[k][1]; + data.shape_hessians[dof_index3_2][q][1][0][0] = + 4.0 * polyx[j][2] * polyz[k][0]; + data.shape_hessians[dof_index3_2][q][1][0][2] = + 4.0 * polyx[j][1] * polyz[k][1]; - data - .shape_hessians[dof_index3_2][q][1][2][0] = - data - .shape_hessians[dof_index3_2][q][1][0][2]; - data - .shape_hessians[dof_index3_2][q][1][2][2] = - 4.0 * polyx[j][0] * polyz[k][2]; + data.shape_hessians[dof_index3_2][q][1][2][0] = + data.shape_hessians[dof_index3_2][q][1][0][2]; + data.shape_hessians[dof_index3_2][q][1][2][2] = + 4.0 * polyx[j][0] * polyz[k][2]; - data - .shape_hessians[dof_index3_3][q][2][0][0] = - 4.0 * polyx[j][2] * polyy[k][0]; - data - .shape_hessians[dof_index3_3][q][2][0][1] = - 4.0 * polyx[j][1] * polyy[k][1]; + data.shape_hessians[dof_index3_3][q][2][0][0] = + 4.0 * polyx[j][2] * polyy[k][0]; + data.shape_hessians[dof_index3_3][q][2][0][1] = + 4.0 * polyx[j][1] * polyy[k][1]; - data - .shape_hessians[dof_index3_3][q][2][1][0] = - data - .shape_hessians[dof_index3_3][q][2][0][1]; - data - .shape_hessians[dof_index3_3][q][2][1][1] = - 4.0 * polyx[j][0] * polyy[k][2]; - } - } - } + data.shape_hessians[dof_index3_3][q][2][1][0] = + data.shape_hessians[dof_index3_3][q][2][0][1]; + data.shape_hessians[dof_index3_3][q][2][1][1] = + 4.0 * polyx[j][0] * polyy[k][2]; } } } - break; - } - default: - { - DEAL_II_NOT_IMPLEMENTED(); } } +} + + + +template +std::unique_ptr::InternalDataBase> +FE_NedelecSZ::get_data( + const UpdateFlags update_flags, + const Mapping & /*mapping*/, + const Quadrature &quadrature, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + & /*output_data*/) const +{ + std::unique_ptr< + typename dealii::FiniteElement::InternalDataBase> + data_ptr = std::make_unique(); + + const unsigned int n_q_points = quadrature.size(); + + std::vector> p_list(n_q_points); + p_list = quadrature.get_points(); + + this->evaluate(p_list, update_flags, data_ptr); + return data_ptr; } @@ -3415,8 +3272,6 @@ FE_NedelecSZ::create_polynomials(const unsigned int degree) IntegratedLegendreSZ::generate_complete_basis(degree + 1); } - - // explicit instantiations #include "fe_nedelec_sz.inst"