From: Julius Witte Date: Thu, 12 Mar 2020 15:56:30 +0000 (+0100) Subject: refactoring ShapeInfo's 1D data into UnivariateShapeData X-Git-Tag: v9.2.0-rc1~349^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3928bfb1ee7f88444301414473a1f86d0195d749;p=dealii.git refactoring ShapeInfo's 1D data into UnivariateShapeData --- 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 dc91f7d444..df97cbe8e1 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -875,7 +875,7 @@ namespace CUDAWrappers cudaError_t cuda_error = cudaMemcpyToSymbol(internal::get_global_shape_values(), - shape_info.shape_values.data(), + shape_info.data.front().shape_values.data(), size_shape_values, 0, cudaMemcpyHostToDevice); @@ -885,7 +885,7 @@ namespace CUDAWrappers { cuda_error = cudaMemcpyToSymbol(internal::get_global_shape_gradients(), - shape_info.shape_gradients.data(), + shape_info.data.front().shape_gradients.data(), 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 79623b3e8a..d559935ff2 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -163,14 +163,17 @@ namespace internal fe_degree + 1, n_q_points_1d, Number>; - Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo : - shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_gradients_eo : - shape_info.shape_gradients, - variant == evaluate_evenodd ? shape_info.shape_hessians_eo : - shape_info.shape_hessians, - shape_info.fe_degree + 1, - shape_info.n_q_points_1d); + Eval eval(variant == evaluate_evenodd ? + shape_info.data.front().shape_values_eo : + shape_info.data.front().shape_values, + variant == evaluate_evenodd ? + shape_info.data.front().shape_gradients_eo : + shape_info.data.front().shape_gradients, + variant == evaluate_evenodd ? + shape_info.data.front().shape_hessians_eo : + shape_info.data.front().shape_hessians, + shape_info.data.front().fe_degree + 1, + shape_info.data.front().n_q_points_1d); const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? @@ -183,9 +186,10 @@ namespace internal if (temp_size == 0) { temp1 = scratch_data; - temp2 = temp1 + - std::max(Utilities::fixed_power(shape_info.fe_degree + 1), - Utilities::fixed_power(shape_info.n_q_points_1d)); + temp2 = temp1 + std::max(Utilities::fixed_power( + shape_info.data.front().fe_degree + 1), + Utilities::fixed_power( + shape_info.data.front().n_q_points_1d)); } else { @@ -197,7 +201,7 @@ namespace internal temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product; const unsigned int dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ? - Utilities::fixed_power(shape_info.fe_degree + 1) : + Utilities::fixed_power(shape_info.data.front().fe_degree + 1) : shape_info.dofs_per_component_on_cell; const Number *values_dofs = values_dofs_actual; if (type == MatrixFreeFunctions::truncated_tensor) @@ -205,7 +209,8 @@ namespace internal Number *values_dofs_tmp = scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell, shape_info.n_q_points)); - const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree; + const int degree = + fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; unsigned int count_p = 0, count_q = 0; for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i) { @@ -434,14 +439,17 @@ namespace internal fe_degree + 1, n_q_points_1d, Number>; - Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo : - shape_info.shape_values, - variant == evaluate_evenodd ? shape_info.shape_gradients_eo : - shape_info.shape_gradients, - variant == evaluate_evenodd ? shape_info.shape_hessians_eo : - shape_info.shape_hessians, - shape_info.fe_degree + 1, - shape_info.n_q_points_1d); + Eval eval(variant == evaluate_evenodd ? + shape_info.data.front().shape_values_eo : + shape_info.data.front().shape_values, + variant == evaluate_evenodd ? + shape_info.data.front().shape_gradients_eo : + shape_info.data.front().shape_gradients, + variant == evaluate_evenodd ? + shape_info.data.front().shape_hessians_eo : + shape_info.data.front().shape_hessians, + shape_info.data.front().fe_degree + 1, + shape_info.data.front().n_q_points_1d); const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? @@ -454,9 +462,10 @@ namespace internal if (temp_size == 0) { temp1 = scratch_data; - temp2 = temp1 + - std::max(Utilities::fixed_power(shape_info.fe_degree + 1), - Utilities::fixed_power(shape_info.n_q_points_1d)); + temp2 = temp1 + std::max(Utilities::fixed_power( + shape_info.data.front().fe_degree + 1), + Utilities::fixed_power( + shape_info.data.front().n_q_points_1d)); } else { @@ -468,7 +477,7 @@ namespace internal temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product; const unsigned int dofs_per_comp = (type == MatrixFreeFunctions::truncated_tensor) ? - Utilities::fixed_power(shape_info.fe_degree + 1) : + Utilities::fixed_power(shape_info.data.front().fe_degree + 1) : shape_info.dofs_per_component_on_cell; // expand dof_values to tensor product for truncated tensor products Number *values_dofs = @@ -612,7 +621,8 @@ namespace internal { values_dofs -= dofs_per_comp * n_components; unsigned int count_p = 0, count_q = 0; - const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree; + const int degree = + fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i) { for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j) @@ -630,7 +640,8 @@ namespace internal count_q += i * (degree + 1); } AssertDimension(count_q, - Utilities::fixed_power(shape_info.fe_degree + 1)); + Utilities::fixed_power( + shape_info.data.front().fe_degree + 1)); } } @@ -1039,8 +1050,9 @@ namespace internal const bool evaluate_gradients, const bool evaluate_hessians) { - AssertDimension(shape_info.shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); + AssertDimension( + shape_info.data.front().shape_gradients_collocation_eo.size(), + (fe_degree + 2) / 2 * (fe_degree + 1)); EvaluatorTensorProduct eval(AlignedVector(), - shape_info.shape_gradients_collocation_eo, - shape_info.shape_hessians_collocation_eo); + shape_info.data.front().shape_gradients_collocation_eo, + shape_info.data.front().shape_hessians_collocation_eo); constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim); for (unsigned int c = 0; c < n_components; c++) @@ -1115,8 +1127,9 @@ namespace internal const bool integrate_gradients, const bool add_into_values_array) { - AssertDimension(shape_info.shape_gradients_collocation_eo.size(), - (fe_degree + 2) / 2 * (fe_degree + 1)); + AssertDimension( + shape_info.data.front().shape_gradients_collocation_eo.size(), + (fe_degree + 2) / 2 * (fe_degree + 1)); EvaluatorTensorProduct eval(AlignedVector(), - shape_info.shape_gradients_collocation_eo, - shape_info.shape_hessians_collocation_eo); + shape_info.data.front().shape_gradients_collocation_eo, + shape_info.data.front().shape_hessians_collocation_eo); constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim); for (unsigned int c = 0; c < n_components; c++) @@ -1241,7 +1254,7 @@ namespace internal n_q_points_1d, 1, Number, - Number>::do_forward(shape_info.shape_values_eo, + Number>::do_forward(shape_info.data.front().shape_values_eo, values_dofs, values_quad); @@ -1292,8 +1305,9 @@ namespace internal "of lower degree, so the evaluation results would be " "wrong. Thus, this class does not permit the desired " "operation.")); - AssertDimension(shape_info.shape_gradients_collocation_eo.size(), - (n_q_points_1d + 1) / 2 * n_q_points_1d); + AssertDimension( + shape_info.data.front().shape_gradients_collocation_eo.size(), + (n_q_points_1d + 1) / 2 * n_q_points_1d); constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim); for (unsigned int c = 0; c < n_components; c++) @@ -1318,7 +1332,7 @@ namespace internal n_q_points_1d, 1, Number, - Number>::do_backward(shape_info.shape_values_eo, + Number>::do_backward(shape_info.data.front().shape_values_eo, add_into_values_array, values_quad, values_dofs); @@ -1360,29 +1374,29 @@ namespace internal { const AlignedVector &val1 = symmetric_evaluate ? - data.shape_values_eo : + data.data.front().shape_values_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_values : - data.values_within_subface[subface_index % 2]); + data.data.front().shape_values : + data.data.front().values_within_subface[subface_index % 2]); const AlignedVector &val2 = symmetric_evaluate ? - data.shape_values_eo : + data.data.front().shape_values_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_values : - data.values_within_subface[subface_index / 2]); + data.data.front().shape_values : + data.data.front().values_within_subface[subface_index / 2]); const AlignedVector &grad1 = symmetric_evaluate ? - data.shape_gradients_eo : + data.data.front().shape_gradients_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_gradients : - data.gradients_within_subface[subface_index % 2]); + data.data.front().shape_gradients : + data.data.front().gradients_within_subface[subface_index % 2]); const AlignedVector &grad2 = symmetric_evaluate ? - data.shape_gradients_eo : + data.data.front().shape_gradients_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_gradients : - data.gradients_within_subface[subface_index / 2]); + data.data.front().shape_gradients : + data.data.front().gradients_within_subface[subface_index / 2]); using Eval = internal::EvaluatorTensorProduct(), - data.fe_degree + 1, - data.n_q_points_1d); + data.data.front().fe_degree + 1, + data.data.front().n_q_points_1d); Eval eval2(val2, grad2, AlignedVector(), - data.fe_degree + 1, - data.n_q_points_1d); + data.data.front().fe_degree + 1, + data.data.front().n_q_points_1d); const unsigned int size_deg = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - (dim > 1 ? Utilities::fixed_power(data.fe_degree + 1) : 1); + (dim > 1 ? + Utilities::fixed_power(data.data.front().fe_degree + 1) : + 1); const unsigned int n_q_points = fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : @@ -1454,9 +1470,10 @@ namespace internal n_q_points_1d, n_q_points_1d, Number> - eval_grad(AlignedVector(), - data.shape_gradients_collocation_eo, - AlignedVector()); + eval_grad( + AlignedVector(), + data.data.front().shape_gradients_collocation_eo, + AlignedVector()); eval_grad.template gradients<0, true, false>( values_quad, gradients_quad); eval_grad.template gradients<1, true, false>( @@ -1521,29 +1538,29 @@ namespace internal { const AlignedVector &val1 = symmetric_evaluate ? - data.shape_values_eo : + data.data.front().shape_values_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_values : - data.values_within_subface[subface_index % 2]); + data.data.front().shape_values : + data.data.front().values_within_subface[subface_index % 2]); const AlignedVector &val2 = symmetric_evaluate ? - data.shape_values_eo : + data.data.front().shape_values_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_values : - data.values_within_subface[subface_index / 2]); + data.data.front().shape_values : + data.data.front().values_within_subface[subface_index / 2]); const AlignedVector &grad1 = symmetric_evaluate ? - data.shape_gradients_eo : + data.data.front().shape_gradients_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_gradients : - data.gradients_within_subface[subface_index % 2]); + data.data.front().shape_gradients : + data.data.front().gradients_within_subface[subface_index % 2]); const AlignedVector &grad2 = symmetric_evaluate ? - data.shape_gradients_eo : + data.data.front().shape_gradients_eo : (subface_index >= GeometryInfo::max_children_per_cell ? - data.shape_gradients : - data.gradients_within_subface[subface_index / 2]); + data.data.front().shape_gradients : + data.data.front().gradients_within_subface[subface_index / 2]); using Eval = internal::EvaluatorTensorProduct; - Eval eval1(val1, grad1, val1, data.fe_degree + 1, data.n_q_points_1d); - Eval eval2(val2, grad2, val1, data.fe_degree + 1, data.n_q_points_1d); + Eval eval1(val1, + grad1, + val1, + data.data.front().fe_degree + 1, + data.data.front().n_q_points_1d); + Eval eval2(val2, + grad2, + val1, + data.data.front().fe_degree + 1, + data.data.front().n_q_points_1d); const unsigned int size_deg = fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - (dim > 1 ? Utilities::fixed_power(data.fe_degree + 1) : 1); + (dim > 1 ? + Utilities::fixed_power(data.data.front().fe_degree + 1) : + 1); const unsigned int n_q_points = fe_degree > -1 ? Utilities::pow(n_q_points_1d, dim - 1) : @@ -1609,9 +1636,10 @@ namespace internal n_q_points_1d, n_q_points_1d, Number> - eval_grad(AlignedVector(), - data.shape_gradients_collocation_eo, - AlignedVector()); + eval_grad( + AlignedVector(), + data.data.front().shape_gradients_collocation_eo, + AlignedVector()); if (integrate_val) eval_grad.template gradients<1, false, true>( gradients_quad + n_q_points, values_quad); @@ -1687,10 +1715,10 @@ namespace internal fe_degree + 1, 0, Number> - evalf(data.shape_data_on_face[face_no % 2], + evalf(data.data.front().shape_data_on_face[face_no % 2], AlignedVector(), AlignedVector(), - data.fe_degree + 1, + data.data.front().fe_degree + 1, 0); const unsigned int in_stride = do_evaluate ? @@ -1865,8 +1893,9 @@ namespace internal fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(data.fe_degree + 1, dim - 1); + fe_degree > -1 ? + static_dofs_per_face : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); // we allocate small amounts of data on the stack to signal the compiler // that this temporary data is only needed for the calculations but the @@ -1972,8 +2001,9 @@ namespace internal fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(data.fe_degree + 1, dim - 1); + fe_degree > -1 ? + static_dofs_per_face : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); constexpr unsigned int stack_array_size_threshold = 100; @@ -2062,8 +2092,9 @@ namespace internal fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? static_dofs_per_face : - Utilities::pow(data.fe_degree + 1, dim - 1); + fe_degree > -1 ? + static_dofs_per_face : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); constexpr unsigned int stack_array_size_threshold = 100; @@ -2079,7 +2110,7 @@ namespace internal // case 1: contiguous and interleaved indices if (((evaluate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2103,7 +2134,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 1 + side]; + data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); const unsigned int *index_array = @@ -2154,7 +2185,7 @@ namespace internal // case 2: contiguous and interleaved indices with fixed stride else if (((evaluate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2174,7 +2205,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 1 + side]; + data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2240,7 +2271,7 @@ namespace internal // case 3: contiguous and interleaved indices with mixed stride else if (((evaluate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2270,7 +2301,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 1 + side]; + data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2379,7 +2410,7 @@ namespace internal // case 4: contiguous indices without interleaving else if (((evaluate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2400,7 +2431,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 1 + side]; + data.data.front().shape_data_on_face[0][fe_degree + 1 + side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2544,8 +2575,9 @@ namespace internal fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) : numbers::invalid_unsigned_int; const unsigned int dofs_per_face = - fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) : - Utilities::pow(data.fe_degree + 1, dim - 1); + fe_degree > -1 ? + Utilities::pow(fe_degree + 1, dim - 1) : + Utilities::pow(data.data.front().fe_degree + 1, dim - 1); constexpr unsigned int stack_array_size_threshold = 100; @@ -2597,7 +2629,7 @@ namespace internal // case 1: contiguous and interleaved indices if (((integrate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2621,7 +2653,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); const unsigned int *index_array = @@ -2675,7 +2707,7 @@ namespace internal // case 2: contiguous and interleaved indices with fixed stride else if (((integrate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2695,7 +2727,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2764,7 +2796,7 @@ namespace internal // case 3: contiguous and interleaved indices with mixed stride else if (((integrate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2794,7 +2826,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); @@ -2899,7 +2931,7 @@ namespace internal // case 4: contiguous indices without interleaving else if (((integrate_gradients == false && - data.nodal_at_cell_boundaries == true) || + data.data.front().nodal_at_cell_boundaries == true) || (data.element_type == internal::MatrixFreeFunctions::tensor_symmetric_hermite && fe_degree > 1)) && @@ -2921,7 +2953,7 @@ namespace internal // right (side==1) are the negative from the value at the left // (side==0), so we only read out one of them. const VectorizedArrayType grad_weight = - data.shape_data_on_face[0][fe_degree + 2 - side]; + data.data.front().shape_data_on_face[0][fe_degree + 2 - side]; AssertDimension(data.face_to_cell_index_hermite.size(1), 2 * dofs_per_face); const unsigned int *index_array = @@ -3062,9 +3094,10 @@ namespace internal fe_degree + 1, fe_degree + 1, Number> - evaluator(AlignedVector(), - AlignedVector(), - fe_eval.get_shape_info().inverse_shape_values_eo); + evaluator( + AlignedVector(), + AlignedVector(), + fe_eval.get_shape_info().data.front().inverse_shape_values_eo); for (unsigned int d = 0; d < n_components; ++d) { diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index 3f865caaf8..8f24381679 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -174,7 +174,7 @@ namespace internal const bool evaluate_gradients, const bool evaluate_hessians) { - const unsigned int runtime_degree = shape_info.fe_degree; + const unsigned int runtime_degree = shape_info.data.front().fe_degree; constexpr unsigned int start_n_q_points = degree + 1; if (runtime_degree == degree) Factory:: @@ -211,7 +211,7 @@ namespace internal const bool integrate_gradients, const bool sum_into_values_array = false) { - const int runtime_degree = shape_info.fe_degree; + const int runtime_degree = shape_info.data.front().fe_degree; constexpr unsigned int start_n_q_points = degree + 1; if (runtime_degree == degree) Factory:: @@ -276,7 +276,7 @@ namespace internal const bool evaluate_gradients, const bool evaluate_hessians) { - const int runtime_n_q_points_1d = shape_info.n_q_points_1d; + const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d; if (runtime_n_q_points_1d == n_q_points_1d) { if (n_q_points_1d == degree + 1 && @@ -349,7 +349,7 @@ namespace internal const bool integrate_gradients, const bool sum_into_values_array) { - const int runtime_n_q_points_1d = shape_info.n_q_points_1d; + const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d; if (runtime_n_q_points_1d == n_q_points_1d) { if (n_q_points_1d == degree + 1 && diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 8d4e14d739..6bb9519e87 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3371,7 +3371,7 @@ FEEvaluationBase:: Assert(scratch_data_array != nullptr, ExcInternalError()); const unsigned int tensor_dofs_per_component = - Utilities::fixed_power(this->data->fe_degree + 1); + Utilities::fixed_power(this->data->data.front().fe_degree + 1); const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; const unsigned int n_quadrature_points = @@ -6677,7 +6677,8 @@ FEEvaluation(fe_degree) != numbers::invalid_unsigned_int && - static_cast(fe_degree) != this->data->fe_degree) || + static_cast(fe_degree) != + this->data->data.front().fe_degree) || n_q_points != this->n_quadrature_points) { std::string message = @@ -6703,7 +6704,8 @@ FEEvaluation(fe_degree) == this->data->fe_degree) + if (static_cast(fe_degree) == + this->data->data.front().fe_degree) { proposed_dof_comp = dof_no; proposed_fe_comp = first_selected_component; @@ -6716,6 +6718,7 @@ FEEvaluationmatrix_info ->get_shape_info(no, 0, nf, this->active_fe_index, 0) + .data.front() .fe_degree == static_cast(fe_degree)) { proposed_dof_comp = no; @@ -6781,7 +6784,8 @@ FEEvaluationn_quadrature_points, 1. / dim)); message += "Wrong template arguments:\n"; message += " Did you mean FEEvaluationdata->fe_degree) + ","; + message += + Utilities::int_to_string(this->data->data.front().fe_degree) + ","; message += Utilities::int_to_string(proposed_n_q_points_1d); message += "," + Utilities::int_to_string(n_components); message += ",Number>(data"; @@ -6793,7 +6797,8 @@ FEEvaluationdata->fe_degree != static_cast(fe_degree)) + if (this->data->data.front().fe_degree != + static_cast(fe_degree)) correct_pos = " ^"; else correct_pos = " "; @@ -6803,7 +6808,8 @@ FEEvaluation(fe_degree) == this->data->fe_degree && + Assert(static_cast(fe_degree) == + this->data->data.front().fe_degree && n_q_points == this->n_quadrature_points, ExcMessage(message)); } @@ -6946,7 +6952,7 @@ FEEvaluationdata->n_q_points_1d : n_q_points_1d; + fe_degree == -1 ? this->data->data.front().n_q_points_1d : n_q_points_1d; // Cartesian mesh: not all quadrature points are stored, only the // diagonal. Hence, need to find the tensor product index and retrieve the diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 7cd7dd18c0..eec6794d26 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1704,6 +1704,7 @@ MatrixFree::initialize_indices( bool all_nodal = true; for (unsigned int c = 0; c < di.n_base_elements; ++c) if (!shape_info(di.global_base_element_offset + c, 0, 0, 0) + .data.front() .nodal_at_cell_boundaries) all_nodal = false; if (all_nodal == false) diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 5f93f4ed2f..e1d17a200f 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -1048,7 +1048,7 @@ namespace MatrixFreeOperators fe_degree, n_components, VectorizedArrayType>:: - apply(fe_eval.get_shape_info().inverse_shape_values_eo, + apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo, inverse_coefficients, n_actual_components, in_array, @@ -1077,7 +1077,7 @@ namespace MatrixFreeOperators n_components, VectorizedArrayType>:: transform_from_q_points_to_basis( - fe_eval.get_shape_info().inverse_shape_values_eo, + fe_eval.get_shape_info().data.front().inverse_shape_values_eo, n_actual_components, in_array, out_array); diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index afc38492be..d143e40441 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -94,46 +94,9 @@ namespace internal tensor_symmetric_plus_dg0 = 5 }; - /** - * The class that stores the shape functions, gradients and Hessians - * evaluated for a tensor product finite element and tensor product - * quadrature formula on the unit cell. Because of this structure, only - * one-dimensional data is stored. - * - * @ingroup matrixfree - * - * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 - */ template - struct ShapeInfo + struct UnivariateShapeData { - /** - * Empty constructor. Does nothing. - */ - ShapeInfo(); - - /** - * Constructor that initializes the data fields using the reinit method. - */ - template - ShapeInfo(const Quadrature<1> & quad, - const FiniteElement &fe, - const unsigned int base_element = 0); - - /** - * Initializes the data fields. Takes a one-dimensional quadrature - * formula and a finite element as arguments and evaluates the shape - * functions, gradients and Hessians on the one-dimensional unit cell. - * This function assumes that the finite element is derived from a one- - * dimensional element by a tensor product and that the zeroth shape - * function in zero evaluates to one. - */ - template - void - reinit(const Quadrature<1> & quad, - const FiniteElement &fe_dim, - const unsigned int base_element = 0); - /** * Return the memory consumption of this class in bytes. */ @@ -145,7 +108,7 @@ namespace internal * will select the most efficient algorithm based on the given element * type. */ - ElementType element_type; + ElementType element_type = ElementType::tensor_general; /** * Stores the shape values of the 1D finite element evaluated on all 1D @@ -275,6 +238,88 @@ namespace internal */ Quadrature<1> quadrature; + /** + * Stores the degree of the element. + */ + unsigned int fe_degree = 0; + + /** + * Stores the number of quadrature points per dimension. + */ + unsigned int n_q_points_1d = 0; + + /** + * Indicates whether the basis functions are nodal in 0 and 1, i.e., the + * end points of the unit cell. + */ + bool nodal_at_cell_boundaries = false; + }; + + + + /** + * The class that stores the shape functions, gradients and Hessians + * evaluated for a tensor product finite element and tensor product + * quadrature formula on the unit cell. Because of this structure, only + * one-dimensional data is stored. + * + * @ingroup matrixfree + * + * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 + */ + template + struct ShapeInfo + { + /** + * Encodes the type of element detected at construction. FEEvaluation + * will select the most efficient algorithm based on the given element + * type. + */ + ElementType element_type; + + /** + * Empty constructor. Does nothing. + */ + ShapeInfo(); + + /** + * Constructor that initializes the data fields using the reinit method. + */ + template + ShapeInfo(const Quadrature<1> & quad, + const FiniteElement &fe, + const unsigned int base_element = 0); + + /** + * Initializes the data fields. Takes a one-dimensional quadrature + * formula and a finite element as arguments and evaluates the shape + * functions, gradients and Hessians on the one-dimensional unit cell. + * This function assumes that the finite element is derived from a one- + * dimensional element by a tensor product and that the zeroth shape + * function in zero evaluates to one. + */ + template + void + reinit(const Quadrature<1> & quad, + const FiniteElement &fe_dim, + const unsigned int base_element = 0); + + /** + * Return data of univariate shape functions which defines the + * dimension @p dimension of tensor product shape functions + * regarding vector component @p component of the underlying + * finite element. + */ + const UnivariateShapeData & + get_shape_data(const unsigned int dimension = 0, + const unsigned int component = 0) const; + + /** + * Return the memory consumption of this class in bytes. + */ + std::size_t + memory_consumption() const; + /** * Renumbering from deal.II's numbering of cell degrees of freedom to * lexicographic numbering used inside the FEEvaluation schemes of the @@ -285,14 +330,28 @@ namespace internal std::vector lexicographic_numbering; /** - * Stores the degree of the element. + * Stores data of univariate shape functions defining the + * underlying tensor product finite element. */ - unsigned int fe_degree; + std::vector> data; /** - * Stores the number of quadrature points per dimension. + * Grants access to univariate shape function data of given + * dimension and vector component. Rows identify dimensions and + * columns identify vector components. + */ + dealii::Table<2, UnivariateShapeData *> data_access; + + /** + * Stores the number of space dimensions. */ - unsigned int n_q_points_1d; + unsigned int n_dimensions; + + /** + * Stores the number of vector components of the underlying + * vector-valued finite element. + */ + unsigned int n_components; /** * Stores the number of quadrature points in @p dim dimensions for a @@ -316,12 +375,6 @@ namespace internal */ unsigned int dofs_per_component_on_face; - /** - * Indicates whether the basis functions are nodal in 0 and 1, i.e., the - * end points of the unit cell. - */ - bool nodal_at_cell_boundaries; - /** * For nodal basis functions with nodes located at the boundary of the * unit cell, face integrals that involve only the values of the shape @@ -405,7 +458,8 @@ namespace internal * also fill the shape_???_eo fields. */ bool - check_1d_shapes_symmetric(const unsigned int n_q_points_1d); + check_1d_shapes_symmetric( + UnivariateShapeData &univariate_shape_data); /** * Check whether symmetric 1D basis functions are such that the shape @@ -414,7 +468,8 @@ namespace internal * that save some operations in the evaluation. */ bool - check_1d_shapes_collocation() const; + check_1d_shapes_collocation( + const UnivariateShapeData &univariate_shape_data) const; }; @@ -427,17 +482,28 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) : element_type(tensor_general) - , fe_degree(0) - , n_q_points_1d(0) + , n_dimensions(0) + , n_components(0) , n_q_points(0) , dofs_per_component_on_cell(0) , n_q_points_face(0) , dofs_per_component_on_face(0) - , nodal_at_cell_boundaries(false) { reinit(quad, fe_in, base_element_number); } + template + inline const UnivariateShapeData & + ShapeInfo::get_shape_data(const unsigned int dimension, + const unsigned int component) const + { + AssertDimension(n_dimensions, data_access.size(0)); + AssertDimension(n_components, data_access.size(1)); + AssertIndexRange(dimension, n_dimensions); + AssertIndexRange(component, n_components); + return *(data_access(dimension, component)); + } + } // end of namespace MatrixFreeFunctions } // end of namespace internal diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index c89ba34f22..da01d69018 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -61,13 +61,12 @@ namespace internal template ShapeInfo::ShapeInfo() : element_type(tensor_general) - , fe_degree(numbers::invalid_unsigned_int) - , n_q_points_1d(0) + , n_dimensions(0) + , n_components(0) , n_q_points(0) , dofs_per_component_on_cell(0) , n_q_points_face(0) , dofs_per_component_on_face(0) - , nodal_at_cell_boundaries(false) {} @@ -79,15 +78,43 @@ namespace internal const FiniteElement &fe_in, const unsigned int base_element_number) { - this->quadrature = quad; const FiniteElement *fe = &fe_in.base_element(base_element_number); + n_dimensions = dim; + n_components = fe_in.n_components(); Assert(fe->n_components() == 1, ExcMessage("FEEvaluation only works for scalar finite elements.")); - fe_degree = fe->degree; - n_q_points_1d = quad.size(); - + // assuming isotropy of dimensions and components + data.resize(1); + UnivariateShapeData *univariate_shape_data = &data.front(); + data_access.reinit(n_dimensions, n_components); + data_access.fill(univariate_shape_data); + univariate_shape_data->quadrature = quad; + univariate_shape_data->fe_degree = fe->degree; + univariate_shape_data->n_q_points_1d = quad.size(); + + // grant write access to common univariate shape data + auto &shape_values = univariate_shape_data->shape_values; + auto &shape_gradients = univariate_shape_data->shape_gradients; + auto &shape_hessians = univariate_shape_data->shape_hessians; + auto &shape_gradients_collocation = + univariate_shape_data->shape_gradients_collocation; + auto &shape_hessians_collocation = + univariate_shape_data->shape_hessians_collocation; + auto &inverse_shape_values = univariate_shape_data->inverse_shape_values; + auto &shape_data_on_face = univariate_shape_data->shape_data_on_face; + auto &values_within_subface = + univariate_shape_data->values_within_subface; + auto &gradients_within_subface = + univariate_shape_data->gradients_within_subface; + auto &hessians_within_subface = + univariate_shape_data->hessians_within_subface; + auto &nodal_at_cell_boundaries = + univariate_shape_data->nodal_at_cell_boundaries; + + const unsigned int fe_degree = fe->degree; + const unsigned int n_q_points_1d = quad.size(); const unsigned int n_dofs_1d = std::min(fe->dofs_per_cell, fe_degree + 1); // renumber (this is necessary for FE_Q, for example, since there the @@ -205,18 +232,18 @@ namespace internal dim > 1 ? Utilities::fixed_power(fe_degree + 1) : 1; const unsigned int array_size = n_dofs_1d * n_q_points_1d; - this->shape_gradients.resize_fast(array_size); - this->shape_values.resize_fast(array_size); - this->shape_hessians.resize_fast(array_size); - - this->shape_data_on_face[0].resize(3 * n_dofs_1d); - this->shape_data_on_face[1].resize(3 * n_dofs_1d); - this->values_within_subface[0].resize(array_size); - this->values_within_subface[1].resize(array_size); - this->gradients_within_subface[0].resize(array_size); - this->gradients_within_subface[1].resize(array_size); - this->hessians_within_subface[0].resize(array_size); - this->hessians_within_subface[1].resize(array_size); + shape_gradients.resize_fast(array_size); + shape_values.resize_fast(array_size); + shape_hessians.resize_fast(array_size); + + shape_data_on_face[0].resize(3 * n_dofs_1d); + shape_data_on_face[1].resize(3 * n_dofs_1d); + values_within_subface[0].resize(array_size); + values_within_subface[1].resize(array_size); + gradients_within_subface[0].resize(array_size); + gradients_within_subface[1].resize(array_size); + hessians_within_subface[0].resize(array_size); + hessians_within_subface[1].resize(array_size); for (unsigned int i = 0; i < n_dofs_1d; ++i) { @@ -255,18 +282,18 @@ namespace internal } // evaluate basis functions on the 1D faces, i.e., in zero and one - Point q_point = unit_point; - q_point[0] = 0; - this->shape_data_on_face[0][i] = fe->shape_value(my_i, q_point); - this->shape_data_on_face[0][i + n_dofs_1d] = + Point q_point = unit_point; + q_point[0] = 0; + shape_data_on_face[0][i] = fe->shape_value(my_i, q_point); + shape_data_on_face[0][i + n_dofs_1d] = fe->shape_grad(my_i, q_point)[0]; - this->shape_data_on_face[0][i + 2 * n_dofs_1d] = + shape_data_on_face[0][i + 2 * n_dofs_1d] = fe->shape_grad_grad(my_i, q_point)[0][0]; - q_point[0] = 1; - this->shape_data_on_face[1][i] = fe->shape_value(my_i, q_point); - this->shape_data_on_face[1][i + n_dofs_1d] = + q_point[0] = 1; + shape_data_on_face[1][i] = fe->shape_value(my_i, q_point); + shape_data_on_face[1][i + n_dofs_1d] = fe->shape_grad(my_i, q_point)[0]; - this->shape_data_on_face[1][i + 2 * n_dofs_1d] = + shape_data_on_face[1][i + 2 * n_dofs_1d] = fe->shape_grad_grad(my_i, q_point)[0][0]; } @@ -391,9 +418,9 @@ namespace internal } if (element_type == tensor_general && - check_1d_shapes_symmetric(n_q_points_1d)) + check_1d_shapes_symmetric(*univariate_shape_data)) { - if (check_1d_shapes_collocation()) + if (check_1d_shapes_collocation(*univariate_shape_data)) element_type = tensor_symmetric_collocation; else element_type = tensor_symmetric; @@ -412,7 +439,7 @@ namespace internal } } else if (element_type == tensor_symmetric_plus_dg0) - check_1d_shapes_symmetric(n_q_points_1d); + check_1d_shapes_symmetric(*univariate_shape_data); nodal_at_cell_boundaries = true; for (unsigned int i = 1; i < n_dofs_1d; ++i) @@ -494,6 +521,9 @@ namespace internal } } } + + // TODO !!! + univariate_shape_data->element_type = this->element_type; } @@ -501,11 +531,31 @@ namespace internal template bool ShapeInfo::check_1d_shapes_symmetric( - const unsigned int n_q_points_1d) + UnivariateShapeData &univariate_shape_data) { if (dofs_per_component_on_cell == 0) return false; + const auto n_q_points_1d = univariate_shape_data.n_q_points_1d; + const auto fe_degree = univariate_shape_data.fe_degree; + auto & shape_values = univariate_shape_data.shape_values; + auto & shape_gradients = univariate_shape_data.shape_gradients; + auto & shape_hessians = univariate_shape_data.shape_hessians; + auto & shape_gradients_collocation = + univariate_shape_data.shape_gradients_collocation; + auto &shape_hessians_collocation = + univariate_shape_data.shape_hessians_collocation; + auto &shape_values_eo = univariate_shape_data.shape_values_eo; + auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo; + auto &shape_hessians_eo = univariate_shape_data.shape_hessians_eo; + auto &shape_gradients_collocation_eo = + univariate_shape_data.shape_gradients_collocation_eo; + auto &shape_hessians_collocation_eo = + univariate_shape_data.shape_hessians_collocation_eo; + auto &inverse_shape_values = univariate_shape_data.inverse_shape_values; + auto &inverse_shape_values_eo = + univariate_shape_data.inverse_shape_values_eo; + const double zero_tol = std::is_same::value == true ? 1e-12 : 1e-7; // symmetry for values @@ -616,11 +666,15 @@ namespace internal template bool - ShapeInfo::check_1d_shapes_collocation() const + ShapeInfo::check_1d_shapes_collocation( + const UnivariateShapeData &univariate_shape_data) const { if (dofs_per_component_on_cell != n_q_points) return false; + const auto fe_degree = univariate_shape_data.fe_degree; + auto & shape_values = univariate_shape_data.shape_values; + const double zero_tol = std::is_same::value == true ? 1e-12 : 1e-7; // check: identity operation for shape values @@ -648,6 +702,16 @@ namespace internal template std::size_t ShapeInfo::memory_consumption() const + { + std::size_t memory = sizeof(*this); + for (const auto &univariate_shape_data : data) + memory += univariate_shape_data.memory_consumption(); + return memory; + } + + template + std::size_t + UnivariateShapeData::memory_consumption() const { std::size_t memory = sizeof(*this); memory += MemoryConsumption::memory_consumption(shape_values); diff --git a/tests/matrix_free/shape_info_inverse.cc b/tests/matrix_free/shape_info_inverse.cc index be98ff8a3a..34f62acfc6 100644 --- a/tests/matrix_free/shape_info_inverse.cc +++ b/tests/matrix_free/shape_info_inverse.cc @@ -43,19 +43,23 @@ test(const FiniteElement &fe, deallog << "Testing " << fe.get_name() << " with " << quadrature_name << "(" << quad.size() << ")" << std::endl; deallog << "shape values: " << std::endl; - for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i) + const auto &univariate_shape_data = shape_info.get_shape_data(0, 0); + for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i) { for (unsigned int q = 0; q < quad.size(); ++q) - deallog << std::setw(15) << shape_info.shape_values[i * quad.size() + q] + deallog << std::setw(15) + << univariate_shape_data.shape_values[i * quad.size() + q] << " "; deallog << std::endl; } deallog << "inverse shape values: " << std::endl; - for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i) + for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i) { for (unsigned int q = 0; q < quad.size(); ++q) - deallog << std::setw(15) - << shape_info.inverse_shape_values[i * quad.size() + q] << " "; + deallog + << std::setw(15) + << univariate_shape_data.inverse_shape_values[i * quad.size() + q] + << " "; deallog << std::endl; } deallog << "inverse shapes' * shapes: " << std::endl; @@ -64,22 +68,24 @@ test(const FiniteElement &fe, for (unsigned int j = 0; j < quad.size(); ++j) { double sum = 0; - for (unsigned int k = 0; k < shape_info.fe_degree + 1; ++k) - sum += shape_info.inverse_shape_values[k * quad.size() + i] * - shape_info.shape_values[k * quad.size() + j]; + for (unsigned int k = 0; k < univariate_shape_data.fe_degree + 1; ++k) + sum += + univariate_shape_data.inverse_shape_values[k * quad.size() + i] * + univariate_shape_data.shape_values[k * quad.size() + j]; deallog << std::setw(15) << sum << " "; } deallog << std::endl; } deallog << "inverse shapes * shapes': " << std::endl; - for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i) + for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i) { - for (unsigned int j = 0; j < shape_info.fe_degree + 1; ++j) + for (unsigned int j = 0; j < univariate_shape_data.fe_degree + 1; ++j) { double sum = 0; for (unsigned int k = 0; k < quad.size(); ++k) - sum += shape_info.inverse_shape_values[i * quad.size() + k] * - shape_info.shape_values[j * quad.size() + k]; + sum += + univariate_shape_data.inverse_shape_values[i * quad.size() + k] * + univariate_shape_data.shape_values[j * quad.size() + k]; deallog << std::setw(15) << sum << " "; } deallog << std::endl;