From: David Wells Date: Mon, 2 Jun 2025 03:04:45 +0000 (-0400) Subject: FE_PolyTensor: Combine the filling functions. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=83be794e817f1893e27ca5146d189a667fa1a2f0;p=dealii.git FE_PolyTensor: Combine the filling functions. These three functions are essentially identical. Furthermore, given the difficulty with orientations (a cell may be a translation of another cell and have faces in a nonstandard orientation) completely ignore the CellSimilarity subsystem. --- diff --git a/include/deal.II/fe/fe_poly_tensor.h b/include/deal.II/fe/fe_poly_tensor.h index 8ad405734f..b31004fe73 100644 --- a/include/deal.II/fe/fe_poly_tensor.h +++ b/include/deal.II/fe/fe_poly_tensor.h @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -207,6 +208,27 @@ public: const unsigned int component) const override; protected: +#ifndef DOXYGEN + class InternalData; +#endif + + /** + * Internal function called by fill_fe_values(), fill_fe_face_values(), + * and fill_fe_subface_values(). + */ + void + compute_fill( + const typename Triangulation::cell_iterator &cell, + const typename QProjector::DataSetDescriptor &offset, + const unsigned int n_q_points, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FE_PolyTensor::InternalData &fe_internal, + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const; + /** * The mapping type to be used to map shape functions from the reference * cell to the mesh cell. If this vector is length one, the same mapping diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc index bd2bc5d866..31f70ffb68 100644 --- a/source/fe/fe_poly_tensor.cc +++ b/source/fe/fe_poly_tensor.cc @@ -527,1308 +527,65 @@ template void FE_PolyTensor::fill_fe_values( const typename Triangulation::cell_iterator &cell, - const CellSimilarity::Similarity cell_similarity, - const Quadrature &quadrature, - const Mapping &mapping, - const typename Mapping::InternalDataBase &mapping_internal, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const typename FiniteElement::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - &output_data) const -{ - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - - const unsigned int n_q_points = quadrature.size(); - - Assert(!(fe_data.update_each & update_values) || - fe_data.shape_values.size()[0] == this->n_dofs_per_cell(), - ExcDimensionMismatch(fe_data.shape_values.size()[0], - this->n_dofs_per_cell())); - Assert(!(fe_data.update_each & update_values) || - fe_data.shape_values.size()[1] == n_q_points, - ExcDimensionMismatch(fe_data.shape_values.size()[1], n_q_points)); - - // TODO: The dof_sign_change only affects Nedelec elements and is not the - // correct thing on complicated meshes for higher order Nedelec elements. - // Something similar to FE_Q should be done to permute dofs and to change the - // dof signs. A static way using tables (as done in the RaviartThomas - // class) is preferable. - std::fill(fe_data.dof_sign_change.begin(), - fe_data.dof_sign_change.end(), - 1.0); - if (fe_data.update_each & update_values) - internal::FE_PolyTensor::get_dof_sign_change_nedelec( - cell, *this, this->mapping_kind, fe_data.dof_sign_change); - - // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d - // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d. - // Also nothing in 3d since we take care of it by using the - // adjust_quad_dof_sign_for_face_orientation_table. - if (fe_data.update_each & update_values) - internal::FE_PolyTensor::get_dof_sign_change_h_div(cell, - *this, - this->mapping_kind, - fe_data.dof_sign_change); - - // What is the first dof_index on a quad? - const unsigned int first_quad_index = this->get_first_quad_index(); - // How many dofs per quad and how many quad dofs do we have at all? - const unsigned int n_dofs_per_quad = this->n_dofs_per_quad(); - const unsigned int n_quad_dofs = - n_dofs_per_quad * GeometryInfo::faces_per_cell; - - for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell(); - ++dof_index) - { - /* - * This assumes that the dofs are ordered by first vertices, lines, quads - * and volume dofs. Note that in 2d this always gives false. - */ - const bool is_quad_dof = - (dim == 2 ? false : - (first_quad_index <= dof_index) && - (dof_index < first_quad_index + n_quad_dofs)); - - // TODO: This hack is not pretty and it is only here to handle the 2d - // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never - // handled by the - // >>if(is_quad_dof){...}<< but still a possible dof sign change must be - // handled, also for line_dofs in 3d such as in Nedelec. In these cases - // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d - // it is handles with a table. This array is allocated in - // fe_poly_tensor.h. - double dof_sign = 1.0; - // under some circumstances fe_data.dof_sign_change is not allocated - if (fe_data.update_each & update_values) - dof_sign = fe_data.dof_sign_change[dof_index]; - - if (is_quad_dof) - { - /* - * Find the face belonging to this dof_index. This is integer - * division. - */ - const unsigned int face_index_from_dof_index = - (dof_index - first_quad_index) / (n_dofs_per_quad); - - const unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad; - - // Correct the dof_sign if necessary - if (adjust_quad_dof_sign_for_face_orientation( - local_quad_dof_index, - face_index_from_dof_index, - cell->combined_face_orientation(face_index_from_dof_index))) - dof_sign = -1.0; - } - - const MappingKind mapping_kind = get_mapping_kind(dof_index); - - const unsigned int first = - output_data.shape_function_to_row_table - [dof_index * this->n_components() + - this->get_nonzero_components(dof_index).first_selected_component()]; - - // update the shape function values as necessary - // - // we only need to do this if the current cell is not a translation of - // the previous one; or, even if it is a translation, if we use - // mappings other than the standard mappings that require us to - // recompute values and derivatives because of possible sign changes - if (fe_data.update_each & update_values && - ((cell_similarity != CellSimilarity::translation) || - ((mapping_kind == mapping_piola) || - (mapping_kind == mapping_raviart_thomas) || - (mapping_kind == mapping_nedelec)))) - { - switch (mapping_kind) - { - case mapping_none: - { - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - fe_data.shape_values[dof_index][k][d]; - break; - } - - case mapping_covariant: - case mapping_contravariant: - { - mapping.transform( - make_array_view(fe_data.shape_values, dof_index), - mapping_kind, - mapping_internal, - make_array_view(fe_data.transformed_shape_values)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - fe_data.transformed_shape_values[k][d]; - - break; - } - - case mapping_raviart_thomas: - case mapping_piola: - { - mapping.transform( - make_array_view(fe_data.shape_values, dof_index), - mapping_piola, - mapping_internal, - make_array_view(fe_data.transformed_shape_values)); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - dof_sign * fe_data.transformed_shape_values[k][d]; - break; - } - - case mapping_nedelec: - { - mapping.transform( - make_array_view(fe_data.shape_values, dof_index), - mapping_covariant, - mapping_internal, - make_array_view(fe_data.transformed_shape_values)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - dof_sign * fe_data.transformed_shape_values[k][d]; - - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - // update gradients. apply the same logic as above - if (fe_data.update_each & update_gradients && - ((cell_similarity != CellSimilarity::translation) || - ((mapping_kind == mapping_piola) || - (mapping_kind == mapping_raviart_thomas) || - (mapping_kind == mapping_nedelec)))) - - { - switch (mapping_kind) - { - case mapping_none: - { - mapping.transform( - make_array_view(fe_data.shape_grads, dof_index), - mapping_covariant, - mapping_internal, - make_array_view(fe_data.transformed_shape_grads)); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - fe_data.transformed_shape_grads[k][d]; - break; - } - case mapping_covariant: - { - mapping.transform( - make_array_view(fe_data.shape_grads, dof_index), - mapping_covariant_gradient, - mapping_internal, - make_array_view(fe_data.transformed_shape_grads)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - fe_data.transformed_shape_grads[k][d] -= - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][d]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - fe_data.transformed_shape_grads[k][d]; - - break; - } - case mapping_contravariant: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k] = - fe_data.shape_grads[dof_index][k]; - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads), - mapping_contravariant_gradient, - mapping_internal, - make_array_view(fe_data.transformed_shape_grads)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - fe_data.transformed_shape_grads[k][d] += - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][d][n]; - - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - fe_data.transformed_shape_grads[k][d]; - - break; - } - case mapping_raviart_thomas: - case mapping_piola: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k] = - fe_data.shape_grads[dof_index][k]; - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads), - mapping_piola_gradient, - mapping_internal, - make_array_view(fe_data.transformed_shape_grads)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - fe_data.transformed_shape_grads[k][d] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_grads[k][d][n]) - - (output_data.shape_values(first + d, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][n]); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - dof_sign * fe_data.transformed_shape_grads[k][d]; - - break; - } - - case mapping_nedelec: - { - // treat the gradients of - // this particular shape - // function at all - // q-points. if Dv is the - // gradient of the shape - // function on the unit - // cell, then - // (J^-T)Dv(J^-1) is the - // value we want to have on - // the real cell. - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k] = - fe_data.shape_grads[dof_index][k]; - - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads), - mapping_covariant_gradient, - mapping_internal, - make_array_view(fe_data.transformed_shape_grads)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - fe_data.transformed_shape_grads[k][d] -= - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][d]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - dof_sign * fe_data.transformed_shape_grads[k][d]; - - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - // update hessians. apply the same logic as above - if (fe_data.update_each & update_hessians && - ((cell_similarity != CellSimilarity::translation) || - ((mapping_kind == mapping_piola) || - (mapping_kind == mapping_raviart_thomas) || - (mapping_kind == mapping_nedelec)))) - - { - switch (mapping_kind) - { - case mapping_none: - { - mapping.transform( - make_array_view(fe_data.shape_grad_grads, dof_index), - mapping_covariant_gradient, - mapping_internal, - make_array_view(fe_data.transformed_shape_hessians)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - fe_data.transformed_shape_hessians[k][d] -= - output_data.shape_gradients[first + d][k][n] * - mapping_data.jacobian_pushed_forward_grads[k][n]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - fe_data.transformed_shape_hessians[k][d]; - - break; - } - case mapping_covariant: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k] = - fe_data.shape_grad_grads[dof_index][k]; - - mapping.transform( - make_array_view( - fe_data.untransformed_shape_hessian_tensors), - mapping_covariant_hessian, - mapping_internal, - make_array_view(fe_data.transformed_shape_hessians)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - fe_data.transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][d][i][j]) + - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][d][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][d]); - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - fe_data.transformed_shape_hessians[k][d]; - - break; - } - case mapping_contravariant: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k] = - fe_data.shape_grad_grads[dof_index][k]; - - mapping.transform( - make_array_view( - fe_data.untransformed_shape_hessian_tensors), - mapping_contravariant_hessian, - mapping_internal, - make_array_view(fe_data.transformed_shape_hessians)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - fe_data.transformed_shape_hessians[k][d][i][j] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][d][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][d][n][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][d][i][n]) - - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]); - for (unsigned int m = 0; m < spacedim; ++m) - fe_data - .transformed_shape_hessians[k][d][i][j] -= - (mapping_data - .jacobian_pushed_forward_grads[k][d][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + n, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][d][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + n, k)); - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - fe_data.transformed_shape_hessians[k][d]; - - break; - } - case mapping_raviart_thomas: - case mapping_piola: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k] = - fe_data.shape_grad_grads[dof_index][k]; - - mapping.transform( - make_array_view( - fe_data.untransformed_shape_hessian_tensors), - mapping_piola_hessian, - mapping_internal, - make_array_view(fe_data.transformed_shape_hessians)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - fe_data.transformed_shape_hessians[k][d][i][j] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][d][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][d][n][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][d][i][n]) - - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]); - - fe_data.transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + d, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][n][i][j]) + - (output_data.shape_gradients[first + d][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][n][j]) + - (output_data.shape_gradients[first + d][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][n][i]); - - for (unsigned int m = 0; m < spacedim; ++m) - { - fe_data - .transformed_shape_hessians[k][d][i][j] -= - (mapping_data - .jacobian_pushed_forward_grads[k][d][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + n, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][d][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + n, k)); - - fe_data - .transformed_shape_hessians[k][d][i][j] += - (mapping_data - .jacobian_pushed_forward_grads[k][n][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + d, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][n][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + d, k)); - } - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - dof_sign * fe_data.transformed_shape_hessians[k][d]; - - break; - } - - case mapping_nedelec: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k] = - fe_data.shape_grad_grads[dof_index][k]; - - mapping.transform( - make_array_view( - fe_data.untransformed_shape_hessian_tensors), - mapping_covariant_hessian, - mapping_internal, - make_array_view(fe_data.transformed_shape_hessians)); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - fe_data.transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][d][i][j]) + - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][d][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][d]); - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - dof_sign * fe_data.transformed_shape_hessians[k][d]; - - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - // third derivatives are not implemented - if (fe_data.update_each & update_3rd_derivatives && - ((cell_similarity != CellSimilarity::translation) || - ((mapping_kind == mapping_piola) || - (mapping_kind == mapping_raviart_thomas) || - (mapping_kind == mapping_nedelec)))) - { - DEAL_II_NOT_IMPLEMENTED(); - } - } -} - - - -template -void -FE_PolyTensor::fill_fe_face_values( - const typename Triangulation::cell_iterator &cell, - const unsigned int face_no, - const hp::QCollection &quadrature, - const Mapping &mapping, - const typename Mapping::InternalDataBase &mapping_internal, - const internal::FEValuesImplementation::MappingRelatedData - &mapping_data, - const typename FiniteElement::InternalDataBase &fe_internal, - dealii::internal::FEValuesImplementation::FiniteElementRelatedData - &output_data) const -{ - AssertDimension(quadrature.size(), 1); - - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible - Assert(dynamic_cast(&fe_internal) != nullptr, - ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - - const unsigned int n_q_points = quadrature[0].size(); - // offset determines which data set - // to take (all data sets for all - // faces are stored contiguously) - - // TODO: The same 'legacy' comments for 2d apply here as well: these classes - // do not handle non-standard orientations in 2d in a way consistent with the - // rest of the library, but are consistent with themselves (see, e.g., the - // fe_conformity_dim_2 tests). - // - // In this case: all of this code was written assuming that QProjector assumed - // that all faces were in the default orientation in 2d, but contains special - // workarounds in case that isn't the case. Hence, to keep those workarounds - // working, we still assume that all faces are in the default orientation. - const auto offset = QProjector::DataSetDescriptor::face( - this->reference_cell(), - face_no, - dim == 2 ? numbers::default_geometric_orientation : - cell->combined_face_orientation(face_no), - n_q_points); - - // TODO: Size assertions - - // TODO: The dof_sign_change only affects Nedelec elements and is not the - // correct thing on complicated meshes for higher order Nedelec elements. - // Something similar to FE_Q should be done to permute dofs and to change the - // dof signs. A static way using tables (as done in the RaviartThomas - // class) is preferable. - std::fill(fe_data.dof_sign_change.begin(), - fe_data.dof_sign_change.end(), - 1.0); - if (fe_data.update_each & update_values) - internal::FE_PolyTensor::get_dof_sign_change_nedelec( - cell, *this, this->mapping_kind, fe_data.dof_sign_change); - - // TODO: This, similarly to the Nedelec case, is just a legacy function in 2d - // and affects only face_dofs of H(div) conformal FEs. It does nothing in 1d. - // Also nothing in 3d since we take care of it by using the - // adjust_quad_dof_sign_for_face_orientation_table. - if (fe_data.update_each & update_values) - internal::FE_PolyTensor::get_dof_sign_change_h_div(cell, - *this, - this->mapping_kind, - fe_data.dof_sign_change); - - // What is the first dof_index on a quad? - const unsigned int first_quad_index = this->get_first_quad_index(); - // How many dofs per quad and how many quad dofs do we have at all? - const unsigned int n_dofs_per_quad = this->n_dofs_per_quad(); - const unsigned int n_quad_dofs = - n_dofs_per_quad * GeometryInfo::faces_per_cell; - - for (unsigned int dof_index = 0; dof_index < this->n_dofs_per_cell(); - ++dof_index) - { - /* - * This assumes that the dofs are ordered by first vertices, lines, quads - * and volume dofs. Note that in 2d this always gives false. - */ - const bool is_quad_dof = - (dim == 2 ? false : - (first_quad_index <= dof_index) && - (dof_index < first_quad_index + n_quad_dofs)); - - // TODO: This hack is not pretty and it is only here to handle the 2d - // case and the Nedelec legacy case. In 2d dof_sign of a face_dof is never - // handled by the - // >>if(is_quad_dof){...}<< but still a possible dof sign change must be - // handled, also for line_dofs in 3d such as in Nedelec. In these cases - // this is encoded in the array fe_data.dof_sign_change[dof_index]. In 3d - // it is handles with a table. This array is allocated in - // fe_poly_tensor.h. - double dof_sign = 1.0; - // under some circumstances fe_data.dof_sign_change is not allocated - if (fe_data.update_each & update_values) - dof_sign = fe_data.dof_sign_change[dof_index]; - - if (is_quad_dof) - { - /* - * Find the face belonging to this dof_index. This is integer - * division. - */ - unsigned int face_index_from_dof_index = - (dof_index - first_quad_index) / (n_dofs_per_quad); - - unsigned int local_quad_dof_index = dof_index % n_dofs_per_quad; - - // Correct the dof_sign if necessary - if (adjust_quad_dof_sign_for_face_orientation( - local_quad_dof_index, - face_index_from_dof_index, - cell->combined_face_orientation(face_index_from_dof_index))) - dof_sign = -1.0; - } - - const MappingKind mapping_kind = get_mapping_kind(dof_index); - - const unsigned int first = - output_data.shape_function_to_row_table - [dof_index * this->n_components() + - this->get_nonzero_components(dof_index).first_selected_component()]; - - if (fe_data.update_each & update_values) - { - switch (mapping_kind) - { - case mapping_none: - { - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - fe_data.shape_values[dof_index][k + offset][d]; - break; - } - - case mapping_covariant: - case mapping_contravariant: - { - const ArrayView> - transformed_shape_values = - make_array_view(fe_data.transformed_shape_values, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_values, - dof_index, - offset, - n_q_points), - mapping_kind, - mapping_internal, - transformed_shape_values); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - transformed_shape_values[k][d]; - - break; - } - case mapping_raviart_thomas: - case mapping_piola: - { - const ArrayView> - transformed_shape_values = - make_array_view(fe_data.transformed_shape_values, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_values, - dof_index, - offset, - n_q_points), - mapping_piola, - mapping_internal, - transformed_shape_values); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - dof_sign * transformed_shape_values[k][d]; - break; - } - - case mapping_nedelec: - { - const ArrayView> - transformed_shape_values = - make_array_view(fe_data.transformed_shape_values, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_values, - dof_index, - offset, - n_q_points), - mapping_covariant, - mapping_internal, - transformed_shape_values); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_values(first + d, k) = - dof_sign * transformed_shape_values[k][d]; - - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - if (fe_data.update_each & update_gradients) - { - switch (mapping_kind) - { - case mapping_none: - { - const ArrayView> transformed_shape_grads = - make_array_view(fe_data.transformed_shape_grads, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_grads, - dof_index, - offset, - n_q_points), - mapping_covariant, - mapping_internal, - transformed_shape_grads); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - transformed_shape_grads[k][d]; - break; - } - - case mapping_covariant: - { - const ArrayView> transformed_shape_grads = - make_array_view(fe_data.transformed_shape_grads, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_grads, - dof_index, - offset, - n_q_points), - mapping_covariant_gradient, - mapping_internal, - transformed_shape_grads); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - transformed_shape_grads[k][d] -= - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][d]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - transformed_shape_grads[k][d]; - break; - } - - case mapping_contravariant: - { - const ArrayView> transformed_shape_grads = - make_array_view(fe_data.transformed_shape_grads, - offset, - n_q_points); - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k + offset] = - fe_data.shape_grads[dof_index][k + offset]; - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads, - offset, - n_q_points), - mapping_contravariant_gradient, - mapping_internal, - transformed_shape_grads); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - transformed_shape_grads[k][d] += - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][d][n]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - transformed_shape_grads[k][d]; - - break; - } - - case mapping_raviart_thomas: - case mapping_piola: - { - const ArrayView> transformed_shape_grads = - make_array_view(fe_data.transformed_shape_grads, - offset, - n_q_points); - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k + offset] = - fe_data.shape_grads[dof_index][k + offset]; - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads, - offset, - n_q_points), - mapping_piola_gradient, - mapping_internal, - transformed_shape_grads); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - transformed_shape_grads[k][d] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_grads[k][d][n]) - - (output_data.shape_values(first + d, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][n]); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - dof_sign * transformed_shape_grads[k][d]; - - break; - } - - case mapping_nedelec: - { - // treat the gradients of - // this particular shape - // function at all - // q-points. if Dv is the - // gradient of the shape - // function on the unit - // cell, then - // (J^-T)Dv(J^-1) is the - // value we want to have on - // the real cell. - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_grads[k + offset] = - fe_data.shape_grads[dof_index][k + offset]; - - const ArrayView> transformed_shape_grads = - make_array_view(fe_data.transformed_shape_grads, - offset, - n_q_points); - mapping.transform( - make_array_view(fe_data.untransformed_shape_grads, - offset, - n_q_points), - mapping_covariant_gradient, - mapping_internal, - transformed_shape_grads); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - transformed_shape_grads[k][d] -= - output_data.shape_values(first + n, k) * - mapping_data.jacobian_pushed_forward_grads[k][n][d]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_gradients[first + d][k] = - dof_sign * transformed_shape_grads[k][d]; - - break; - } - - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } - - if (fe_data.update_each & update_hessians) - { - switch (mapping_kind) - { - case mapping_none: - { - const ArrayView> - transformed_shape_hessians = - make_array_view(fe_data.transformed_shape_hessians, - offset, - n_q_points); - mapping.transform(make_array_view(fe_data.shape_grad_grads, - dof_index, - offset, - n_q_points), - mapping_covariant_gradient, - mapping_internal, - transformed_shape_hessians); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - transformed_shape_hessians[k][d] -= - output_data.shape_gradients[first + d][k][n] * - mapping_data.jacobian_pushed_forward_grads[k][n]; - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - transformed_shape_hessians[k][d]; - - break; - } - case mapping_covariant: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k + offset] = - fe_data.shape_grad_grads[dof_index][k + offset]; - - const ArrayView> - transformed_shape_hessians = - make_array_view(fe_data.transformed_shape_hessians, - offset, - n_q_points); - mapping.transform( - make_array_view(fe_data.untransformed_shape_hessian_tensors, - offset, - n_q_points), - mapping_covariant_hessian, - mapping_internal, - transformed_shape_hessians); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][d][i][j]) + - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][d][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][d]); - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - transformed_shape_hessians[k][d]; - - break; - } - - case mapping_contravariant: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k + offset] = - fe_data.shape_grad_grads[dof_index][k + offset]; - - const ArrayView> - transformed_shape_hessians = - make_array_view(fe_data.transformed_shape_hessians, - offset, - n_q_points); - mapping.transform( - make_array_view(fe_data.untransformed_shape_hessian_tensors, - offset, - n_q_points), - mapping_contravariant_hessian, - mapping_internal, - transformed_shape_hessians); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - transformed_shape_hessians[k][d][i][j] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][d][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][d][n][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][d][i][n]) - - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]); - for (unsigned int m = 0; m < spacedim; ++m) - transformed_shape_hessians[k][d][i][j] -= - (mapping_data - .jacobian_pushed_forward_grads[k][d][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + n, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][d][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + n, k)); - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - transformed_shape_hessians[k][d]; - - break; - } - - case mapping_raviart_thomas: - case mapping_piola: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k + offset] = - fe_data.shape_grad_grads[dof_index][k + offset]; - - const ArrayView> - transformed_shape_hessians = - make_array_view(fe_data.transformed_shape_hessians, - offset, - n_q_points); - mapping.transform( - make_array_view(fe_data.untransformed_shape_hessian_tensors, - offset, - n_q_points), - mapping_piola_hessian, - mapping_internal, - transformed_shape_hessians); - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - transformed_shape_hessians[k][d][i][j] += - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][d][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][d][n][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][d][i][n]) - - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]); - - transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + d, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][n][i][j]) + - (output_data.shape_gradients[first + d][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][n][j]) + - (output_data.shape_gradients[first + d][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][n][i]); - - for (unsigned int m = 0; m < spacedim; ++m) - { - transformed_shape_hessians[k][d][i][j] -= - (mapping_data - .jacobian_pushed_forward_grads[k][d][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + n, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][d][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + n, k)); - - transformed_shape_hessians[k][d][i][j] += - (mapping_data - .jacobian_pushed_forward_grads[k][n][i] - [m] * - mapping_data - .jacobian_pushed_forward_grads[k][m][n] - [j] * - output_data.shape_values(first + d, k)) + - (mapping_data - .jacobian_pushed_forward_grads[k][n][m] - [j] * - mapping_data - .jacobian_pushed_forward_grads[k][m][i] - [n] * - output_data.shape_values(first + d, k)); - } - } - - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - dof_sign * transformed_shape_hessians[k][d]; - - break; - } - - case mapping_nedelec: - { - for (unsigned int k = 0; k < n_q_points; ++k) - fe_data.untransformed_shape_hessian_tensors[k + offset] = - fe_data.shape_grad_grads[dof_index][k + offset]; - - const ArrayView> - transformed_shape_hessians = - make_array_view(fe_data.transformed_shape_hessians, - offset, - n_q_points); - mapping.transform( - make_array_view(fe_data.untransformed_shape_hessian_tensors, - offset, - n_q_points), - mapping_covariant_hessian, - mapping_internal, - transformed_shape_hessians); + const CellSimilarity::Similarity /*cell_similarity*/, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + Assert(dynamic_cast(&fe_internal) != nullptr, + ExcInternalError()); - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < spacedim; ++d) - for (unsigned int n = 0; n < spacedim; ++n) - for (unsigned int i = 0; i < spacedim; ++i) - for (unsigned int j = 0; j < spacedim; ++j) - { - transformed_shape_hessians[k][d][i][j] -= - (output_data.shape_values(first + n, k) * - mapping_data - .jacobian_pushed_forward_2nd_derivatives - [k][n][d][i][j]) + - (output_data.shape_gradients[first + d][k][n] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][j]) + - (output_data.shape_gradients[first + n][k][i] * - mapping_data - .jacobian_pushed_forward_grads[k][n][d][j]) + - (output_data.shape_gradients[first + n][k][j] * - mapping_data - .jacobian_pushed_forward_grads[k][n][i][d]); - } + compute_fill(cell, + QProjector::DataSetDescriptor::cell(), + quadrature.size(), + mapping, + mapping_internal, + mapping_data, + static_cast(fe_internal), + output_data); +} - for (unsigned int k = 0; k < n_q_points; ++k) - for (unsigned int d = 0; d < dim; ++d) - output_data.shape_hessians[first + d][k] = - dof_sign * transformed_shape_hessians[k][d]; - break; - } - default: - DEAL_II_NOT_IMPLEMENTED(); - } - } +template +void +FE_PolyTensor::fill_fe_face_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const hp::QCollection &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + AssertDimension(quadrature.size(), 1); + Assert(dynamic_cast(&fe_internal) != nullptr, + ExcInternalError()); - // third derivatives are not implemented - if (fe_data.update_each & update_3rd_derivatives) - { - DEAL_II_NOT_IMPLEMENTED(); - } - } + compute_fill(cell, + QProjector::DataSetDescriptor::face( + this->reference_cell(), + face_no, + // TODO: fix the custom implementation of orientation + dim == 2 ? numbers::default_geometric_orientation : + cell->combined_face_orientation(face_no), + quadrature[0].size()), + quadrature[0].size(), + mapping, + mapping_internal, + mapping_data, + static_cast(fe_internal), + output_data); } @@ -1849,29 +606,47 @@ FE_PolyTensor::fill_fe_subface_values( spacedim> &output_data) const { - // convert data object to internal - // data for this class. fails with - // an exception if that is not - // possible Assert(dynamic_cast(&fe_internal) != nullptr, ExcInternalError()); - const InternalData &fe_data = static_cast(fe_internal); - const unsigned int n_q_points = quadrature.size(); + compute_fill(cell, + QProjector::DataSetDescriptor::subface( + this->reference_cell(), + face_no, + sub_no, + // TODO: fix the custom implementation of orientation + dim == 2 ? numbers::default_geometric_orientation : + cell->combined_face_orientation(face_no), + quadrature.size(), + cell->subface_case(face_no)), + quadrature.size(), + mapping, + mapping_internal, + mapping_data, + static_cast(fe_internal), + output_data); +} + - // offset determines which data set - // to take (all data sets for all - // sub-faces are stored contiguously) - const auto offset = QProjector::DataSetDescriptor::subface( - this->reference_cell(), - face_no, - sub_no, - dim == 2 ? numbers::default_geometric_orientation : - cell->combined_face_orientation(face_no), - n_q_points, - cell->subface_case(face_no)); - // TODO: Size assertions +template +void +FE_PolyTensor::compute_fill( + const typename Triangulation::cell_iterator &cell, + const typename QProjector::DataSetDescriptor &offset, + const unsigned int n_q_points, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const internal::FEValuesImplementation::MappingRelatedData + &mapping_data, + const typename FE_PolyTensor::InternalData &fe_data, + internal::FEValuesImplementation::FiniteElementRelatedData + &output_data) const +{ + Assert(!(fe_data.update_each & update_values) || + fe_data.shape_values.size()[0] == this->n_dofs_per_cell(), + ExcDimensionMismatch(fe_data.shape_values.size()[0], + this->n_dofs_per_cell())); // TODO: The dof_sign_change only affects Nedelec elements and is not the // correct thing on complicated meshes for higher order Nedelec elements. @@ -1989,7 +764,6 @@ FE_PolyTensor::fill_fe_subface_values( break; } - case mapping_raviart_thomas: case mapping_piola: { @@ -1998,7 +772,6 @@ FE_PolyTensor::fill_fe_subface_values( make_array_view(fe_data.transformed_shape_values, offset, n_q_points); - mapping.transform(make_array_view(fe_data.shape_values, dof_index, offset, @@ -2020,7 +793,6 @@ FE_PolyTensor::fill_fe_subface_values( make_array_view(fe_data.transformed_shape_values, offset, n_q_points); - mapping.transform(make_array_view(fe_data.shape_values, dof_index, offset, @@ -2087,7 +859,6 @@ FE_PolyTensor::fill_fe_subface_values( for (unsigned int d = 0; d < dim; ++d) output_data.shape_gradients[first + d][k] = transformed_shape_grads[k][d]; - break; } @@ -2096,7 +867,6 @@ FE_PolyTensor::fill_fe_subface_values( for (unsigned int k = 0; k < n_q_points; ++k) fe_data.untransformed_shape_grads[k + offset] = fe_data.shape_grads[dof_index][k + offset]; - mapping.transform( make_array_view(fe_data.untransformed_shape_grads, offset, @@ -2126,7 +896,6 @@ FE_PolyTensor::fill_fe_subface_values( for (unsigned int k = 0; k < n_q_points; ++k) fe_data.untransformed_shape_grads[k + offset] = fe_data.shape_grads[dof_index][k + offset]; - mapping.transform( make_array_view(fe_data.untransformed_shape_grads, offset, @@ -2155,6 +924,7 @@ FE_PolyTensor::fill_fe_subface_values( case mapping_nedelec: { + // treat the gradients of // this particular shape // function at all // q-points. if Dv is the @@ -2395,6 +1165,7 @@ FE_PolyTensor::fill_fe_subface_values( (output_data.shape_gradients[first + d][k][j] * mapping_data .jacobian_pushed_forward_grads[k][n][n][i]); + for (unsigned int m = 0; m < spacedim; ++m) { transformed_shape_hessians[k][d][i][j] -=