From 1e6c934750958ca4e7a89200cf46fd54706e5184 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 22 Aug 2017 19:00:17 +0200 Subject: [PATCH] Use one function for maybe_compute/update_* --- source/fe/mapping_q_generic.cc | 570 ++++++++++++++++----------------- 1 file changed, 276 insertions(+), 294 deletions(-) diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 08c6b24fbb..72ab60eaa9 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -1583,89 +1583,182 @@ namespace internal } /** - * Compute the locations of quadrature points on the object described by - * the first argument (and the cell for which the mapping support points - * have already been set), but only if the update_flags of the @p data - * argument indicate so. + * In case the quadrature formula is a tensor product, this is a replacement + * for maybe_compute_q_points(), maybe_update_Jacobians() and + * maybe_update_jacobian_grads() */ template void - maybe_compute_q_points - (const typename QProjector::DataSetDescriptor data_set, + maybe_update_q_points_Jacobians_and_grads_tensor + (const CellSimilarity::Similarity cell_similarity, const typename dealii::MappingQGeneric::InternalData &data, - std::vector > &quadrature_points) + std::vector > &quadrature_points, + std::vector > &jacobian_grads) { const UpdateFlags update_flags = data.update_each; - if (update_flags & update_quadrature_points) + const unsigned int n_shape_values = data.n_shape_functions; + const unsigned int n_q_points = data.shape_info.n_q_points; + const unsigned int vec_length = VectorizedArray::n_array_elements; + const unsigned int n_comp = 1+ (spacedim-1)/vec_length; + const unsigned int n_hessians = (dim*(dim+1))/2; + + VectorizedArray *values_dofs_ptr[n_comp]; + VectorizedArray *values_quad_ptr[n_comp]; + VectorizedArray *gradients_quad_ptr[n_comp][dim]; + VectorizedArray *hessians_quad_ptr[n_comp][n_hessians]; + + const bool evaluate_values = update_flags & update_quadrature_points; + const bool evaluate_gradients= (cell_similarity != CellSimilarity::translation) + &&(update_flags & update_contravariant_transformation); + const bool evaluate_hessians = (cell_similarity != CellSimilarity::translation) + &&(update_flags & update_jacobian_grads); + + Assert (!evaluate_values || n_q_points > 0, ExcInternalError()); + Assert (!evaluate_values || n_q_points == quadrature_points.size(), + ExcDimensionMismatch(n_q_points, quadrature_points.size())); + Assert (!evaluate_gradients || data.n_shape_functions > 0, ExcInternalError()); + Assert (!evaluate_gradients || n_q_points == data.contravariant.size(), + ExcDimensionMismatch(n_q_points, data.contravariant.size())); + Assert (!evaluate_hessians || n_q_points == jacobian_grads.size(), + ExcDimensionMismatch(n_q_points, jacobian_grads.size())); + + // prepare arrays + if (evaluate_values || evaluate_gradients || evaluate_hessians) { - if (dim>1 && data.tensor_product_quadrature) - { - Assert(data.shape_info.n_q_points > 0, ExcInternalError()); + data.values_dofs.resize(n_comp*n_shape_values); + data.values_quad.resize(n_comp*n_q_points); + data.gradients_quad.resize (n_comp*n_q_points*dim); + + const std::vector &renumber_to_lexicographic + = data.shape_info.lexicographic_numbering; + for (unsigned int i=0; i::n_array_elements; - const unsigned int n_comp = 1+ (spacedim-1)/vec_length; + for (unsigned int c=0; c *values_dofs_ptr[n_comp]; - data.values_quad.resize(n_comp*n_q_points); - VectorizedArray *values_quad_ptr[n_comp]; - // Some evaluators need to write into the gradients array. - data.gradients_quad.resize (n_comp*n_q_points*dim); - VectorizedArray *gradients_quad_ptr[n_comp][dim]; + // do the actual tensorized evaluation + SelectEvaluator::evaluate + (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]), + &(gradients_quad_ptr[0]), &(hessians_quad_ptr[0]), &(data.scratch[0]), + evaluate_values, evaluate_gradients, evaluate_hessians); + } - for (unsigned int c=0; c &renumber_to_lexicographic - = data.shape_info.lexicographic_numbering; - for (unsigned int i=0; i()); + // We need to reinterpret the data after evaluate has been applied. + for (unsigned int out_comp=0; out_comp::evaluate - (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]), - &(gradients_quad_ptr[0]), nullptr, &(data.scratch[0]), true, false, false); + if (update_flags & update_volume_elements) + if (cell_similarity != CellSimilarity::translation) + for (unsigned int point=0; point result = (shape[0] * - data.mapping_support_points[0]); - for (unsigned int k=1; k + void + maybe_compute_q_points + (const typename QProjector::DataSetDescriptor data_set, + const typename dealii::MappingQGeneric::InternalData &data, + std::vector > &quadrature_points) + { + const UpdateFlags update_flags = data.update_each; + + if (update_flags & update_quadrature_points) + for (unsigned int point=0; point result = (shape[0] * + data.mapping_support_points[0]); + for (unsigned int k=1; k 0, ExcInternalError()); - if (dim>1 && data.tensor_product_quadrature) - { - const unsigned int n_shape_values = data.n_shape_functions; - const unsigned int vec_length = VectorizedArray::n_array_elements; - const unsigned int n_comp = 1+ (spacedim-1)/vec_length; - - Assert (data.shape_info.n_q_points == data.contravariant.size(), - ExcDimensionMismatch(data.shape_info.n_q_points, data.contravariant.size())); - - data.values_dofs.resize(n_comp*n_shape_values); - VectorizedArray *values_dofs_ptr[n_comp]; - data.gradients_quad.resize (n_comp*n_q_points*dim); - VectorizedArray *gradients_quad_ptr[n_comp][dim]; - // Some evaluators need to write into the values array. - data.values_quad.resize(n_comp*n_q_points); - VectorizedArray *values_quad_ptr[n_comp]; - - // transform data appropriately - const std::vector &renumber_to_lexicographic - = data.shape_info.lexicographic_numbering; - for (unsigned int i=0; i *supp_pts = + &data.mapping_support_points[0]; - for (unsigned int c=0; c::evaluate - (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]), - &(gradients_quad_ptr[0]), nullptr, &(data.scratch[0]), false, true, false); - - // We need to reinterpret the data after evaluate has been applied. - for (unsigned int out_comp=0; out_comp 0, ExcInternalError()); - const Tensor<1,spacedim> *supp_pts = - &data.mapping_support_points[0]; + const Tensor<1,dim> *data_derv = + &data.derivative(point+data_set, 0); - for (unsigned int point=0; point *data_derv = - &data.derivative(point+data_set, 0); - - double result [spacedim][dim]; + double result [spacedim][dim]; - // peel away part of sum to avoid zeroing the - // entries and adding for the first time - for (unsigned int i=0; i1 && data.tensor_product_quadrature) - { - const unsigned int n_shape_values = data.n_shape_functions; - const unsigned int vec_length = VectorizedArray::n_array_elements; - const unsigned int n_comp = 1+ (spacedim-1)/vec_length; - const unsigned int n_hessians = (dim*(dim+1))/2; - - Assert (data.shape_info.n_q_points == jacobian_grads.size(), - ExcDimensionMismatch(data.shape_info.n_q_points, jacobian_grads.size())); - - data.values_dofs.resize(n_comp*n_shape_values); - VectorizedArray *values_dofs_ptr[n_comp]; - data.hessians_quad.resize(n_comp*n_q_points*n_hessians); - VectorizedArray *hessians_quad_ptr[n_comp][n_hessians]; - // Some evaluators need to write into the gradients array - // and into the values array. - data.gradients_quad.resize (n_comp*n_q_points*dim); - VectorizedArray *gradients_quad_ptr[n_comp][dim]; - data.values_quad.resize(n_comp*n_q_points); - VectorizedArray *values_quad_ptr[n_comp]; - - // transform data appropriately - const std::vector &renumber_to_lexicographic - = data.shape_info.lexicographic_numbering; - for (unsigned int i=0; i::evaluate - (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]), - &(gradients_quad_ptr[0]), &(hessians_quad_ptr[0]), - &(data.scratch[0]), false, false, true); - - constexpr int desymmetrize_3d [6][2] = {{0,0},{1,1},{2,2},{0,1},{0,2},{1,2}}; - constexpr int desymmetrize_2d [3][2] = {{0,0},{1,1},{0,1}}; - - // We need to reinterpret the data after evaluate has been applied. - for (unsigned int out_comp=0; out_comp *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i *second = + &data.second_derivative(point+data_set, 0); + double result [spacedim][dim][dim]; + for (unsigned int i=0; i::cell_iterator &cell, const CellSimilarity::Similarity computed_cell_similarity = (polynomial_degree == 1 ? cell_similarity : CellSimilarity::none); - internal::MappingQGeneric::maybe_compute_q_points - (QProjector::DataSetDescriptor::cell (), + if (dim>1 && data.tensor_product_quadrature) + { + internal::MappingQGeneric::maybe_update_q_points_Jacobians_and_grads_tensor + (computed_cell_similarity, + data, + output_data.quadrature_points, + output_data.jacobian_grads); + } + else + { + internal::MappingQGeneric::maybe_compute_q_points + (QProjector::DataSetDescriptor::cell (), + data, + output_data.quadrature_points); + + internal::MappingQGeneric::maybe_update_Jacobians + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data); + + internal::MappingQGeneric::maybe_update_jacobian_grads + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_grads); + } + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_grads + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_grads); + + internal::MappingQGeneric::maybe_update_jacobian_2nd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_2nd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_2nd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), + data, + output_data.jacobian_pushed_forward_2nd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_3rd_derivatives + (computed_cell_similarity, + QProjector::DataSetDescriptor::cell (), data, - output_data.quadrature_points); - internal::MappingQGeneric::maybe_update_Jacobians + output_data.jacobian_3rd_derivatives); + + internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_3rd_derivatives (computed_cell_similarity, QProjector::DataSetDescriptor::cell (), - data); + data, + output_data.jacobian_pushed_forward_3rd_derivatives); const UpdateFlags update_flags = data.update_each; const std::vector &weights=quadrature.get_weights(); @@ -2962,42 +2969,6 @@ fill_fe_values (const typename Triangulation::cell_iterator &cell, output_data.inverse_jacobians[point] = data.covariant[point].transpose(); } - internal::MappingQGeneric::maybe_update_jacobian_grads - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_grads); - - internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_grads - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_grads); - - internal::MappingQGeneric::maybe_update_jacobian_2nd_derivatives - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_2nd_derivatives); - - internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_2nd_derivatives - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_2nd_derivatives); - - internal::MappingQGeneric::maybe_update_jacobian_3rd_derivatives - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_3rd_derivatives); - - internal::MappingQGeneric::maybe_update_jacobian_pushed_forward_3rd_derivatives - (computed_cell_similarity, - QProjector::DataSetDescriptor::cell (), - data, - output_data.jacobian_pushed_forward_3rd_derivatives); - return computed_cell_similarity; } @@ -3183,16 +3154,27 @@ namespace internal const typename dealii::MappingQGeneric::InternalData &data, internal::FEValues::MappingRelatedData &output_data) { - maybe_compute_q_points (data_set, - data, - output_data.quadrature_points); - maybe_update_Jacobians (CellSimilarity::none, - data_set, - data); - maybe_update_jacobian_grads (CellSimilarity::none, - data_set, - data, - output_data.jacobian_grads); + if (dim>1 && data.tensor_product_quadrature) + { + maybe_update_q_points_Jacobians_and_grads_tensor + (CellSimilarity::none, + data, + output_data.quadrature_points, + output_data.jacobian_grads); + } + else + { + maybe_compute_q_points (data_set, + data, + output_data.quadrature_points); + maybe_update_Jacobians (CellSimilarity::none, + data_set, + data); + maybe_update_jacobian_grads (CellSimilarity::none, + data_set, + data, + output_data.jacobian_grads); + } maybe_update_jacobian_pushed_forward_grads (CellSimilarity::none, data_set, data, -- 2.39.5