From df9e9e4a78c7a1bc26eba3a9f421ff4dac23df02 Mon Sep 17 00:00:00 2001 From: Maximilian Bergbauer Date: Tue, 21 Mar 2023 10:49:55 +0100 Subject: [PATCH] Remove unit_points member, access them through MappingInfo instead --- .../deal.II/matrix_free/fe_point_evaluation.h | 104 +++++++++++------- 1 file changed, 62 insertions(+), 42 deletions(-) diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h index 76e714638d..1f55b1487a 100644 --- a/include/deal.II/matrix_free/fe_point_evaluation.h +++ b/include/deal.II/matrix_free/fe_point_evaluation.h @@ -759,14 +759,14 @@ private: unsigned int current_face_number; /** - * The reference points specified at reinit(). + * Bool indicating if fast path is chosen. */ - ArrayView> unit_points; + bool fast_path; /** - * Bool indicating if fast path is chosen. + * Bool indicating if class is reinitialized and data vectors a resized. */ - bool fast_path; + bool is_reinitialized; }; // ----------------------- template and inline function ---------------------- @@ -788,6 +788,7 @@ FEPointEvaluation::FEPointEvaluation( , mapping_info(mapping_info_on_the_fly.get()) , current_cell_index(numbers::invalid_unsigned_int) , current_face_number(numbers::invalid_unsigned_int) + , is_reinitialized(false) { setup(first_selected_component); } @@ -806,6 +807,7 @@ FEPointEvaluation::FEPointEvaluation( , mapping_info(&mapping_info) , current_cell_index(numbers::invalid_unsigned_int) , current_face_number(numbers::invalid_unsigned_int) + , is_reinitialized(false) { setup(first_selected_component); } @@ -904,14 +906,14 @@ FEPointEvaluation::reinit( fe_values->reinit(cell); } - this->unit_points = unit_points; - const_cast(n_q_points) = unit_points.size(); if (update_flags & update_values) values.resize(n_q_points, numbers::signaling_nan()); if (update_flags & update_gradients) gradients.resize(n_q_points, numbers::signaling_nan()); + + is_reinitialized = true; } @@ -927,6 +929,13 @@ FEPointEvaluation::reinit( const_cast(n_q_points) = mapping_info->get_unit_points(current_cell_index, current_face_number) .size(); + + if (update_flags & update_values) + values.resize(n_q_points, numbers::signaling_nan()); + if (update_flags & update_gradients) + gradients.resize(n_q_points, numbers::signaling_nan()); + + is_reinitialized = true; } @@ -943,6 +952,13 @@ FEPointEvaluation::reinit( const_cast(n_q_points) = mapping_info->get_unit_points(current_cell_index, current_face_number) .size(); + + if (update_flags & update_values) + values.resize(n_q_points, numbers::signaling_nan()); + if (update_flags & update_gradients) + gradients.resize(n_q_points, numbers::signaling_nan()); + + is_reinitialized = true; } @@ -953,20 +969,18 @@ FEPointEvaluation::evaluate( const ArrayView & solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flag) { - const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr; - if (precomputed_mapping) + if (!is_reinitialized) { - unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - + const_cast(n_q_points) = + mapping_info->get_unit_points(current_cell_index, current_face_number) + .size(); if (update_flags & update_values) - values.resize(unit_points.size(), numbers::signaling_nan()); + values.resize(n_q_points, numbers::signaling_nan()); if (update_flags & update_gradients) - gradients.resize(unit_points.size(), - numbers::signaling_nan()); + gradients.resize(n_q_points, numbers::signaling_nan()); } - if (unit_points.size() == 0) + if (n_q_points == 0) return; AssertDimension(solution_values.size(), fe->dofs_per_cell); @@ -989,9 +1003,12 @@ FEPointEvaluation::evaluate( // unit gradients are currently only implemented with the fast tensor // path - unit_gradients.resize(unit_points.size(), + unit_gradients.resize(n_q_points, numbers::signaling_nan()); + const auto unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); + const std::size_t n_points = unit_points.size(); const std::size_t n_lanes = VectorizedArray::size(); for (unsigned int i = 0; i < n_points; i += n_lanes) @@ -1049,7 +1066,7 @@ FEPointEvaluation::evaluate( if (evaluation_flag & EvaluationFlags::values) { - values.resize(unit_points.size()); + values.resize(n_q_points); std::fill(values.begin(), values.end(), value_type()); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { @@ -1057,12 +1074,12 @@ FEPointEvaluation::evaluate( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) internal::FEPointEvaluation:: EvaluatorTypeTraits::access( values[q], d) += fe_values->shape_value(i, q) * value; else if (nonzero_shape_function_component[i][d]) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) internal::FEPointEvaluation:: EvaluatorTypeTraits::access( values[q], d) += @@ -1072,7 +1089,7 @@ FEPointEvaluation::evaluate( if (evaluation_flag & EvaluationFlags::gradients) { - gradients.resize(unit_points.size()); + gradients.resize(n_q_points); std::fill(gradients.begin(), gradients.end(), gradient_type()); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { @@ -1080,12 +1097,12 @@ FEPointEvaluation::evaluate( for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) internal::FEPointEvaluation:: EvaluatorTypeTraits::access( gradients[q], d) += fe_values->shape_grad(i, q) * value; else if (nonzero_shape_function_component[i][d]) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) internal::FEPointEvaluation:: EvaluatorTypeTraits::access( gradients[q], d) += @@ -1103,20 +1120,18 @@ FEPointEvaluation::integrate( const ArrayView & solution_values, const EvaluationFlags::EvaluationFlags &integration_flags) { - const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr; - if (precomputed_mapping) + if (!is_reinitialized) { - unit_points = - mapping_info->get_unit_points(current_cell_index, current_face_number); - + const_cast(n_q_points) = + mapping_info->get_unit_points(current_cell_index, current_face_number) + .size(); if (update_flags & update_values) - values.resize(unit_points.size(), numbers::signaling_nan()); + values.resize(n_q_points, numbers::signaling_nan()); if (update_flags & update_gradients) - gradients.resize(unit_points.size(), - numbers::signaling_nan()); + gradients.resize(n_q_points, numbers::signaling_nan()); } - if (unit_points.size() == 0) // no evaluation points provided + if (n_q_points == 0) // no evaluation points provided { std::fill(solution_values.begin(), solution_values.end(), 0.0); return; @@ -1130,9 +1145,9 @@ FEPointEvaluation::integrate( // fast path with tensor product integration if (integration_flags & EvaluationFlags::values) - AssertIndexRange(unit_points.size(), values.size() + 1); + AssertIndexRange(n_q_points, values.size() + 1); if (integration_flags & EvaluationFlags::gradients) - AssertIndexRange(unit_points.size(), gradients.size() + 1); + AssertIndexRange(n_q_points, gradients.size() + 1); if (solution_renumbered_vectorized.size() != dofs_per_component) solution_renumbered_vectorized.resize(dofs_per_component); @@ -1143,6 +1158,9 @@ FEPointEvaluation::integrate( n_components, VectorizedArray>::value_type()); + const auto unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); + const std::size_t n_points = unit_points.size(); const std::size_t n_lanes = VectorizedArray::size(); for (unsigned int i = 0; i < n_points; i += n_lanes) @@ -1219,20 +1237,20 @@ FEPointEvaluation::integrate( if (integration_flags & EvaluationFlags::values) { - AssertIndexRange(unit_points.size(), values.size() + 1); + AssertIndexRange(n_q_points, values.size() + 1); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) solution_values[i] += fe_values->shape_value(i, q) * internal::FEPointEvaluation:: EvaluatorTypeTraits::access( values[q], d); else if (nonzero_shape_function_component[i][d]) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) solution_values[i] += fe_values->shape_value_component(i, q, d) * internal::FEPointEvaluation:: @@ -1243,20 +1261,20 @@ FEPointEvaluation::integrate( if (integration_flags & EvaluationFlags::gradients) { - AssertIndexRange(unit_points.size(), gradients.size() + 1); + AssertIndexRange(n_q_points, gradients.size() + 1); for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i) { for (unsigned int d = 0; d < n_components; ++d) if (nonzero_shape_function_component[i][d] && (fe->is_primitive(i) || fe->is_primitive())) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) solution_values[i] += fe_values->shape_grad(i, q) * internal::FEPointEvaluation:: EvaluatorTypeTraits::access( gradients[q], d); else if (nonzero_shape_function_component[i][d]) - for (unsigned int q = 0; q < unit_points.size(); ++q) + for (unsigned int q = 0; q < n_q_points; ++q) solution_values[i] += fe_values->shape_grad_component(i, q, d) * internal::FEPointEvaluation:: @@ -1315,7 +1333,7 @@ FEPointEvaluation::submit_value( const value_type & value, const unsigned int point_index) { - AssertIndexRange(point_index, unit_points.size()); + AssertIndexRange(point_index, n_q_points); values[point_index] = value; } @@ -1327,7 +1345,7 @@ FEPointEvaluation::submit_gradient( const gradient_type &gradient, const unsigned int point_index) { - AssertIndexRange(point_index, unit_points.size()); + AssertIndexRange(point_index, n_q_points); gradients[point_index] = gradient; } @@ -1402,7 +1420,9 @@ inline Point FEPointEvaluation::unit_point( const unsigned int point_index) const { - AssertIndexRange(point_index, unit_points.size()); + AssertIndexRange(point_index, n_q_points); + const auto unit_points = + mapping_info->get_unit_points(current_cell_index, current_face_number); return unit_points[point_index]; } -- 2.39.5