From 2b0457cd91e0ecc40654f5cda0a3954c4cbecc06 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 4 Aug 2017 23:36:04 +0200 Subject: [PATCH] Use tensorized evaluation in MappingQGeneric if supported --- include/deal.II/fe/mapping_q_generic.h | 44 ++++ source/fe/mapping_q_generic.cc | 341 +++++++++++++++++++++---- 2 files changed, 332 insertions(+), 53 deletions(-) diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index a98193885e..15684b7328 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -21,9 +21,11 @@ #include #include #include +#include #include #include #include +#include #include #include @@ -448,6 +450,48 @@ public: */ QGaussLobatto<1> line_support_points; + /** + * In case the quadrature rule given represents a tensor product + * we need to store the evaluations of the 1d polynomials at the + * the 1d quadrature quadrature points. That is what this variable is for. + */ + internal::MatrixFreeFunctions::ShapeInfo> shape_info; + + /** + * In case the quadrature rule given represents a tensor product + * we need to store temporary data in this object. + */ + mutable AlignedVector > scratch; + + /** + * In case the quadrature rule given represents a tensor product + * the values at the mapped support points are stored in this object. + */ + mutable AlignedVector > values_dofs; + + /** + * In case the quadrature rule given represents a tensor product + * the values at the quadrature points are stored in this object. + */ + mutable AlignedVector > values_quad; + + /** + * In case the quadrature rule given represents a tensor product + * the gradients at the quadrature points are stored in this object. + */ + mutable AlignedVector > gradients_quad; + + /** + * In case the quadrature rule given represents a tensor product + * the hessians at the quadrature points are stored in this object. + */ + mutable AlignedVector > hessians_quad; + + /** + * Indicates whether the given Quadrature object is a tensor product. + */ + bool tensor_product_quadrature; + /** * Tensors of covariant transformation at each of the quadrature points. * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} * diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index f1728e381d..f95378bea4 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -35,6 +35,8 @@ #include #include #include +#include +#include #include #include @@ -680,8 +682,57 @@ initialize (const UpdateFlags update_flags, (update_jacobian_3rd_derivatives | update_jacobian_pushed_forward_3rd_derivatives) ) shape_fourth_derivatives.resize(n_shape_functions * n_q_points); + const std::vector > &ref_q_points = q.get_points(); // now also fill the various fields with their correct values - compute_shape_function_values (q.get_points()); + compute_shape_function_values (ref_q_points); + + tensor_product_quadrature = q.is_tensor_product(); + + if (dim>1) + { + // find out if the one-dimensional formula is the same + // in all directions + if (tensor_product_quadrature) + { + const std::array, dim> quad_array = q.get_tensor_basis(); + for (unsigned int i=1; i> &points_1 = quad_array[i-1].get_points(); + const std::vector> &points_2 = quad_array[i].get_points(); + const std::vector &weights_1 = quad_array[i-1].get_weights(); + const std::vector &weights_2 = quad_array[i].get_weights(); + for (unsigned int j=0; j1.e-10 + || std::abs(weights_1[j]-weights_2[j])>1.e-10) + tensor_product_quadrature = false; + break; + } + } + } + + if (tensor_product_quadrature) + { + const FE_Q fe(polynomial_degree); + shape_info.reinit(q.get_tensor_basis()[0], fe); + + const unsigned int n_shape_values = fe.n_dofs_per_cell(); + const unsigned int max_size = std::max(n_q_points,n_shape_values); + const unsigned int vec_length = VectorizedArray::n_array_elements; + const unsigned int n_comp = 1+ (spacedim-1)/vec_length; + + scratch.resize((dim-1)*max_size); + values_dofs.resize(n_comp*n_shape_values); + } + } + } } @@ -695,6 +746,22 @@ initialize_face (const UpdateFlags update_flags, { initialize (update_flags, q, n_original_q_points); + if (dim>1 && tensor_product_quadrature) + { + const unsigned int facedim = dim > 1 ? dim-1 : 1; + const FE_Q fe(polynomial_degree); + shape_info.reinit(q.get_tensor_basis()[0], fe); + + const unsigned int n_shape_values = fe.n_dofs_per_cell(); + const unsigned int n_q_points = q.size(); + const unsigned int max_size = std::max(n_q_points,n_shape_values); + const unsigned int vec_length = VectorizedArray::n_array_elements; + const unsigned int n_comp = 1+ (spacedim-1)/vec_length; + + scratch.resize((dim-1)*max_size); + values_dofs.resize(n_comp*n_shape_values); + } + if (dim > 1) { if (this->update_each & (update_boundary_forms | @@ -1527,20 +1594,69 @@ namespace internal if (update_flags & update_quadrature_points) { - for (unsigned int point=0; point1 && data.tensor_product_quadrature) { - const double *shape = &data.shape(point+data_set,0); - Point result = (shape[0] * - data.mapping_support_points[0]); - for (unsigned int k=1; k 0, ExcInternalError()); + + const unsigned int n_shape_values = data.n_shape_functions; + const unsigned int n_q_points = quadrature_points.size(); + 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 == quadrature_points.size(), + ExcDimensionMismatch(data.shape_info.n_q_points, quadrature_points.size())); + + data.values_dofs.resize(n_comp*n_shape_values); + VectorizedArray *values_dofs_ptr[n_comp]; + data.values_quad.resize(n_comp*n_q_points); + VectorizedArray *values_quad_ptr[n_comp]; + + for (unsigned int c=0; c &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]), nullptr, nullptr, + &(data.scratch[0]), true, false, false); + + for (unsigned int out_comp=0; out_comp result = (shape[0] * + data.mapping_support_points[0]); + for (unsigned int k=1; k()); Assert (data.n_shape_functions > 0, ExcInternalError()); - const Tensor<1,spacedim> *supp_pts = - &data.mapping_support_points[0]; - for (unsigned int point=0; point1 && data.tensor_product_quadrature) { - const Tensor<1,dim> *data_derv = - &data.derivative(point+data_set, 0); + 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]; + + // 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]), nullptr, &(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]; + + for (unsigned int point=0; point *data_derv = + &data.derivative(point+data_set, 0); + + 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 Tensor<2,dim> *second = - &data.second_derivative(point+data_set, 0); - double result [spacedim][dim][dim]; - for (unsigned int i=0; i::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]; + + // 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]), nullptr, nullptr, &(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