From ea622a62f3633b0520620612422db41c96525a48 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 13 Feb 2026 18:34:01 -0500 Subject: [PATCH 1/4] handle out of bound mode --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 103 ++++++++++++++++---- src/pcms/field.h | 21 ++++ test/CMakeLists.txt | 3 +- test/test_omega_h_field2_outofbounds.cpp | 79 +++++++++++++++ 4 files changed, 185 insertions(+), 21 deletions(-) create mode 100644 test/test_omega_h_field2_outofbounds.cpp diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index ea45bf10..ed7f3960 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -158,49 +158,97 @@ struct OmegaHField2LocalizationHint { OmegaHField2LocalizationHint( Omega_h::Mesh& mesh, - Kokkos::View search_results) - : offsets_("", mesh.nelems() + 1), - coordinates_("", search_results.size(), mesh.dim() + 1), - indices_("", search_results.size()) + Kokkos::View search_results, + OutOfBoundsMode mode, + Real fill_value) + : mode_(mode), + fill_value_(fill_value), + num_valid_(0), + num_missing_(0) { - Kokkos::View elem_counts("", mesh.nelems()); - + // First pass: count valid and invalid points + std::vector valid_point_indices; + std::vector missing_point_indices; + for (size_t i = 0; i < search_results.size(); ++i) { auto [dim, elem_idx, coord] = search_results(i); - elem_counts[elem_idx] += 1; + bool is_valid = (static_cast(dim) == mesh.dim()) && + (elem_idx >= 0) && (elem_idx < mesh.nelems()); + + if (is_valid) { + valid_point_indices.push_back(i); + } else { + missing_point_indices.push_back(i); + } + } + + num_valid_ = valid_point_indices.size(); + num_missing_ = missing_point_indices.size(); + + // Handle missing points based on mode + if (num_missing_ > 0 && mode_ == OutOfBoundsMode::ERROR) { + PCMS_ALWAYS_ASSERT(false && "Points found outside mesh domain"); } + if (num_missing_ > 0 && mode_ == OutOfBoundsMode::CLAMP) { + PCMS_ALWAYS_ASSERT(false && "CLAMP mode not implemented yet"); + } + + // Allocate arrays for valid points only + offsets_ = Kokkos::View("offsets", mesh.nelems() + 1); + coordinates_ = Kokkos::View("coordinates", num_valid_, mesh.dim() + 1); + indices_ = Kokkos::View("indices", num_valid_); + + // Store missing point indices + if (num_missing_ > 0) { + missing_indices_ = Kokkos::View("missing_indices", num_missing_); + for (size_t i = 0; i < num_missing_; ++i) { + missing_indices_(i) = static_cast(missing_point_indices[i]); + } + } + + // Count points per element (valid points only) + Kokkos::View elem_counts("elem_counts", mesh.nelems()); + for (size_t i = 0; i < num_valid_; ++i) { + auto [dim, elem_idx, coord] = search_results(valid_point_indices[i]); + elem_counts[elem_idx] += 1; + } + + // Compute offsets LO total; - ComputeOffsetsFunctor functor(offsets_, elem_counts); Kokkos::parallel_scan( "ComputeOffsets", Kokkos::RangePolicy(0, mesh.nelems()), functor, total); offsets_(mesh.nelems()) = total; - - for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, coord] = search_results(i); - // currently don't handle case where point is on a boundary - PCMS_ALWAYS_ASSERT(static_cast(dim) == mesh.dim()); - // element should be inside the domain (positive) - PCMS_ALWAYS_ASSERT(elem_idx >= 0 && elem_idx < mesh.nelems()); + + // Fill coordinates and indices for valid points + for (size_t i = 0; i < num_valid_; ++i) { + size_t orig_idx = valid_point_indices[i]; + auto [dim, elem_idx, coord] = search_results(orig_idx); elem_counts(elem_idx) -= 1; LO index = offsets_(elem_idx) + elem_counts(elem_idx); for (int j = 0; j < (mesh.dim() + 1); ++j) { coordinates_(index, j) = coord[j]; } - // coordinates_(index, mesh.dim()) = coord[0]; - indices_(index) = i; + indices_(index) = static_cast(orig_idx); } } + OutOfBoundsMode mode_; + Real fill_value_; + size_t num_valid_; + size_t num_missing_; + // offsets is the number of points in each element Kokkos::View offsets_; // coordinates are the parametric coordinates of each point Kokkos::View coordinates_; - // indices are the index of the original point + // indices are the index of the original point (for valid points) Kokkos::View indices_; + // indices of points not found in mesh + Kokkos::View missing_indices_; }; /* @@ -210,7 +258,7 @@ OmegaHField2::OmegaHField2(const OmegaHFieldLayout& layout) : layout_(layout), mesh_(layout.GetMesh()), search_(mesh_, 10, 10), - dof_holder_data_("", static_cast(layout.OwnedSize())) + dof_holder_data_("dof_holder_data", static_cast(layout.OwnedSize())) { auto nodes_per_dim = layout.GetNodesPerDim(); if (nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) { @@ -302,7 +350,8 @@ LocalizationHint OmegaHField2::GetLocalizationHint( Kokkos::View results_h( "results_h", results.size()); Kokkos::deep_copy(results_h, results); - auto hint = std::make_shared(mesh_, results_h); + auto hint = std::make_shared( + mesh_, results_h, out_of_bounds_mode_, fill_value_); return LocalizationHint{hint}; } @@ -332,11 +381,24 @@ void OmegaHField2::Evaluate(LocalizationHint location, "eval_results_h", eval_results.extent(0), eval_results.extent(1)); deep_copy_mismatch_layouts(eval_results_h, eval_results); Rank1View values = results.GetValues(); + + // Copy results for valid points Kokkos::parallel_for( "CopyEvalResultsToValues", Kokkos::RangePolicy( 0, eval_results_h.extent(0)), KOKKOS_LAMBDA(LO i) { values[hint.indices_(i)] = eval_results_h(i, 0); }); + + // Handle missing points based on mode + if (hint.num_missing_ > 0 && hint.mode_ == OutOfBoundsMode::FILL) { + Kokkos::parallel_for( + "FillMissingValues", + Kokkos::RangePolicy( + 0, hint.num_missing_), + KOKKOS_LAMBDA(LO i) { + values[hint.missing_indices_(i)] = hint.fill_value_; + }); + } } void OmegaHField2::EvaluateGradient( @@ -388,4 +450,5 @@ void OmegaHField2::Deserialize( SetDOFHolderData(pcms::make_const_array_view(sorted_buffer)); } + } // namespace pcms diff --git a/src/pcms/field.h b/src/pcms/field.h index 15b40c18..efb59467 100644 --- a/src/pcms/field.h +++ b/src/pcms/field.h @@ -76,6 +76,13 @@ struct LocalizationHint std::shared_ptr data = nullptr; }; +enum class OutOfBoundsMode +{ + ERROR, // Throw error when points are out of bounds + FILL, // Fill out-of-bounds points with a fill value + CLAMP // Clamp to nearest boundary cell (extrapolate) +}; + class FieldLayout; /* @@ -127,7 +134,21 @@ class FieldT Rank1View buffer, Rank1View permutation) = 0; + // Out-of-bounds handling + void SetOutOfBoundsMode(OutOfBoundsMode mode, Real fill_value = 0.0) + { + out_of_bounds_mode_ = mode; + fill_value_ = fill_value; + } + + OutOfBoundsMode GetOutOfBoundsMode() const { return out_of_bounds_mode_; } + Real GetFillValue() const { return fill_value_; } + virtual ~FieldT() noexcept = default; + +protected: + OutOfBoundsMode out_of_bounds_mode_ = OutOfBoundsMode::ERROR; + Real fill_value_ = 0.0; }; // Should statically instantiate types using FieldPtr = diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4325dab2..9ebdd032 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -394,7 +394,8 @@ if(Catch2_FOUND) test_normalisation.cpp test_spline_interpolator.cpp test_svd_serial.cpp - test_interpolation_class.cpp) + test_interpolation_class.cpp + test_omega_h_field2_outofbounds.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) target_link_libraries(unit_tests PUBLIC Catch2::Catch2 pcms::core diff --git a/test/test_omega_h_field2_outofbounds.cpp b/test/test_omega_h_field2_outofbounds.cpp new file mode 100644 index 00000000..e4517ed0 --- /dev/null +++ b/test/test_omega_h_field2_outofbounds.cpp @@ -0,0 +1,79 @@ +#include +#include +#include +#include +#include +#include "pcms/adapter/omega_h/omega_h_field2.h" +#include "pcms/create_field.h" +#include +#include + +using pcms::Real; + +TEST_CASE("omega_h_field2 out of bounds FILL mode") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + // Create a 1x1 box mesh (coords from 0 to 1) + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false); + auto layout = + pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + const auto nverts = mesh.nents(0); + auto mesh_coords = mesh.coords(); + + // Set up a simple linear field + auto f = KOKKOS_LAMBDA(Real x, Real y) { return x + y; }; + Omega_h::Write test_f(nverts); + Omega_h::parallel_for( + nverts, OMEGA_H_LAMBDA(int i) { + Real x = mesh_coords[2 * i + 0]; + Real y = mesh_coords[2 * i + 1]; + test_f[i] = f(x, y); + }); + Omega_h::HostWrite test_f_host(test_f); + auto field = layout->CreateField(); + field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); + + // Set FILL mode with fill value of -999.0 + Real fill_value = -999.0; + field->SetOutOfBoundsMode(pcms::OutOfBoundsMode::FILL, fill_value); + + // Test points - mix of inside and outside + std::vector coords = { + 0.5, 0.5, // inside - should evaluate normally + 1.5, 0.5, // outside (x > 1) - should return fill_value + 0.5, -0.1, // outside (y < 0) - should return fill_value + 0.3, 0.7, // inside - should evaluate normally + -0.1, 0.5, // outside (x < 0) - should return fill_value + }; + + std::vector evaluation(coords.size() / 2); + pcms::Rank1View eval_view{evaluation.data(), + evaluation.size()}; + pcms::Rank2View coords_view( + coords.data(), coords.size() / 2, 2); + pcms::FieldDataView data_view( + eval_view, field->GetCoordinateSystem()); + pcms::CoordinateView coordinate_view{ + field->GetCoordinateSystem(), coords_view}; + + auto locale = field->GetLocalizationHint(coordinate_view); + field->Evaluate(locale, data_view); + + // Check results + // Point 0: (0.5, 0.5) - inside, should be close to f(0.5, 0.5) = 1.0 + REQUIRE(std::abs(evaluation[0] - 1.0) < 0.1); + + // Point 1: (1.5, 0.5) - outside, should be fill_value + REQUIRE(evaluation[1] == fill_value); + + // Point 2: (0.5, -0.1) - outside, should be fill_value + REQUIRE(evaluation[2] == fill_value); + + // Point 3: (0.3, 0.7) - inside, should be close to f(0.3, 0.7) = 1.0 + REQUIRE(std::abs(evaluation[3] - 1.0) < 0.1); + + // Point 4: (-0.1, 0.5) - outside, should be fill_value + REQUIRE(evaluation[4] == fill_value); +} \ No newline at end of file From 71984dc2c0e090f469c6b42f17df2008928a1918 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 13 Feb 2026 19:44:20 -0500 Subject: [PATCH 2/4] reformat --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 51 ++++++++++----------- src/pcms/field.h | 8 ++-- test/test_omega_h_field2_outofbounds.cpp | 25 +++++----- 3 files changed, 43 insertions(+), 41 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index ed7f3960..e1abcf5c 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -159,32 +159,28 @@ struct OmegaHField2LocalizationHint OmegaHField2LocalizationHint( Omega_h::Mesh& mesh, Kokkos::View search_results, - OutOfBoundsMode mode, - Real fill_value) - : mode_(mode), - fill_value_(fill_value), - num_valid_(0), - num_missing_(0) + OutOfBoundsMode mode, Real fill_value) + : mode_(mode), fill_value_(fill_value), num_valid_(0), num_missing_(0) { // First pass: count valid and invalid points std::vector valid_point_indices; std::vector missing_point_indices; - + for (size_t i = 0; i < search_results.size(); ++i) { auto [dim, elem_idx, coord] = search_results(i); - bool is_valid = (static_cast(dim) == mesh.dim()) && + bool is_valid = (static_cast(dim) == mesh.dim()) && (elem_idx >= 0) && (elem_idx < mesh.nelems()); - + if (is_valid) { valid_point_indices.push_back(i); } else { missing_point_indices.push_back(i); } } - + num_valid_ = valid_point_indices.size(); num_missing_ = missing_point_indices.size(); - + // Handle missing points based on mode if (num_missing_ > 0 && mode_ == OutOfBoundsMode::ERROR) { PCMS_ALWAYS_ASSERT(false && "Points found outside mesh domain"); @@ -193,27 +189,30 @@ struct OmegaHField2LocalizationHint if (num_missing_ > 0 && mode_ == OutOfBoundsMode::CLAMP) { PCMS_ALWAYS_ASSERT(false && "CLAMP mode not implemented yet"); } - + // Allocate arrays for valid points only offsets_ = Kokkos::View("offsets", mesh.nelems() + 1); - coordinates_ = Kokkos::View("coordinates", num_valid_, mesh.dim() + 1); + coordinates_ = Kokkos::View( + "coordinates", num_valid_, mesh.dim() + 1); indices_ = Kokkos::View("indices", num_valid_); - + // Store missing point indices if (num_missing_ > 0) { - missing_indices_ = Kokkos::View("missing_indices", num_missing_); + missing_indices_ = + Kokkos::View("missing_indices", num_missing_); for (size_t i = 0; i < num_missing_; ++i) { missing_indices_(i) = static_cast(missing_point_indices[i]); } } - + // Count points per element (valid points only) - Kokkos::View elem_counts("elem_counts", mesh.nelems()); + Kokkos::View elem_counts("elem_counts", + mesh.nelems()); for (size_t i = 0; i < num_valid_; ++i) { auto [dim, elem_idx, coord] = search_results(valid_point_indices[i]); elem_counts[elem_idx] += 1; } - + // Compute offsets LO total; ComputeOffsetsFunctor functor(offsets_, elem_counts); @@ -222,7 +221,7 @@ struct OmegaHField2LocalizationHint Kokkos::RangePolicy(0, mesh.nelems()), functor, total); offsets_(mesh.nelems()) = total; - + // Fill coordinates and indices for valid points for (size_t i = 0; i < num_valid_; ++i) { size_t orig_idx = valid_point_indices[i]; @@ -240,7 +239,7 @@ struct OmegaHField2LocalizationHint Real fill_value_; size_t num_valid_; size_t num_missing_; - + // offsets is the number of points in each element Kokkos::View offsets_; // coordinates are the parametric coordinates of each point @@ -381,22 +380,22 @@ void OmegaHField2::Evaluate(LocalizationHint location, "eval_results_h", eval_results.extent(0), eval_results.extent(1)); deep_copy_mismatch_layouts(eval_results_h, eval_results); Rank1View values = results.GetValues(); - + // Copy results for valid points Kokkos::parallel_for( "CopyEvalResultsToValues", Kokkos::RangePolicy( 0, eval_results_h.extent(0)), KOKKOS_LAMBDA(LO i) { values[hint.indices_(i)] = eval_results_h(i, 0); }); - + // Handle missing points based on mode if (hint.num_missing_ > 0 && hint.mode_ == OutOfBoundsMode::FILL) { Kokkos::parallel_for( "FillMissingValues", - Kokkos::RangePolicy( - 0, hint.num_missing_), - KOKKOS_LAMBDA(LO i) { - values[hint.missing_indices_(i)] = hint.fill_value_; + Kokkos::RangePolicy(0, + hint.num_missing_), + KOKKOS_LAMBDA(LO i) { + values[hint.missing_indices_(i)] = hint.fill_value_; }); } } diff --git a/src/pcms/field.h b/src/pcms/field.h index efb59467..f4ee8951 100644 --- a/src/pcms/field.h +++ b/src/pcms/field.h @@ -78,9 +78,9 @@ struct LocalizationHint enum class OutOfBoundsMode { - ERROR, // Throw error when points are out of bounds - FILL, // Fill out-of-bounds points with a fill value - CLAMP // Clamp to nearest boundary cell (extrapolate) + ERROR, // Throw error when points are out of bounds + FILL, // Fill out-of-bounds points with a fill value + CLAMP // Clamp to nearest boundary cell (extrapolate) }; class FieldLayout; @@ -140,7 +140,7 @@ class FieldT out_of_bounds_mode_ = mode; fill_value_ = fill_value; } - + OutOfBoundsMode GetOutOfBoundsMode() const { return out_of_bounds_mode_; } Real GetFillValue() const { return fill_value_; } diff --git a/test/test_omega_h_field2_outofbounds.cpp b/test/test_omega_h_field2_outofbounds.cpp index e4517ed0..ec7b399a 100644 --- a/test/test_omega_h_field2_outofbounds.cpp +++ b/test/test_omega_h_field2_outofbounds.cpp @@ -21,9 +21,12 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); const auto nverts = mesh.nents(0); auto mesh_coords = mesh.coords(); - + // Set up a simple linear field - auto f = KOKKOS_LAMBDA(Real x, Real y) { return x + y; }; + auto f = KOKKOS_LAMBDA(Real x, Real y) + { + return x + y; + }; Omega_h::Write test_f(nverts); Omega_h::parallel_for( nverts, OMEGA_H_LAMBDA(int i) { @@ -41,11 +44,11 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") // Test points - mix of inside and outside std::vector coords = { - 0.5, 0.5, // inside - should evaluate normally - 1.5, 0.5, // outside (x > 1) - should return fill_value - 0.5, -0.1, // outside (y < 0) - should return fill_value - 0.3, 0.7, // inside - should evaluate normally - -0.1, 0.5, // outside (x < 0) - should return fill_value + 0.5, 0.5, // inside - should evaluate normally + 1.5, 0.5, // outside (x > 1) - should return fill_value + 0.5, -0.1, // outside (y < 0) - should return fill_value + 0.3, 0.7, // inside - should evaluate normally + -0.1, 0.5, // outside (x < 0) - should return fill_value }; std::vector evaluation(coords.size() / 2); @@ -64,16 +67,16 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") // Check results // Point 0: (0.5, 0.5) - inside, should be close to f(0.5, 0.5) = 1.0 REQUIRE(std::abs(evaluation[0] - 1.0) < 0.1); - + // Point 1: (1.5, 0.5) - outside, should be fill_value REQUIRE(evaluation[1] == fill_value); - + // Point 2: (0.5, -0.1) - outside, should be fill_value REQUIRE(evaluation[2] == fill_value); - + // Point 3: (0.3, 0.7) - inside, should be close to f(0.3, 0.7) = 1.0 REQUIRE(std::abs(evaluation[3] - 1.0) < 0.1); - + // Point 4: (-0.1, 0.5) - outside, should be fill_value REQUIRE(evaluation[4] == fill_value); } \ No newline at end of file From ef55cbcd67b759b70c47f1adc88ed59109d6fd01 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Tue, 17 Feb 2026 13:51:55 -0500 Subject: [PATCH 3/4] enum NEAREST_BOUNDARY and detect missing --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 35 ++++++++++++--------- src/pcms/field.h | 6 ++-- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index e1abcf5c..e4e03d10 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -166,15 +166,26 @@ struct OmegaHField2LocalizationHint std::vector valid_point_indices; std::vector missing_point_indices; - for (size_t i = 0; i < search_results.size(); ++i) { - auto [dim, elem_idx, coord] = search_results(i); - bool is_valid = (static_cast(dim) == mesh.dim()) && - (elem_idx >= 0) && (elem_idx < mesh.nelems()); - - if (is_valid) { + if (mode_ == OutOfBoundsMode::ERROR) { + // Error mode - throw error immediately if any point is out of bounds + for (size_t i = 0; i < search_results.size(); ++i) { + auto [dim, elem_idx, coord] = search_results(i); + bool is_missing = + (static_cast(dim) != mesh.dim()) || (elem_idx < 0); + PCMS_ALWAYS_ASSERT(!is_missing && "Points found outside mesh domain"); valid_point_indices.push_back(i); - } else { - missing_point_indices.push_back(i); + } + } else { + // Other modes - collect valid and missing points separately + for (size_t i = 0; i < search_results.size(); ++i) { + auto [dim, elem_idx, coord] = search_results(i); + bool is_missing = + (static_cast(dim) != mesh.dim()) || (elem_idx < 0); + if (is_missing) { + missing_point_indices.push_back(i); + } else { + valid_point_indices.push_back(i); + } } } @@ -182,12 +193,8 @@ struct OmegaHField2LocalizationHint num_missing_ = missing_point_indices.size(); // Handle missing points based on mode - if (num_missing_ > 0 && mode_ == OutOfBoundsMode::ERROR) { - PCMS_ALWAYS_ASSERT(false && "Points found outside mesh domain"); - } - - if (num_missing_ > 0 && mode_ == OutOfBoundsMode::CLAMP) { - PCMS_ALWAYS_ASSERT(false && "CLAMP mode not implemented yet"); + if (num_missing_ > 0 && mode_ == OutOfBoundsMode::NEAREST_BOUNDARY) { + PCMS_ALWAYS_ASSERT(false && "NEAREST_BOUNDARY mode not implemented yet"); } // Allocate arrays for valid points only diff --git a/src/pcms/field.h b/src/pcms/field.h index f4ee8951..93389ca7 100644 --- a/src/pcms/field.h +++ b/src/pcms/field.h @@ -78,9 +78,9 @@ struct LocalizationHint enum class OutOfBoundsMode { - ERROR, // Throw error when points are out of bounds - FILL, // Fill out-of-bounds points with a fill value - CLAMP // Clamp to nearest boundary cell (extrapolate) + ERROR, // Throw error when points are out of bounds + FILL, // Fill out-of-bounds points with a fill value + NEAREST_BOUNDARY // Map to nearest boundary cell (extrapolate) }; class FieldLayout; From 0193d8ca97b7faee85dab99824544537049e6897 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 20 Feb 2026 17:54:34 -0500 Subject: [PATCH 4/4] remove fill_value from localization --- src/pcms/adapter/omega_h/omega_h_field2.cpp | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/pcms/adapter/omega_h/omega_h_field2.cpp b/src/pcms/adapter/omega_h/omega_h_field2.cpp index e4e03d10..e63858cd 100644 --- a/src/pcms/adapter/omega_h/omega_h_field2.cpp +++ b/src/pcms/adapter/omega_h/omega_h_field2.cpp @@ -159,8 +159,8 @@ struct OmegaHField2LocalizationHint OmegaHField2LocalizationHint( Omega_h::Mesh& mesh, Kokkos::View search_results, - OutOfBoundsMode mode, Real fill_value) - : mode_(mode), fill_value_(fill_value), num_valid_(0), num_missing_(0) + OutOfBoundsMode mode) + : mode_(mode), num_valid_(0), num_missing_(0) { // First pass: count valid and invalid points std::vector valid_point_indices; @@ -243,7 +243,6 @@ struct OmegaHField2LocalizationHint } OutOfBoundsMode mode_; - Real fill_value_; size_t num_valid_; size_t num_missing_; @@ -357,7 +356,7 @@ LocalizationHint OmegaHField2::GetLocalizationHint( "results_h", results.size()); Kokkos::deep_copy(results_h, results); auto hint = std::make_shared( - mesh_, results_h, out_of_bounds_mode_, fill_value_); + mesh_, results_h, out_of_bounds_mode_); return LocalizationHint{hint}; } @@ -397,13 +396,12 @@ void OmegaHField2::Evaluate(LocalizationHint location, // Handle missing points based on mode if (hint.num_missing_ > 0 && hint.mode_ == OutOfBoundsMode::FILL) { + auto fill_val = fill_value_; Kokkos::parallel_for( "FillMissingValues", Kokkos::RangePolicy(0, hint.num_missing_), - KOKKOS_LAMBDA(LO i) { - values[hint.missing_indices_(i)] = hint.fill_value_; - }); + KOKKOS_LAMBDA(LO i) { values[hint.missing_indices_(i)] = fill_val; }); } }