From 4d2e2f4f1749c169452486b8172548bff7327d65 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 3 Aug 2017 23:16:35 +0200 Subject: [PATCH] Do not hardcode VectorizedArray in MatrixFreeFunctions::ShapeInfo. Make class templated on generic type Number that can be VectorizedArray. --- .../matrix_free/cuda_matrix_free.templates.h | 4 +- .../deal.II/matrix_free/evaluation_kernels.h | 24 ++--- include/deal.II/matrix_free/fe_evaluation.h | 16 ++-- include/deal.II/matrix_free/matrix_free.h | 22 ++--- include/deal.II/matrix_free/shape_info.h | 89 +++++++++++++----- .../matrix_free/shape_info.templates.h | 91 +++++++++++-------- source/matrix_free/matrix_free.cc | 3 + 7 files changed, 154 insertions(+), 95 deletions(-) diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h index 8c93531625..b90a7bb3db 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -461,7 +461,7 @@ namespace CUDAWrappers unsigned int size_shape_values = n_dofs_1d*n_q_points_1d*sizeof(Number); cudaError_t cuda_error = cudaMemcpyToSymbol(internal::global_shape_values, - &shape_info.shape_values_number[0], + &shape_info.shape_values[0], size_shape_values, 0, cudaMemcpyHostToDevice); @@ -470,7 +470,7 @@ namespace CUDAWrappers if (update_flags & update_gradients) { cuda_error = cudaMemcpyToSymbol(internal::global_shape_gradients, - &shape_info.shape_gradient_number[0], + &shape_info.shape_gradient[0], size_shape_values, 0, cudaMemcpyHostToDevice); diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 13de640678..d7cb55e021 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -100,7 +100,7 @@ namespace internal struct FEEvaluationImpl { static - void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs_actual[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -111,7 +111,7 @@ namespace internal const bool evaluate_hessians); static - void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs_actual[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -126,7 +126,7 @@ namespace internal inline void FEEvaluationImpl - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs_actual[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -351,7 +351,7 @@ namespace internal inline void FEEvaluationImpl - ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs_actual[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -554,7 +554,7 @@ namespace internal struct FEEvaluationImplTransformToCollocation { static - void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -565,7 +565,7 @@ namespace internal const bool evaluate_hessians); static - void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -578,7 +578,7 @@ namespace internal inline void FEEvaluationImplTransformToCollocation - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -662,7 +662,7 @@ namespace internal inline void FEEvaluationImplTransformToCollocation - ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -741,7 +741,7 @@ namespace internal struct FEEvaluationImplCollocation { static - void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -752,7 +752,7 @@ namespace internal const bool evaluate_hessians); static - void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + void integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -765,7 +765,7 @@ namespace internal inline void FEEvaluationImplCollocation - ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], @@ -829,7 +829,7 @@ namespace internal inline void FEEvaluationImplCollocation - ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, VectorizedArray *values_dofs[], VectorizedArray *values_quad[], VectorizedArray *gradients_quad[][dim], diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 3058fb2891..eb000af737 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -153,8 +153,8 @@ public: /** * Return a reference to the ShapeInfo object currently in use. */ - const internal::MatrixFreeFunctions::ShapeInfo & - get_shape_info() const; + const internal::MatrixFreeFunctions::ShapeInfo> & + get_shape_info() const; /** * Fills the JxW values currently used. @@ -808,7 +808,7 @@ protected: * product. Also contained in matrix_info, but it simplifies code if we * store a reference to it. */ - const internal::MatrixFreeFunctions::ShapeInfo *data; + const internal::MatrixFreeFunctions::ShapeInfo> *data; /** * A pointer to the Cartesian Jacobian information of the present cell. Only @@ -2168,7 +2168,7 @@ FEEvaluationBase dof_info (nullptr), mapping_info (nullptr), // select the correct base element from the given FE component - data (new internal::MatrixFreeFunctions::ShapeInfo(quadrature, fe, fe.component_to_base_index(first_selected_component).first)), + data (new internal::MatrixFreeFunctions::ShapeInfo>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)), cartesian_data (nullptr), jacobian (nullptr), J_value (nullptr), @@ -2230,7 +2230,7 @@ FEEvaluationBase dof_info (other.dof_info), mapping_info (other.mapping_info), data (other.matrix_info == nullptr ? - new internal::MatrixFreeFunctions::ShapeInfo(*other.data) : + new internal::MatrixFreeFunctions::ShapeInfo>(*other.data) : other.data), cartesian_data (nullptr), jacobian (nullptr), @@ -2301,7 +2301,7 @@ FEEvaluationBase mapping_info = other.mapping_info; if (other.matrix_info == nullptr) { - data = new internal::MatrixFreeFunctions::ShapeInfo(*other.data); + data = new internal::MatrixFreeFunctions::ShapeInfo>(*other.data); scratch_data_array = new AlignedVector >(); } else @@ -2552,8 +2552,8 @@ FEEvaluationBase::get_cell_type () const template inline -const internal::MatrixFreeFunctions::ShapeInfo & -FEEvaluationBase::get_shape_info() const +const internal::MatrixFreeFunctions::ShapeInfo> & + FEEvaluationBase::get_shape_info() const { Assert(data != nullptr, ExcInternalError()); return *data; diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 8a9d79b752..7f1c9af072 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -841,11 +841,11 @@ public: /** * Return the unit cell information for given hp index. */ - const internal::MatrixFreeFunctions::ShapeInfo & - get_shape_info (const unsigned int fe_component = 0, - const unsigned int quad_index = 0, - const unsigned int hp_active_fe_index = 0, - const unsigned int hp_active_quad_index = 0) const; + const internal::MatrixFreeFunctions::ShapeInfo> & + get_shape_info (const unsigned int fe_component = 0, + const unsigned int quad_index = 0, + const unsigned int hp_active_fe_index = 0, + const unsigned int hp_active_quad_index = 0) const; /** * Obtains a scratch data object for internal use. Make sure to release it @@ -980,7 +980,7 @@ private: /** * Contains shape value information on the unit cell. */ - Table<4,internal::MatrixFreeFunctions::ShapeInfo > shape_info; + Table<4,internal::MatrixFreeFunctions::ShapeInfo>> shape_info; /** * Describes how the cells are gone through. With the cell level (first @@ -1452,11 +1452,11 @@ MatrixFree::get_ghost_set(const unsigned int dof_index) const template inline -const internal::MatrixFreeFunctions::ShapeInfo & -MatrixFree::get_shape_info (const unsigned int index_fe, - const unsigned int index_quad, - const unsigned int active_fe_index, - const unsigned int active_quad_index) const +const internal::MatrixFreeFunctions::ShapeInfo> & + MatrixFree::get_shape_info (const unsigned int index_fe, + const unsigned int index_quad, + const unsigned int active_fe_index, + const unsigned int active_quad_index) const { AssertIndexRange (index_fe, shape_info.size(0)); AssertIndexRange (index_quad, shape_info.size(1)); diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index e89de80f49..424ebd69b3 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -20,7 +20,6 @@ #include #include -#include #include #include @@ -132,7 +131,7 @@ namespace internal * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. */ - AlignedVector > shape_values; + AlignedVector shape_values; /** * Stores the shape gradients of the 1D finite element evaluated on all @@ -141,7 +140,7 @@ namespace internal * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. */ - AlignedVector > shape_gradients; + AlignedVector shape_gradients; /** * Stores the shape Hessians of the 1D finite element evaluated on all @@ -150,67 +149,67 @@ namespace internal * this array is n_dofs_1d * n_q_points_1d and quadrature * points are the index running fastest. */ - AlignedVector > shape_hessians; + AlignedVector shape_hessians; /** * Stores the shape values in a different format, namely the so-called * even-odd scheme where the symmetries in shape_values are used for * faster evaluation. */ - AlignedVector > shape_values_eo; + AlignedVector shape_values_eo; /** * Stores the shape gradients in a different format, namely the so- * called even-odd scheme where the symmetries in shape_gradients are * used for faster evaluation. */ - AlignedVector > shape_gradients_eo; + AlignedVector shape_gradients_eo; /** * Stores the shape second derivatives in a different format, namely the * so-called even-odd scheme where the symmetries in shape_hessians are * used for faster evaluation. */ - AlignedVector > shape_hessians_eo; + AlignedVector shape_hessians_eo; /** * Stores the shape gradients of the shape function space associated to * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For * faster evaluation only the even-odd format is necessary. */ - AlignedVector > shape_gradients_collocation_eo; + AlignedVector shape_gradients_collocation_eo; /** * Stores the shape hessians of the shape function space associated to * the quadrature (collocation), given by FE_DGQ<1>(Quadrature<1>). For * faster evaluation only the even-odd format is necessary. */ - AlignedVector > shape_hessians_collocation_eo; + AlignedVector shape_hessians_collocation_eo; /** * Collects all data of 1D shape values evaluated at the point 0 and 1 * (the vertices) in one data structure. Sorting is first the values, * then gradients, then second derivatives. */ - AlignedVector > shape_data_on_face[2]; + AlignedVector shape_data_on_face[2]; /** * Stores one-dimensional values of shape functions on subface. Since * there are two subfaces, store two variants. */ - AlignedVector > values_within_subface[2]; + AlignedVector values_within_subface[2]; /** * Stores one-dimensional gradients of shape functions on subface. Since * there are two subfaces, store two variants. */ - AlignedVector > gradients_within_subface[2]; + AlignedVector gradients_within_subface[2]; /** * Stores one-dimensional gradients of shape functions on subface. Since * there are two subfaces, store two variants. */ - AlignedVector > hessians_within_subface[2]; + AlignedVector hessians_within_subface[2]; /** * Renumbering from deal.II's numbering of cell degrees of freedom to @@ -259,31 +258,77 @@ namespace internal bool nodal_at_cell_boundaries; /** - * For nodal cells, we might get around by simply loading the indices to - * the degrees of freedom that act on a particular face, rather than the - * whole set of indices that is then interpolated down to the - * element. This array stores this indirect addressing. + * For nodal basis functions with nodes located at the boundary of the + * unit cell, face integrals that involve only the values of the shape + * functions (approximations of first derivatives in DG) do not need to + * load all degrees of freedom of the cell but rather only the degrees + * of freedom located on the face. While it would also be possible to + * compute these indices on the fly, we choose to simplify the code and + * store the indirect addressing in this class. * * The first table index runs through the faces of a cell, and the * second runs through the nodal degrees of freedom of the face, using * @p dofs_per_face entries. * + * The indices stored in this member variable are as follows. Consider + * for example a 2D element of degree 3 with the following degrees of + * freedom in lexicographic numbering: + * @code + * 12 13 14 15 + * 8 9 10 11 + * 4 5 6 7 + * 0 1 2 3 + * @endcode + * + * The first row stores the indices on the face with index 0, i.e., the + * numbers 0, 4, 8, 12, the second row holds the indices + * 3, 7, 11, 15 for face 1, the third row holds the indices + * 0, 1, 2, 3 for face 2, and the last (fourth) row holds + * the indices 12, 13, 14, 15. Similarly, the indices are + * stored in 3D. (Note that the y faces in 3D use indices reversed in + * terms of the lexicographic numbers due to the orientation of the + * coordinate system.) + * * @note This object is only filled in case @p nodal_at_cell_boundaries * evaluates to @p true. */ dealii::Table<2,unsigned int> face_to_cell_index_nodal; /** - * For Hermite-type basis functions, the @p face_to_cell_index_nodal for - * the values on a face of the cell is used together with a respective - * slot in the derivative. In the lexicographic ordering, this index is - * in the next "layer" of degrees of freedom. This array stores the - * indirect addressing of both the values and the gradient. + * The @p face_to_cell_index_nodal provides a shortcut for the + * evaluation of values on the faces. For Hermite-type basis functions, + * one can go one step further and also use shortcuts to get derivatives + * more cheaply where only two layers of degrees of freedom contribute + * to the derivative on the face. In the lexicographic ordering, the + * respective indices is in the next "layer" of degrees of freedom as + * compared to the nodal values. This array stores the indirect + * addressing of both the values and the gradient. * * The first table index runs through the faces of a cell, and the * second runs through the pairs of the nodal degrees of freedom of the * face and the derivatives, using 2*dofs_per_face entries. * + * The indices stored in this member variable are as follows. Consider + * for example a 2D element of degree 3 with the following degrees of + * freedom in lexicographic numbering: + * @code + * 20 21 22 23 24 + * 15 16 17 18 19 + * 10 11 12 13 14 + * 5 6 7 8 9 + * 0 1 2 3 4 + * @endcode + * + * The first row stores the indices for values and gradients on the face + * with index 0, i.e., the numbers 0, 1, 5, 6, 10, 11, 15, 16, 20, + * 21, the second row holds the indices 4, 3, 9, 8, 14, 13, + * 19, 18, 24, 23 for face 1, the third row holds the indices + * 0, 5, 1, 6, 2, 7, 3, 8, 4, 9 for face 2, and the last + * (fourth) row holds the indices 20, 15, 21, 16, 22, 17, 23, 18, + * 24, 19. Similarly, the indices are stored in 3D. (Note that + * the y faces in 3D use indices reversed in terms of the lexicographic + * numbers due to the orientation of the coordinate system.) + * * @note This object is only filled in case @p element_type evaluates to * @p tensor_symmetric_hermite. */ diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index cefbd43e06..9bc6bea4fe 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -40,6 +40,21 @@ namespace internal // ----------------- actual ShapeInfo functions -------------------- + namespace + { + template + Number get_first_array_element(const Number a) + { + return a; + } + + template + Number get_first_array_element(const VectorizedArray a) + { + return a[0]; + } + } + template ShapeInfo::ShapeInfo () : @@ -159,8 +174,8 @@ namespace internal if (fe->has_support_points()) unit_point = fe->get_unit_support_points()[scalar_lexicographic[0]]; Assert(fe->dofs_per_cell == 0 || - std::fabs(fe->shape_value(scalar_lexicographic[0], - unit_point)-1) < 1e-13, + std::abs(fe->shape_value(scalar_lexicographic[0], + unit_point)-1) < 1e-13, ExcInternalError("Could not decode 1D shape functions for the " "element " + fe->get_name())); } @@ -191,9 +206,6 @@ namespace internal const unsigned int my_i = scalar_lexicographic[i]; for (unsigned int q=0; q::n_array_elements copies for the - // shape information and non-vectorized fields Point q_point = unit_point; q_point[0] = quad.get_points()[q][0]; @@ -271,10 +283,10 @@ namespace internal // check if we are a Hermite type element_type = tensor_symmetric_hermite; for (unsigned int i=1; ishape_data_on_face[0][i][0]) > 1e-12) + if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-12) element_type = tensor_symmetric; for (unsigned int i=2; ishape_data_on_face[0][n_dofs_1d+i][0]) > 1e-12) + if (std::abs(get_first_array_element(shape_data_on_face[0][n_dofs_1d+i])) > 1e-12) element_type = tensor_symmetric; } } @@ -283,8 +295,8 @@ namespace internal nodal_at_cell_boundaries = true; for (unsigned int i=1; ishape_data_on_face[0][i][0]) > 1e-13 || - std::abs(this->shape_data_on_face[1][i-1][0]) > 1e-13) + if (std::abs(get_first_array_element(shape_data_on_face[0][i])) > 1e-13 || + std::abs(get_first_array_element(shape_data_on_face[1][i-1])) > 1e-13) nodal_at_cell_boundaries = false; if (nodal_at_cell_boundaries == true) @@ -371,11 +383,10 @@ namespace internal const unsigned int n_dofs_1d = fe_degree + 1; for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; j + if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j] - + shape_values[(n_dofs_1d-i)*n_q_points_1d-j-1])) > std::max(zero_tol, zero_tol* - std::abs(shape_values[i*n_q_points_1d+j][0]))) + std::abs(get_first_array_element(shape_values[i*n_q_points_1d+j])))) return false; // shape values should be zero at x=0.5 for all basis functions except @@ -383,11 +394,11 @@ namespace internal if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1) { for (unsigned int i=0; i zero_tol) + if (std::abs(get_first_array_element(shape_values[i*n_q_points_1d+ + n_q_points_1d/2])) > zero_tol) return false; - if (std::fabs(shape_values[(n_dofs_1d/2)*n_q_points_1d+ - n_q_points_1d/2][0]-1.)> zero_tol) + if (std::abs(get_first_array_element(shape_values[(n_dofs_1d/2)*n_q_points_1d+ + n_q_points_1d/2])-1.)> zero_tol) return false; } @@ -397,13 +408,14 @@ namespace internal const double zero_tol_gradient = zero_tol * std::sqrt(fe_degree+1.)*(fe_degree+1); for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; j zero_tol_gradient) + if (std::abs(get_first_array_element(shape_gradients[i*n_q_points_1d+j] + + shape_gradients[(n_dofs_1d-i)*n_q_points_1d- + j-1])) > zero_tol_gradient) return false; if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1) - if (std::fabs(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ - (n_q_points_1d/2)][0]) > zero_tol_gradient) + if (std::abs(get_first_array_element(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ + (n_q_points_1d/2)])) + > zero_tol_gradient) return false; // symmetry for Hessian. Multiply tolerance by degree^3 of the element @@ -411,9 +423,9 @@ namespace internal const double zero_tol_hessian = zero_tol * (fe_degree+1)*(fe_degree+1)*(fe_degree+1); for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) for (unsigned int j=0; j zero_tol_hessian) + if (std::abs(get_first_array_element(shape_hessians[i*n_q_points_1d+j] - + shape_hessians[(n_dofs_1d-i)*n_q_points_1d- + j-1])) > zero_tol_hessian) return false; const unsigned int stride = (n_q_points_1d+1)/2; @@ -424,25 +436,25 @@ namespace internal for (unsigned int q=0; qzero_tol) + if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j]))>zero_tol) return false; } else { - if (std::fabs(shape_values[i*n_points_1d+ - j][0]-1.)>zero_tol) + if (std::abs(get_first_array_element(shape_values[i*n_points_1d+j])-1.)>zero_tol) return false; } return true; diff --git a/source/matrix_free/matrix_free.cc b/source/matrix_free/matrix_free.cc index d6a07ec0f4..ea2c2e421d 100644 --- a/source/matrix_free/matrix_free.cc +++ b/source/matrix_free/matrix_free.cc @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -26,5 +27,7 @@ DEAL_II_NAMESPACE_OPEN template struct internal::MatrixFreeFunctions::ShapeInfo; template struct internal::MatrixFreeFunctions::ShapeInfo; +template struct internal::MatrixFreeFunctions::ShapeInfo>; +template struct internal::MatrixFreeFunctions::ShapeInfo>; DEAL_II_NAMESPACE_CLOSE -- 2.39.5