Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

FEPointEvaluation: Remove unit_points member, access them through MappingInfo instead #14926

Merged
merged 2 commits into from
Mar 25, 2023
Merged

FEPointEvaluation: Remove unit_points member, access them through MappingInfo instead #14926

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
df9e9e4
Remove unit_points member, access them through MappingInfo instead
bergbauer Mar 21, 2023
8ce64c0
Use reinit() function internal
bergbauer Mar 21, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
104 changes: 54 additions & 50 deletions include/deal.II/matrix_free/fe_point_evaluation.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -759,14 +759,14 @@ class FEPointEvaluation
unsigned int current_face_number;

/**
* The reference points specified at reinit().
* Bool indicating if fast path is chosen.
*/
ArrayView<const Point<dim>> 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 ----------------------
Expand All @@ -788,6 +788,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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);
}
Expand All @@ -806,6 +807,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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);
}
Expand Down Expand Up @@ -904,14 +906,14 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
fe_values->reinit(cell);
}

this->unit_points = unit_points;

const_cast<unsigned int &>(n_q_points) = unit_points.size();

if (update_flags & update_values)
values.resize(n_q_points, numbers::signaling_nan<value_type>());
if (update_flags & update_gradients)
gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());

is_reinitialized = true;
}


Expand All @@ -927,6 +929,13 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
const_cast<unsigned int &>(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<value_type>());
if (update_flags & update_gradients)
gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());

is_reinitialized = true;
}


Expand All @@ -943,6 +952,13 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::reinit(
const_cast<unsigned int &>(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<value_type>());
if (update_flags & update_gradients)
gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());

is_reinitialized = true;
}


Expand All @@ -953,20 +969,10 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
const ArrayView<const Number> & solution_values,
const EvaluationFlags::EvaluationFlags &evaluation_flag)
{
const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
if (precomputed_mapping)
{
unit_points =
mapping_info->get_unit_points(current_cell_index, current_face_number);

if (update_flags & update_values)
values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
if (update_flags & update_gradients)
gradients.resize(unit_points.size(),
numbers::signaling_nan<gradient_type>());
}
if (!is_reinitialized)
reinit(numbers::invalid_unsigned_int);
Comment on lines -956 to +973
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gassmoeller Could this be the problematic change which leads to #14984 ? Could you try to call FEPointEvaluation::reinit(numbers::invalid_unsigned_int) manually after NonMatching::MappingInfo::reinit()?


if (unit_points.size() == 0)
if (n_q_points == 0)
return;

AssertDimension(solution_values.size(), fe->dofs_per_cell);
Expand All @@ -989,9 +995,12 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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<gradient_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<Number>::size();
for (unsigned int i = 0; i < n_points; i += n_lanes)
Expand Down Expand Up @@ -1049,20 +1058,20 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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)
{
const Number value = solution_values[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)
internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components, Number>::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<dim, n_components, Number>::access(
values[q], d) +=
Expand All @@ -1072,20 +1081,20 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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)
{
const Number value = solution_values[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)
internal::FEPointEvaluation::
EvaluatorTypeTraits<dim, n_components, Number>::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<dim, n_components, Number>::access(
gradients[q], d) +=
Expand All @@ -1103,20 +1112,10 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
const ArrayView<Number> & solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags)
{
const bool precomputed_mapping = mapping_info_on_the_fly.get() == nullptr;
if (precomputed_mapping)
{
unit_points =
mapping_info->get_unit_points(current_cell_index, current_face_number);
if (!is_reinitialized)
reinit(numbers::invalid_unsigned_int);

if (update_flags & update_values)
values.resize(unit_points.size(), numbers::signaling_nan<value_type>());
if (update_flags & update_gradients)
gradients.resize(unit_points.size(),
numbers::signaling_nan<gradient_type>());
}

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;
Expand All @@ -1130,9 +1129,9 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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);
Expand All @@ -1143,6 +1142,9 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
n_components,
VectorizedArray<Number>>::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<Number>::size();
for (unsigned int i = 0; i < n_points; i += n_lanes)
Expand Down Expand Up @@ -1219,20 +1221,20 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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<dim, n_components, Number>::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::
Expand All @@ -1243,20 +1245,20 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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<dim, n_components, Number>::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::
Expand Down Expand Up @@ -1315,7 +1317,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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;
}

Expand All @@ -1327,7 +1329,7 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::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;
}

Expand Down Expand Up @@ -1402,7 +1404,9 @@ inline Point<dim>
FEPointEvaluation<n_components, dim, spacedim, Number>::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];
}

Expand Down