From 4c2d75b1c39cfdea40d2d05d239d0a8257ba6106 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Tue, 13 Jan 2026 18:30:42 -0500 Subject: [PATCH 1/3] create uniform grid field --- .../mesh_to_uniform_grid_workflow.cpp | 337 +++++++++ src/CMakeLists.txt | 8 +- .../uniform_grid/uniform_grid_field.cpp | 232 ++++++ .../adapter/uniform_grid/uniform_grid_field.h | 54 ++ .../uniform_grid_field_layout.cpp | 162 ++++ .../uniform_grid/uniform_grid_field_layout.h | 54 ++ src/pcms/create_field.cpp | 70 ++ src/pcms/create_field.h | 63 ++ src/pcms/uniform_grid.h | 56 ++ test/CMakeLists.txt | 1 + test/test_uniform_grid.cpp | 120 +++ test/test_uniform_grid_field.cpp | 708 ++++++++++++++++++ 12 files changed, 1864 insertions(+), 1 deletion(-) create mode 100644 examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp create mode 100644 src/pcms/adapter/uniform_grid/uniform_grid_field.cpp create mode 100644 src/pcms/adapter/uniform_grid/uniform_grid_field.h create mode 100644 src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp create mode 100644 src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h create mode 100644 test/test_uniform_grid_field.cpp diff --git a/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp b/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp new file mode 100644 index 00000000..4f944f5c --- /dev/null +++ b/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp @@ -0,0 +1,337 @@ +/** + * @file mesh_to_uniform_grid_workflow.cpp + * @brief Complete workflow: Transfer field data from Omega_h mesh to uniform + * grid with mask + * + * This example demonstrates a complete workflow starting from: + * 1. An Omega_h mesh with field data + * 2. Creating a uniform grid that covers the mesh domain + * 3. Creating a binary mask field (inside/outside mesh) + * 4. Transferring field data from mesh to uniform grid via interpolation + * 5. Using the mask to identify valid grid regions + * 6. Evaluating the grid field at arbitrary points + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Helper function to print field statistics +void print_field_stats(const std::string& name, + const std::vector& data) +{ + auto min = *std::min_element(data.begin(), data.end()); + auto max = *std::max_element(data.begin(), data.end()); + auto sum = std::accumulate(data.begin(), data.end(), 0.0); + + std::cout << name << " statistics:\n"; + std::cout << " Size: " << data.size() << "\n"; + std::cout << " Min: " << min << "\n"; + std::cout << " Max: " << max << "\n"; +} + +// Helper function to visualize mask field (2D only) +void print_mask_visualization(const pcms::UniformGrid<2>& grid, + const std::vector& mask, + int sample_size = 20) +{ + // Only visualize if grid is reasonable size + if (grid.divisions[0] > sample_size || grid.divisions[1] > sample_size) { + std::cout << "(Grid too large for visualization, skipping...)\n"; + return; + } + + std::cout << "\nMask field visualization (1=inside, 0=outside):\n"; + std::cout << std::string(grid.divisions[0] * 2 + 1, '-') << "\n"; + + // Print from top to bottom + for (int j = grid.divisions[1] - 1; j >= 0; --j) { + std::cout << "|"; + for (int i = 0; i < grid.divisions[0]; ++i) { + auto cell_id = grid.GetCellIndex({i, j}); + std::cout << mask[cell_id] << "|"; + } + std::cout << "\n"; + } + std::cout << std::string(grid.divisions[0] * 2 + 1, '-') << "\n"; +} + +int main(int argc, char** argv) +{ + // Initialize Omega_h + auto lib = Omega_h::Library(&argc, &argv); + auto world = lib.world(); + + std::cout << "========================================================\n"; + std::cout << "Mesh-to-UniformGrid Field Transfer Workflow\n"; + std::cout << "========================================================\n\n"; + + // ============================================================================ + // STEP 1: Create an Omega_h mesh + // ============================================================================ + std::cout << "STEP 1: Creating Omega_h mesh\n"; + std::cout << "------------------------------\n"; + + // Create a 2D mesh: domain [0, 2] x [0, 2] with triangular elements + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 2.0, 0.0, 8, 8, 0, false); + + std::cout << "Created 2D mesh:\n"; + std::cout << " Domain: [0, 2] x [0, 2]\n"; + std::cout << " Elements: " << mesh.nelems() << "\n"; + std::cout << " Vertices: " << mesh.nverts() << "\n\n"; + + // ============================================================================ + // STEP 2: Create an Omega_h field with mathematical data + // ============================================================================ + std::cout << "STEP 2: Creating Omega_h field with data\n"; + std::cout << "-----------------------------------------\n"; + + // Create a field layout on the mesh (P1 linear elements) + auto omega_h_layout = + pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto omega_h_field = omega_h_layout->CreateField(); + + // Get mesh vertex coordinates + auto mesh_coords = omega_h_layout->GetDOFHolderCoordinates(); + auto mesh_coords_data = mesh_coords.GetCoordinates(); + int num_mesh_nodes = omega_h_layout->GetNumOwnedDofHolder(); + + // Initialize field with a mathematical function: f(x,y) = sin(πx) * cos(πy) + std::cout << "Setting field values: f(x,y) = sin(πx) * cos(πy)\n"; + std::vector mesh_field_data(num_mesh_nodes); + for (int i = 0; i < num_mesh_nodes; ++i) { + pcms::Real x = mesh_coords_data(i, 0); + pcms::Real y = mesh_coords_data(i, 1); + mesh_field_data[i] = std::sin(M_PI * x) * std::cos(M_PI * y); + } + + auto mesh_data_view = + pcms::Rank1View( + mesh_field_data.data(), mesh_field_data.size()); + omega_h_field->SetDOFHolderData(mesh_data_view); + + print_field_stats("Omega_h field", mesh_field_data); + std::cout << "\n"; + + // ============================================================================ + // STEP 3: Create uniform grid from mesh + // ============================================================================ + std::cout << "STEP 3: Creating uniform grid from mesh\n"; + std::cout << "----------------------------------------\n"; + + // Create a 10x10 uniform grid covering the mesh bounding box + auto grid = pcms::CreateUniformGridFromMesh<2>(mesh, {10, 10}); + + std::cout << "Created uniform grid:\n"; + std::cout << " Number of cells: " << grid.GetNumCells() << "\n"; + std::cout << " Number of vertices: " + << (grid.divisions[0] + 1) * (grid.divisions[1] + 1) << "\n"; + std::cout << " Domain: [" << grid.bot_left[0] << ", " + << grid.bot_left[0] + grid.edge_length[0] << "] x [" + << grid.bot_left[1] << ", " + << grid.bot_left[1] + grid.edge_length[1] << "]\n"; + std::cout << " Cell size: " + << grid.edge_length[0] / grid.divisions[0] << " x " + << grid.edge_length[1] / grid.divisions[1] << "\n\n"; + + // ============================================================================ + // STEP 4: Create binary mask field (inside/outside mesh) + // ============================================================================ + std::cout << "STEP 4: Creating binary mask field\n"; + std::cout << "-----------------------------------\n"; + + auto mask_field = pcms::CreateUniformGridBinaryField<2>(mesh, {10, 10}); + + int inside_count = std::accumulate(mask_field.begin(), mask_field.end(), 0); + int outside_count = mask_field.size() - inside_count; + + std::cout << "Mask field created:\n"; + std::cout << " Total cells: " << mask_field.size() << "\n"; + std::cout << " Inside mesh: " << inside_count << " cells (" + << 100.0 * inside_count / mask_field.size() << "%)\n"; + std::cout << " Outside mesh: " << outside_count << " cells (" + << 100.0 * outside_count / mask_field.size() << "%)\n"; + + print_mask_visualization(grid, mask_field); + std::cout << "\n"; + + // ============================================================================ + // STEP 5: Create uniform grid field and transfer data + // ============================================================================ + std::cout << "STEP 5: Transferring field data to uniform grid\n"; + std::cout << "------------------------------------------------\n"; + + // Create uniform grid field layout + pcms::UniformGridFieldLayout<2> ug_layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto ug_field = ug_layout.CreateField(); + + std::cout << "Grid field layout:\n"; + std::cout << " DOF holders (vertices): " << ug_layout.GetNumOwnedDofHolder() + << "\n"; + std::cout << " Components: " << ug_layout.GetNumComponents() + << "\n\n"; + + // Get uniform grid vertex coordinates for evaluation + auto ug_coords = ug_layout.GetDOFHolderCoordinates(); + + // Evaluate omega_h field at uniform grid vertices + std::cout << "Interpolating mesh field to grid vertices...\n"; + std::vector ug_field_data(ug_layout.GetNumOwnedDofHolder()); + auto ug_data_view = pcms::Rank1View( + ug_field_data.data(), ug_field_data.size()); + pcms::FieldDataView field_data_view{ + ug_data_view, omega_h_field->GetCoordinateSystem()}; + + auto localization_hint = omega_h_field->GetLocalizationHint(ug_coords); + omega_h_field->Evaluate(localization_hint, field_data_view); + + // Set the evaluated data on the uniform grid field + auto ug_data_const_view = + pcms::Rank1View( + ug_field_data.data(), ug_field_data.size()); + ug_field->SetDOFHolderData(ug_data_const_view); + + print_field_stats("Uniform grid field", ug_field_data); + std::cout << "\n"; + + // ============================================================================ + // STEP 6: Verify the transferred data + // ============================================================================ + std::cout << "STEP 6: Verifying data transfer accuracy\n"; + std::cout << "-----------------------------------------\n"; + + auto transferred_data = ug_field->GetDOFHolderData(); + auto ug_coords_data = ug_coords.GetCoordinates(); + + // Check several sample points + std::vector sample_indices = {0, 5, 10, 50, 100, 120}; // Sample vertices + double max_error = 0.0; + + std::cout << "Sample verification points:\n"; + std::cout << std::fixed << std::setprecision(6); + for (int idx : sample_indices) { + if (idx >= ug_layout.GetNumOwnedDofHolder()) + continue; + + pcms::Real x = ug_coords_data(idx, 0); + pcms::Real y = ug_coords_data(idx, 1); + pcms::Real expected = std::sin(M_PI * x) * std::cos(M_PI * y); + pcms::Real actual = transferred_data[idx]; + double error = std::abs(actual - expected); + max_error = std::max(max_error, error); + + std::cout << " Vertex " << std::setw(3) << idx << " at (" << std::setw(8) + << x << ", " << std::setw(8) << y + << "): " << "expected=" << std::setw(10) << expected + << ", actual=" << std::setw(10) << actual + << ", error=" << std::scientific << error << std::fixed << "\n"; + } + + std::cout << "\nMaximum interpolation error: " << std::scientific << max_error + << std::fixed << "\n\n"; + + // ============================================================================ + // STEP 7: Use mask field to process only valid grid regions + // ============================================================================ + std::cout << "STEP 7: Processing valid grid regions using mask\n"; + std::cout << "-------------------------------------------------\n"; + + double sum_inside = 0.0; + int count_inside = 0; + + for (int j = 0; j < grid.divisions[1]; ++j) { + for (int i = 0; i < grid.divisions[0]; ++i) { + pcms::LO cell_id = grid.GetCellIndex({i, j}); + + if (mask_field[cell_id] == 1) { + // This cell is inside the mesh - we can safely work with it + auto cell_bbox = grid.GetCellBBOX(cell_id); + + // For demonstration, accumulate field values at cell centers + // In practice, you would evaluate the field or perform other operations + count_inside++; + + // Get the four vertices of this cell (bottom-left, bottom-right, + // top-left, top-right) + int v_bl = j * (grid.divisions[0] + 1) + i; + int v_br = j * (grid.divisions[0] + 1) + (i + 1); + int v_tl = (j + 1) * (grid.divisions[0] + 1) + i; + int v_tr = (j + 1) * (grid.divisions[0] + 1) + (i + 1); + + // Average the field values at the four corners + double cell_avg = (transferred_data[v_bl] + transferred_data[v_br] + + transferred_data[v_tl] + transferred_data[v_tr]) / + 4.0; + sum_inside += cell_avg; + } + } + } + + std::cout << "Processed " << count_inside << " valid cells (inside mesh)\n"; + std::cout << "Average field value in valid region: " + << sum_inside / count_inside << "\n\n"; + + // ============================================================================ + // STEP 8: Evaluate grid field at arbitrary points + // ============================================================================ + std::cout << "STEP 8: Evaluating grid field at arbitrary points\n"; + std::cout << "--------------------------------------------------\n"; + + // Define some test points + std::vector eval_points = { + 0.5, 0.5, // Point 1 + 1.0, 1.0, // Point 2 (center) + 1.5, 1.5, // Point 3 + 0.25, 1.75 // Point 4 + }; + + auto eval_coords_view = + pcms::Rank2View(eval_points.data(), + 4, 2); + auto eval_coord_view = pcms::CoordinateView( + pcms::CoordinateSystem::Cartesian, eval_coords_view); + + auto eval_hint = ug_field->GetLocalizationHint(eval_coord_view); + + std::vector eval_results(4); + auto eval_results_view = + pcms::Rank1View(eval_results.data(), 4); + auto eval_field_view = pcms::FieldDataView( + eval_results_view, pcms::CoordinateSystem::Cartesian); + + ug_field->Evaluate(eval_hint, eval_field_view); + + std::cout << "Evaluation results:\n"; + for (int i = 0; i < 4; ++i) { + pcms::Real x = eval_points[2 * i]; + pcms::Real y = eval_points[2 * i + 1]; + pcms::Real expected = std::sin(M_PI * x) * std::cos(M_PI * y); + pcms::Real actual = eval_results[i]; + + std::cout << " Point " << i + 1 << " (" << std::setw(8) << x << ", " + << std::setw(8) << y << "): " << "interpolated=" << std::setw(10) + << actual << ", exact=" << std::setw(10) << expected + << ", error=" << std::scientific << std::abs(actual - expected) + << std::fixed << "\n"; + } + + std::cout << "\n"; + + // ============================================================================ + // End + // ============================================================================ + std::cout << "Workflow Complete!\n"; + + return 0; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0d174940..3c6addad 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -54,7 +54,11 @@ if(PCMS_ENABLE_XGC) list(APPEND PCMS_HEADERS pcms/adapter/xgc/xgc_reverse_classification.h) endif() if(PCMS_ENABLE_OMEGA_H) - list(APPEND PCMS_SOURCES pcms/point_search.cpp) + list(APPEND PCMS_SOURCES + pcms/point_search.cpp + pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp + pcms/adapter/uniform_grid/uniform_grid_field.cpp + ) list( APPEND PCMS_HEADERS @@ -63,6 +67,8 @@ if(PCMS_ENABLE_OMEGA_H) pcms/transfer_field2.h pcms/uniform_grid.h pcms/point_search.h + pcms/adapter/uniform_grid/uniform_grid_field_layout.h + pcms/adapter/uniform_grid/uniform_grid_field.h ) endif() diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp new file mode 100644 index 00000000..3ae3c3a4 --- /dev/null +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp @@ -0,0 +1,232 @@ +#include "uniform_grid_field.h" +#include "pcms/profile.h" +#include "pcms/interpolator/linear_interpolant.hpp" +#include "pcms/interpolator/multidimarray.hpp" +#include "pcms/assert.h" +#include + +namespace pcms +{ + +// Localization hint structure for UniformGrid +template +struct UniformGridFieldLocalizationHint +{ + UniformGridFieldLocalizationHint( + Kokkos::View cell_indices, + Kokkos::View coordinates) + : cell_indices_(cell_indices), coordinates_(coordinates) + { + } + + // Cell indices for each point as 1d array + Kokkos::View cell_indices_; + // Coordinates of each point + Kokkos::View coordinates_; +}; + +template +UniformGridField::UniformGridField( + const UniformGridFieldLayout& layout) + : layout_(layout), + grid_(layout.GetGrid()), + dof_holder_data_("dof_holder_data", static_cast(layout.OwnedSize())) +{ + PCMS_FUNCTION_TIMER; +} + +template +Rank1View UniformGridField::GetDOFHolderData() + const +{ + PCMS_FUNCTION_TIMER; + return make_const_array_view(dof_holder_data_); +} + +template +void UniformGridField::SetDOFHolderData( + Rank1View data) +{ + PCMS_FUNCTION_TIMER; + PCMS_ALWAYS_ASSERT(data.size() == dof_holder_data_.size()); + + for (size_t i = 0; i < data.size(); ++i) { + dof_holder_data_[i] = data[i]; + } +} + +template +LocalizationHint UniformGridField::GetLocalizationHint( + CoordinateView coordinate_view) const +{ + PCMS_FUNCTION_TIMER; + + if (coordinate_view.GetCoordinateSystem() != + layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + throw std::runtime_error( + "Coordinate system mismatch in GetLocalizationHint"); + } + + auto coordinates = coordinate_view.GetCoordinates(); + LO num_points = coordinates.extent(0); + + Kokkos::View cell_indices("cell_indices", num_points); + Kokkos::View coords_copy("coords_copy", num_points, + Dim); + + // Find which cell each point belongs to + for (LO i = 0; i < num_points; ++i) { + Omega_h::Vector point; + for (unsigned d = 0; d < Dim; ++d) { + point[d] = coordinates(i, d); + coords_copy(i, d) = coordinates(i, d); + } + cell_indices[i] = grid_.ClosestCellID(point); + } + + auto hint = std::make_shared>( + cell_indices, coords_copy); + return LocalizationHint{hint}; +} + +template +void UniformGridField::Evaluate( + LocalizationHint location, FieldDataView results) const +{ + PCMS_FUNCTION_TIMER; + + if (results.GetCoordinateSystem() != + layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + throw std::runtime_error("Coordinate system mismatch in Evaluate"); + } + + UniformGridFieldLocalizationHint hint = + *reinterpret_cast*>( + location.data.get()); + + auto values = dof_holder_data_; + auto coordinates = hint.coordinates_; + auto cell_indices = hint.cell_indices_; + LO num_points = coordinates.extent(0); + + Kokkos::View cell_dimensioned_indices( + "cell_dimensioned_indices", num_points, Dim); + // Convert cell indices to multi-dimensional indices + for (LO i = 0; i < num_points; ++i) { + auto dim_indices = grid_.GetDimensionedIndex(cell_indices[i]); + for (unsigned d = 0; d < Dim; ++d) { + cell_dimensioned_indices(i, d) = dim_indices[d]; + } + } + + // Dimensions for vertex grid (m+1 vertices per dimension for m cells) + auto cell_divisions = layout_.GetGrid().divisions; + IntVecView dimensions_view("dimensions", Dim); + auto dimensions_view_host = Kokkos::create_mirror_view(dimensions_view); + for (unsigned d = 0; d < Dim; ++d) { + dimensions_view_host(d) = cell_divisions[d] + 1; + } + Kokkos::deep_copy(dimensions_view, dimensions_view_host); + + // Compute parametric coordinates for each point + RealMatView parametric_coords("parametric_coords", num_points, Dim); + auto parametric_coords_host = Kokkos::create_mirror_view(parametric_coords); + + for (LO i = 0; i < num_points; ++i) { + auto cell_bbox = grid_.GetCellBBOX(cell_indices[i]); + for (unsigned d = 0; d < Dim; ++d) { + Real coord = coordinates(i, d); + Real cell_min = cell_bbox.center[d] - cell_bbox.half_width[d]; + Real cell_max = cell_bbox.center[d] + cell_bbox.half_width[d]; + parametric_coords_host(i, d) = (coord - cell_min) / (cell_max - cell_min); + } + } + + Kokkos::deep_copy(parametric_coords, parametric_coords_host); + + // Convert cell indices and values to the required view types + IntMatView cell_indices_interp("cell_indices_interp", num_points, Dim); + auto cell_indices_interp_host = + Kokkos::create_mirror_view(cell_indices_interp); + for (LO i = 0; i < num_points; ++i) { + for (unsigned d = 0; d < Dim; ++d) { + cell_indices_interp_host(i, d) = cell_dimensioned_indices(i, d); + } + } + Kokkos::deep_copy(cell_indices_interp, cell_indices_interp_host); + + RealVecView values_interp("values_interp", values.extent(0)); + Kokkos::deep_copy(values_interp, values); + + auto interpolator = RegularGridInterpolator( + parametric_coords, values_interp, cell_indices_interp, dimensions_view); + auto evaluated_values = results.GetValues(); + auto interpolated_values = interpolator.linear_interpolation(); + auto interpolated_values_host = + Kokkos::create_mirror_view(interpolated_values); + Kokkos::deep_copy(interpolated_values_host, interpolated_values); + for (LO i = 0; i < evaluated_values.size(); ++i) { + evaluated_values[i] = interpolated_values_host[i]; + } +} + +template +void UniformGridField::EvaluateGradient( + FieldDataView) +{ + throw std::runtime_error("Not implemented"); +} + +template +const FieldLayout& UniformGridField::GetLayout() const +{ + return layout_; +} + +template +bool UniformGridField::CanEvaluateGradient() +{ + return false; +} + +template +int UniformGridField::Serialize( + Rank1View buffer, + Rank1View permutation) const +{ + PCMS_FUNCTION_TIMER; + + const auto array_h = GetDOFHolderData(); + if (buffer.size() > 0) { + PCMS_ALWAYS_ASSERT(buffer.size() == array_h.size()); + for (LO i = 0; i < array_h.size(); ++i) { + buffer[permutation[i]] = array_h[i]; + } + } + return array_h.size(); +} + +template +void UniformGridField::Deserialize( + Rank1View buffer, + Rank1View permutation) +{ + PCMS_FUNCTION_TIMER; + + Kokkos::View sorted_buffer("sorted_buffer", + permutation.size()); + auto owned = layout_.GetOwned(); + + for (LO i = 0; i < sorted_buffer.size(); ++i) { + PCMS_ALWAYS_ASSERT(owned[i]); + sorted_buffer[i] = buffer[permutation[i]]; + } + + SetDOFHolderData(pcms::make_const_array_view(sorted_buffer)); +} + +// Explicit template instantiations +template class UniformGridField<2>; +template class UniformGridField<3>; + +} // namespace pcms diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.h b/src/pcms/adapter/uniform_grid/uniform_grid_field.h new file mode 100644 index 00000000..49006b0d --- /dev/null +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.h @@ -0,0 +1,54 @@ +#ifndef PCMS_UNIFORM_GRID_FIELD_H +#define PCMS_UNIFORM_GRID_FIELD_H + +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/types.h" +#include "pcms/field.h" +#include "pcms/coordinate_system.h" +#include "pcms/uniform_grid.h" +#include + +namespace pcms +{ +template +class UniformGridField : public FieldT +{ +public: + UniformGridField(const UniformGridFieldLayout& layout); + + LocalizationHint GetLocalizationHint( + CoordinateView coordinate_view) const override; + + void Evaluate(LocalizationHint location, + FieldDataView results) const override; + + void EvaluateGradient(FieldDataView results) override; + + const FieldLayout& GetLayout() const override; + + bool CanEvaluateGradient() override; + + int Serialize(Rank1View buffer, + Rank1View permutation) + const override; + + void Deserialize( + Rank1View buffer, + Rank1View permutation) override; + + Rank1View GetDOFHolderData() const override; + void SetDOFHolderData(Rank1View data) override; + + ~UniformGridField() noexcept = default; + +private: + const UniformGridFieldLayout& layout_; + UniformGrid& grid_; + Kokkos::View dof_holder_data_; +}; + +using UniformGridField2D = UniformGridField<2>; + +} // namespace pcms + +#endif // PCMS_UNIFORM_GRID_FIELD_H diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp new file mode 100644 index 00000000..e8a52315 --- /dev/null +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp @@ -0,0 +1,162 @@ +#include "uniform_grid_field_layout.h" +#include "uniform_grid_field.h" +#include "pcms/profile.h" +#include + +namespace pcms +{ + +template +UniformGridFieldLayout::UniformGridFieldLayout( + UniformGrid& grid, int num_components, + CoordinateSystem coordinate_system) + : grid_(grid), + num_components_(num_components), + coordinate_system_(coordinate_system), + gids_("gids", GetNumVertices()), + dof_holder_coords_("dof_holder_coords", GetNumVertices(), Dim), + owned_("owned", GetNumVertices()) +{ + PCMS_FUNCTION_TIMER; + + LO num_vertices = GetNumVertices(); + + // Initialize global IDs and ownership + for (LO i = 0; i < num_vertices; ++i) { + gids_[i] = static_cast(i); + owned_[i] = true; + } + + // Initialize DOF holder coordinates at grid vertices + Real vertex_spacing[Dim]; + for (unsigned d = 0; d < Dim; ++d) { + vertex_spacing[d] = grid_.edge_length[d] / grid_.divisions[d]; + } + + if constexpr (Dim == 2) { + LO vertex_idx = 0; + for (LO j = 0; j <= grid_.divisions[1]; ++j) { + for (LO i = 0; i <= grid_.divisions[0]; ++i) { + dof_holder_coords_(vertex_idx, 0) = + grid_.bot_left[0] + i * vertex_spacing[0]; + dof_holder_coords_(vertex_idx, 1) = + grid_.bot_left[1] + j * vertex_spacing[1]; + ++vertex_idx; + } + } + } else if constexpr (Dim == 3) { + LO vertex_idx = 0; + for (LO k = 0; k <= grid_.divisions[2]; ++k) { + for (LO j = 0; j <= grid_.divisions[1]; ++j) { + for (LO i = 0; i <= grid_.divisions[0]; ++i) { + dof_holder_coords_(vertex_idx, 0) = + grid_.bot_left[0] + i * vertex_spacing[0]; + dof_holder_coords_(vertex_idx, 1) = + grid_.bot_left[1] + j * vertex_spacing[1]; + dof_holder_coords_(vertex_idx, 2) = + grid_.bot_left[2] + k * vertex_spacing[2]; + ++vertex_idx; + } + } + } + } +} + +template +std::unique_ptr> UniformGridFieldLayout::CreateField() const +{ + return std::make_unique>(*this); +} + +template +int UniformGridFieldLayout::GetNumComponents() const +{ + return num_components_; +} + +template +LO UniformGridFieldLayout::GetNumOwnedDofHolder() const +{ + return GetNumVertices(); +} + +template +GO UniformGridFieldLayout::GetNumGlobalDofHolder() const +{ + return GetNumVertices(); +} + +template +Rank1View UniformGridFieldLayout::GetOwned() + const +{ + return make_const_array_view(owned_); +} + +template +GlobalIDView UniformGridFieldLayout::GetGids() const +{ + return GlobalIDView(gids_.data(), gids_.size()); +} + +template +CoordinateView +UniformGridFieldLayout::GetDOFHolderCoordinates() const +{ + Rank2View coords_view( + dof_holder_coords_.data(), dof_holder_coords_.extent(0), Dim); + return CoordinateView{coordinate_system_, coords_view}; +} + +template +bool UniformGridFieldLayout::IsDistributed() +{ + return false; +} + +template +UniformGrid& UniformGridFieldLayout::GetGrid() const +{ + return grid_; +} + +template +LO UniformGridFieldLayout::GetNumCells() const +{ + return grid_.GetNumCells(); +} + +template +LO UniformGridFieldLayout::GetNumVertices() const +{ + LO num_vertices = 1; + for (unsigned d = 0; d < Dim; ++d) { + num_vertices *= (grid_.divisions[d] + 1); + } + return num_vertices; +} + +template +EntOffsetsArray UniformGridFieldLayout::GetEntOffsets() const +{ + EntOffsetsArray offsets{}; + offsets[0] = 0; + offsets[1] = grid_.GetNumCells(); + offsets[2] = grid_.GetNumCells(); + offsets[3] = grid_.GetNumCells(); + offsets[4] = grid_.GetNumCells(); + return offsets; +} + +template +ReversePartitionMap2 UniformGridFieldLayout::GetReversePartitionMap( + const redev::Partition& partition) const +{ + throw std::runtime_error("Unimplemented"); +} + +// Explicit template instantiations +template class UniformGridFieldLayout<2>; +template class UniformGridFieldLayout<3>; + +} // namespace pcms diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h new file mode 100644 index 00000000..c87c7733 --- /dev/null +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h @@ -0,0 +1,54 @@ +#ifndef PCMS_UNIFORM_GRID_FIELD_LAYOUT_H +#define PCMS_UNIFORM_GRID_FIELD_LAYOUT_H + +#include "pcms/arrays.h" +#include "pcms/field_layout.h" +#include "pcms/coordinate_system.h" +#include "pcms/field.h" +#include "pcms/uniform_grid.h" + +#include + +namespace pcms +{ +template +class UniformGridFieldLayout : public FieldLayout +{ +public: + UniformGridFieldLayout(UniformGrid& grid, int num_components, + CoordinateSystem coordinate_system); + + std::unique_ptr> CreateField() const override; + + int GetNumComponents() const override; + LO GetNumOwnedDofHolder() const override; + GO GetNumGlobalDofHolder() const override; + + Rank1View GetOwned() const override; + GlobalIDView GetGids() const override; + CoordinateView GetDOFHolderCoordinates() const override; + + bool IsDistributed() override; + + EntOffsetsArray GetEntOffsets() const override; + + ReversePartitionMap2 GetReversePartitionMap( + const redev::Partition& partition) const override; + + UniformGrid& GetGrid() const; + LO GetNumCells() const; + LO GetNumVertices() const; + +private: + UniformGrid& grid_; + int num_components_; + CoordinateSystem coordinate_system_; + Kokkos::View gids_; + Kokkos::View dof_holder_coords_; + Kokkos::View owned_; +}; + +using UniformGridFieldLayout2D = UniformGridFieldLayout<2>; + +} // namespace pcms +#endif // PCMS_UNIFORM_GRID_FIELD_LAYOUT_H diff --git a/src/pcms/create_field.cpp b/src/pcms/create_field.cpp index 95db5d8f..f9fb3809 100644 --- a/src/pcms/create_field.cpp +++ b/src/pcms/create_field.cpp @@ -1,7 +1,10 @@ #include "create_field.h" #include "adapter/omega_h/omega_h_field2.h" #include "adapter/omega_h/omega_h_field_layout.h" +#include "uniform_grid.h" +#include "point_search.h" +#include #include namespace pcms @@ -24,4 +27,71 @@ std::unique_ptr CreateLagrangeLayout( num_components, coordinate_system); } +template <> +std::vector CreateUniformGridBinaryFieldFromGrid<2>( + Omega_h::Mesh& mesh, const UniformGrid<2>& grid) +{ + constexpr unsigned dim = 2; + + // Get total number of cells + const LO num_cells = grid.GetNumCells(); + + // Initialize result vector + std::vector binary_field(num_cells, 0); + + // Create GridPointSearch for point-in-mesh queries + GridPointSearch point_search(mesh, grid.divisions[0], grid.divisions[1]); + + // Create array of grid cell center points + Kokkos::View cell_centers("cell centers", num_cells); + auto cell_centers_h = Kokkos::create_mirror_view(cell_centers); + + // Fill cell center coordinates + for (LO i = 0; i < num_cells; ++i) { + auto bbox = grid.GetCellBBOX(i); + for (unsigned d = 0; d < dim; ++d) { + cell_centers_h(i, d) = bbox.center[d]; + } + } + + // Copy to device + Kokkos::deep_copy(cell_centers, cell_centers_h); + + // Perform point search + auto results = point_search(cell_centers); + + // Copy results back to host + auto results_h = Kokkos::create_mirror_view(results); + Kokkos::deep_copy(results_h, results); + + // Process results: set field to 1 if point is inside mesh (tri_id >= 0) + for (LO i = 0; i < num_cells; ++i) { + binary_field[i] = (results_h(i).tri_id >= 0) ? 1 : 0; + } + + return binary_field; +} + +template <> +std::vector CreateUniformGridBinaryField<2>( + Omega_h::Mesh& mesh, const std::array& divisions) +{ + constexpr unsigned dim = 2; + + // Create the uniform grid from the mesh + auto grid = CreateUniformGridFromMesh(mesh, divisions); + + // Delegate to the grid-based implementation + return CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid); +} + +template <> +std::vector CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh, + LO cells_per_dim) +{ + std::array divisions; + divisions.fill(cells_per_dim); + return CreateUniformGridBinaryField<2>(mesh, divisions); +} + } // namespace pcms diff --git a/src/pcms/create_field.h b/src/pcms/create_field.h index 8cf7ba35..b595d1bf 100644 --- a/src/pcms/create_field.h +++ b/src/pcms/create_field.h @@ -5,15 +5,78 @@ #include "field_layout.h" #include "field.h" #include "coordinate_system.h" +#include "types.h" +#include "uniform_grid.h" #include +#include #include +#include namespace pcms { std::unique_ptr CreateLagrangeLayout( Omega_h::Mesh& mesh, int order, int num_components, CoordinateSystem coordinate_system); + +/** + * \brief Create a binary field on a uniform grid indicating inside/outside mesh + * + * This function creates a uniform grid from the mesh and generates a binary + * field where each grid cell is assigned: + * - 1 if the cell center is inside the original mesh + * - 0 if the cell center is outside the original mesh + * + * Uses GridPointSearch to determine if a point lies within the mesh domain. + * Currently only supports 2D meshes. + * + * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported) + * \param mesh The Omega_h mesh to create the grid from + * \param divisions Array specifying the number of cells in each dimension + * \return std::vector Binary field values (0 or 1) for each grid cell + * + * \note The returned vector has length equal to the number of grid cells. + * Cell values are ordered according to the grid's internal indexing. + */ +template +std::vector CreateUniformGridBinaryField( + Omega_h::Mesh& mesh, const std::array& divisions); + +/** + * \brief Create a binary field with equal divisions in all dimensions + * + * Convenience function that creates a uniform grid with the same number of + * cells in each dimension and generates the binary inside/outside field. + * + * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported) + * \param mesh The Omega_h mesh to create the grid from + * \param cells_per_dim Number of cells per dimension (same for all dimensions) + * \return std::vector Binary field values (0 or 1) for each grid cell + */ +template +std::vector CreateUniformGridBinaryField(Omega_h::Mesh& mesh, + LO cells_per_dim); + +/** + * \brief Create a binary field on a given uniform grid + * + * This function takes a pre-defined uniform grid and generates a binary field + * where each grid cell is assigned: + * - 1 if the cell center is inside the mesh + * - 0 if the cell center is outside the mesh + * + * This allows testing with custom grids that may extend beyond the mesh + * boundaries. + * + * \tparam dim Spatial dimension (currently only dim=2 is supported) + * \param mesh The Omega_h mesh to test against + * \param grid The uniform grid to evaluate + * \return std::vector Binary field values (0 or 1) for each grid cell + */ +template +std::vector CreateUniformGridBinaryFieldFromGrid( + Omega_h::Mesh& mesh, const UniformGrid& grid); + } // namespace pcms #endif // CREATE_FIELD_H_ diff --git a/src/pcms/uniform_grid.h b/src/pcms/uniform_grid.h index 5ec13d19..f1320f1d 100644 --- a/src/pcms/uniform_grid.h +++ b/src/pcms/uniform_grid.h @@ -2,6 +2,8 @@ #define PCMS_COUPLING_UNIFORM_GRID_H #include "pcms/bounding_box.h" #include "Omega_h_vector.hpp" +#include "Omega_h_bbox.hpp" +#include "Omega_h_mesh.hpp" #include namespace pcms { @@ -119,6 +121,60 @@ struct UniformGrid using Uniform2DGrid = UniformGrid<>; +/** + * \brief Create a uniform grid layout from an Omega_h mesh as a bounding box + * + * This function computes the bounding box of the given mesh and creates a + * uniform Cartesian grid that covers the entire mesh domain. The grid cells + * are uniformly distributed across each dimension. + * + * \tparam dim Spatial dimension of the mesh (2 or 3) + * \param mesh The Omega_h mesh to create the grid from + * \param divisions Array specifying the number of cells in each dimension + * \return UniformGrid A uniform grid covering the mesh bounding box + * + */ +template +UniformGrid CreateUniformGridFromMesh(Omega_h::Mesh& mesh, + const std::array& divisions) +{ + // Get the bounding box of the mesh + auto bbox = Omega_h::get_bounding_box(&mesh); + + // Calculate edge lengths and bottom-left corner + std::array edge_length; + std::array bot_left; + + for (unsigned i = 0; i < dim; ++i) { + bot_left[i] = bbox.min[i]; + edge_length[i] = bbox.max[i] - bbox.min[i]; + } + + return UniformGrid{ + .edge_length = edge_length, .bot_left = bot_left, .divisions = divisions}; +} + +/** + * \brief Create a uniform grid with equal divisions in all dimensions + * + * This is a convenience function that creates a uniform grid with the same + * number of cells in each dimension. + * + * \tparam dim Spatial dimension of the mesh (2 or 3) + * \param mesh The Omega_h mesh to create the grid from + * \param cells_per_dim Number of cells per dimension (same for all dimensions) + * \return UniformGrid A uniform grid covering the mesh bounding box + * + */ +template +UniformGrid CreateUniformGridFromMesh(Omega_h::Mesh& mesh, + LO cells_per_dim) +{ + std::array divisions; + divisions.fill(cells_per_dim); + return CreateUniformGridFromMesh(mesh, divisions); +} + } // namespace pcms #endif // PCMS_COUPLING_UNIFORM_GRID_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4325dab2..bfe767af 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -394,6 +394,7 @@ if(Catch2_FOUND) test_normalisation.cpp test_spline_interpolator.cpp test_svd_serial.cpp + test_uniform_grid_field.cpp test_interpolation_class.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) diff --git a/test/test_uniform_grid.cpp b/test/test_uniform_grid.cpp index 9cdc3f39..df145357 100644 --- a/test/test_uniform_grid.cpp +++ b/test/test_uniform_grid.cpp @@ -1,8 +1,12 @@ #include #include #include +#include +#include +using pcms::CreateUniformGridFromMesh; using pcms::Uniform2DGrid; +using pcms::UniformGrid; TEST_CASE("uniform grid") { @@ -50,3 +54,119 @@ TEST_CASE("uniform grid") REQUIRE(119 == uniform_grid.GetCellIndex({11, 9})); } } + +TEST_CASE("Create uniform grid from omega_h mesh") +{ + // Initialize Omega_h library + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + // Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, false); + + SECTION("Create grid with specified divisions per dimension") + { + // Create a 10x12 uniform grid over the mesh + auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 12}); + + // Verify the grid properties + REQUIRE(grid.GetNumCells() == 120); + REQUIRE(grid.edge_length[0] == Catch::Approx(1.0)); + REQUIRE(grid.edge_length[1] == Catch::Approx(1.0)); + REQUIRE(grid.bot_left[0] == Catch::Approx(0.0)); + REQUIRE(grid.bot_left[1] == Catch::Approx(0.0)); + REQUIRE(grid.divisions[0] == 10); + REQUIRE(grid.divisions[1] == 12); + + // Test cell bounding boxes + auto bbox_first = grid.GetCellBBOX(0); + REQUIRE(bbox_first.center[0] == Catch::Approx(0.05)); // 1.0/10/2 + REQUIRE(bbox_first.center[1] == Catch::Approx(0.5 / 12.0)); + + auto bbox_last = grid.GetCellBBOX(119); + REQUIRE(bbox_last.center[0] == Catch::Approx(0.95)); + REQUIRE(bbox_last.center[1] == Catch::Approx(11.5 / 12.0)); + } + + SECTION("Create grid with equal divisions") + { + // Create a 15x15 uniform grid + auto grid = CreateUniformGridFromMesh<2>(mesh, 15); + + REQUIRE(grid.GetNumCells() == 225); + REQUIRE(grid.divisions[0] == 15); + REQUIRE(grid.divisions[1] == 15); + + // Verify all cells have equal size + auto cell_width = grid.edge_length[0] / grid.divisions[0]; + auto cell_height = grid.edge_length[1] / grid.divisions[1]; + REQUIRE(cell_width == Catch::Approx(1.0 / 15.0)); + REQUIRE(cell_height == Catch::Approx(1.0 / 15.0)); + } + + SECTION("Test with non-unit domain") + { + // Create a larger domain: 5.0 x 3.0 + auto large_mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 5.0, 3.0, 0.0, + 10, 6, 0, false); + + auto grid = CreateUniformGridFromMesh<2>(large_mesh, {20, 30}); + + REQUIRE(grid.GetNumCells() == 600); + REQUIRE(grid.edge_length[0] == Catch::Approx(5.0)); + REQUIRE(grid.edge_length[1] == Catch::Approx(3.0)); + REQUIRE(grid.bot_left[0] == Catch::Approx(0.0)); + REQUIRE(grid.bot_left[1] == Catch::Approx(0.0)); + + // Test cell sizes + auto bbox = grid.GetCellBBOX(0); + REQUIRE(bbox.half_width[0] == Catch::Approx(5.0 / 20.0 / 2.0)); + REQUIRE(bbox.half_width[1] == Catch::Approx(3.0 / 30.0 / 2.0)); + } + + SECTION("Test point location in grid cells") + { + auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 10}); + + // Test corner points + REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{0.0, 0.0}) == 0); + REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{1.0, 1.0}) == 99); + + // Test middle point + REQUIRE(grid.ClosestCellID(Omega_h::Vector<2>{0.5, 0.5}) == 55); + } +} + +TEST_CASE("Create uniform grid from 3D mesh") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + // Create a 3D box mesh: 2.0 x 3.0 x 1.5 domain + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 3.0, 1.5, 4, 6, 3, false); + + SECTION("Create 3D grid with specified divisions") + { + auto grid = CreateUniformGridFromMesh<3>(mesh, {10, 15, 8}); + + REQUIRE(grid.GetNumCells() == 1200); // 10 * 15 * 8 + REQUIRE(grid.edge_length[0] == Catch::Approx(2.0)); + REQUIRE(grid.edge_length[1] == Catch::Approx(3.0)); + REQUIRE(grid.edge_length[2] == Catch::Approx(1.5)); + REQUIRE(grid.divisions[0] == 10); + REQUIRE(grid.divisions[1] == 15); + REQUIRE(grid.divisions[2] == 8); + } + + SECTION("Create 3D grid with equal divisions") + { + auto grid = CreateUniformGridFromMesh<3>(mesh, 12); + + REQUIRE(grid.GetNumCells() == 1728); // 12^3 + REQUIRE(grid.divisions[0] == 12); + REQUIRE(grid.divisions[1] == 12); + REQUIRE(grid.divisions[2] == 12); + } +} diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp new file mode 100644 index 00000000..c1d11f05 --- /dev/null +++ b/test/test_uniform_grid_field.cpp @@ -0,0 +1,708 @@ +#include +#include +#include +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field.h" +#include "pcms/uniform_grid.h" +#include "Omega_h_library.hpp" +#include "Omega_h_build.hpp" +#include "pcms/transfer_field2.h" +#include "pcms/adapter/omega_h/omega_h_field_layout.h" +#include "pcms/adapter/omega_h/omega_h_field2.h" +#include "pcms/create_field.h" +#include "pcms/arrays.h" +#include + +using pcms::CreateUniformGridBinaryField; +using pcms::CreateUniformGridFromMesh; + +TEST_CASE("UniformGrid field creation") +{ + // Create a simple 2D uniform grid + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {5, 5}; + + // Create field layout with 1 component (scalar field) + pcms::UniformGridFieldLayout<2> layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + + REQUIRE(layout.GetNumComponents() == 1); + REQUIRE(layout.GetNumOwnedDofHolder() == 36); // (5+1)x(5+1) = 36 vertices + REQUIRE(layout.GetNumGlobalDofHolder() == 36); + REQUIRE_FALSE(layout.IsDistributed()); + + // Create field + auto field = layout.CreateField(); + REQUIRE(field != nullptr); +} + +TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") +{ + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {4, 4}; // 4x4 = 16 cells + + pcms::UniformGridFieldLayout<2> layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto field = layout.CreateField(); + + // Initialize field data with vertex indices (5x5 = 25 vertices for 4x4 cells) + std::vector data(25); + for (size_t i = 0; i < 25; ++i) { + data[i] = static_cast(i); + } + + auto data_view = pcms::Rank1View( + data.data(), data.size()); + field->SetDOFHolderData(data_view); + + // Retrieve data + auto retrieved = field->GetDOFHolderData(); + REQUIRE(retrieved.size() == 25); + + for (size_t i = 0; i < 25; ++i) { + REQUIRE(retrieved[i] == static_cast(i)); + } +} + +TEST_CASE("UniformGrid field evaluation - piecewise constant") +{ + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {2, 2}; // 2x2 = 4 cells + + pcms::UniformGridFieldLayout<2> layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto field = layout.CreateField(); + + // Set vertex values for a 2x2 cell grid (3x3 = 9 vertices) + // Vertex layout: + // v6---v7---v8 + // | 2 | 3 | + // v3---v4---v5 + // | 0 | 1 | + // v0---v1---v2 + std::vector data = { + 1.0, 1.5, 2.0, // v0, v1, v2 (bottom row, y=0) + 2.0, 2.5, 3.0, // v3, v4, v5 (middle row, y=5) + 3.0, 3.5, 4.0 // v6, v7, v8 (top row, y=10) + }; + auto data_view = pcms::Rank1View( + data.data(), data.size()); + field->SetDOFHolderData(data_view); + + // Evaluate at cell centers + std::vector eval_coords = { + 2.5, 2.5, // Cell 0 center + 7.5, 2.5, // Cell 1 center + 2.5, 7.5, // Cell 2 center + 7.5, 7.5 // Cell 3 center + }; + + auto coords_view = pcms::Rank2View( + eval_coords.data(), 4, 2); + auto coord_view = pcms::CoordinateView( + pcms::CoordinateSystem::Cartesian, coords_view); + + auto hint = field->GetLocalizationHint(coord_view); + + std::vector results(4); + auto results_view = + pcms::Rank1View(results.data(), 4); + auto results_field_view = + pcms::FieldDataView( + results_view, pcms::CoordinateSystem::Cartesian); + + field->Evaluate(hint, results_field_view); + + // Check results - interpolated from vertices + // Cell 0 center (2.5, 2.5): avg of v0,v1,v3,v4 = (1.0+1.5+2.0+2.5)/4 = 1.75 + // Cell 1 center (7.5, 2.5): avg of v1,v2,v4,v5 = (1.5+2.0+2.5+3.0)/4 = 2.25 + // Cell 2 center (2.5, 7.5): avg of v3,v4,v6,v7 = (2.0+2.5+3.0+3.5)/4 = 2.75 + // Cell 3 center (7.5, 7.5): avg of v4,v5,v7,v8 = (2.5+3.0+3.5+4.0)/4 = 3.25 + REQUIRE(std::abs(results[0] - 1.75) < 1e-10); + REQUIRE(std::abs(results[1] - 2.25) < 1e-10); + REQUIRE(std::abs(results[2] - 2.75) < 1e-10); + REQUIRE(std::abs(results[3] - 3.25) < 1e-10); +} + +TEST_CASE("UniformGrid field serialization") +{ + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {3, 3}; // 9 cells + + pcms::UniformGridFieldLayout<2> layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto field = layout.CreateField(); + + // Set field data (4x4 = 16 vertices for 3x3 cells) + std::vector data(16); + for (size_t i = 0; i < 16; ++i) { + data[i] = static_cast(i * 10); + } + + auto data_view = pcms::Rank1View( + data.data(), data.size()); + field->SetDOFHolderData(data_view); + + // Create identity permutation + std::vector permutation(16); + for (size_t i = 0; i < 16; ++i) { + permutation[i] = i; + } + + auto perm_view = pcms::Rank1View( + permutation.data(), 16); + + // Serialize + std::vector buffer(16); + auto buffer_view = + pcms::Rank1View(buffer.data(), 16); + int size = field->Serialize(buffer_view, perm_view); + + REQUIRE(size == 16); + + // Verify serialized data + for (size_t i = 0; i < 16; ++i) { + REQUIRE(buffer[i] == data[i]); + } + + // Create new field and deserialize + auto field2 = layout.CreateField(); + auto buffer_const_view = + pcms::Rank1View(buffer.data(), 16); + field2->Deserialize(buffer_const_view, perm_view); + + // Verify deserialized data + auto retrieved = field2->GetDOFHolderData(); + for (size_t i = 0; i < 16; ++i) { + REQUIRE(retrieved[i] == data[i]); + } +} + +TEST_CASE("UniformGrid field copy") +{ + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {2, 2}; // 2x2 grid + + pcms::UniformGridFieldLayout<2> layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto field = layout.CreateField(); + + // Set vertex values with f(x,y) = x + y at 3x3 vertex positions + // Vertices at: (0,0), (5,0), (10,0), (0,5), (5,5), (10,5), (0,10), (5,10), + // (10,10) + std::vector data = { + 0.0, 5.0, 10.0, // y=0: v0(0,0)=0, v1(5,0)=5, v2(10,0)=10 + 5.0, 10.0, 15.0, // y=5: v3(0,5)=5, v4(5,5)=10, v5(10,5)=15 + 10.0, 15.0, 20.0 // y=10: v6(0,10)=10, v7(5,10)=15, v8(10,10)=20 + }; + auto data_view = pcms::Rank1View( + data.data(), data.size()); + field->SetDOFHolderData(data_view); + + // Test 1: Evaluate at cell centers + std::vector eval_coords = { + 2.5, 2.5, // Cell 0 center + 7.5, 2.5, // Cell 1 center + 2.5, 7.5, // Cell 2 center + 7.5, 7.5 // Cell 3 center + }; + + auto coords_view = pcms::Rank2View( + eval_coords.data(), 4, 2); + auto coord_view = pcms::CoordinateView( + pcms::CoordinateSystem::Cartesian, coords_view); + + auto hint = field->GetLocalizationHint(coord_view); + + std::vector results(coord_view.GetCoordinates().size() / 2); + auto results_view = + pcms::Rank1View(results.data(), 4); + auto results_field_view = + pcms::FieldDataView( + results_view, pcms::CoordinateSystem::Cartesian); + + field->Evaluate(hint, results_field_view); + + REQUIRE(std::abs(results[0] - 5.0) < 1e-10); + REQUIRE(std::abs(results[1] - 10.0) < 1e-10); + REQUIRE(std::abs(results[2] - 10.0) < 1e-10); + REQUIRE(std::abs(results[3] - 15.0) < 1e-10); + + // Test 2: test copy_field2 + auto field2 = layout.CreateField(); + pcms::copy_field2(*field, *field2); + + auto copied_data = field2->GetDOFHolderData(); + REQUIRE(copied_data.size() == data.size()); + for (size_t i = 0; i < data.size(); ++i) { + REQUIRE(copied_data[i] == data[i]); + } +} + +TEST_CASE("Transfer from OmegaH field to UniformGrid field") +{ + Omega_h::Library lib; + + // Create a simple omega_h mesh (2x2 box) + auto mesh = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 2, + 2, 0, false); + + // Create OmegaH field layout with linear elements + auto omega_h_layout = + pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto omega_h_field = omega_h_layout->CreateField(); + + // Initialize omega_h field with a simple function f(x,y) = x + 2*y + // Get coordinates from the layout + auto coords = omega_h_layout->GetDOFHolderCoordinates(); + auto coords_data = coords.GetCoordinates(); + int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); + + std::vector omega_h_data(num_nodes); + for (int i = 0; i < num_nodes; ++i) { + pcms::Real x = coords_data(i, 0); + pcms::Real y = coords_data(i, 1); + omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y + } + + auto omega_h_data_view = + pcms::Rank1View( + omega_h_data.data(), omega_h_data.size()); + omega_h_field->SetDOFHolderData(omega_h_data_view); + + // Create a uniform grid field covering the same domain [0,1] x [0,1] + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {1.0, 1.0}; + grid.divisions = {2, 2}; // 4x4 grid for finer resolution + + pcms::UniformGridFieldLayout<2> ug_layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto ug_field = ug_layout.CreateField(); + + // Transfer from omega_h field to uniform grid field using interpolation + auto coords_interpolation = ug_layout.GetDOFHolderCoordinates(); + std::vector evaluation( + coords_interpolation.GetCoordinates().size() / 2); + auto evaluation_view = pcms::Rank1View( + evaluation.data(), evaluation.size()); + pcms::FieldDataView data_view{ + evaluation_view, omega_h_field->GetCoordinateSystem()}; + auto locale = omega_h_field->GetLocalizationHint(coords_interpolation); + omega_h_field->Evaluate(locale, data_view); + auto evaluation_view_const = + pcms::Rank1View(evaluation.data(), + evaluation.size()); + ug_field->SetDOFHolderData(evaluation_view_const); + + // Verify the transferred data at uniform grid vertices + auto transferred_data = ug_field->GetDOFHolderData(); + auto ug_coords = ug_layout.GetDOFHolderCoordinates(); + auto ug_coords_data = ug_coords.GetCoordinates(); + int num_ug_nodes = ug_layout.GetNumOwnedDofHolder(); // 5x5 = 25 vertices + + // Check a few sample points + for (int i = 0; i < num_ug_nodes; ++i) { + pcms::Real x = ug_coords_data(i, 0); + pcms::Real y = ug_coords_data(i, 1); + pcms::Real expected = x + 2.0 * y; + pcms::Real actual = transferred_data[i]; + + // Allow some tolerance for interpolation + REQUIRE(std::abs(actual - expected) < 1e-6); + } +} + +TEST_CASE("Create binary field from uniform grid") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + SECTION("Simple 2D box mesh - all cells inside") + { + // Create a mesh that fills the domain [0,1] x [0,1] + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, + 10, 0, false); + + // Create a 5x5 grid (coarser than mesh) + auto field = CreateUniformGridBinaryField<2>(mesh, {5, 5}); + + REQUIRE(field.size() == 25); + + // All grid cells should be inside the mesh + int sum = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(sum == 25); + } + + SECTION("Binary field with custom divisions") + { + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 8, 8, 0, false); + + auto field = CreateUniformGridBinaryField<2>(mesh, {10, 8}); + + REQUIRE(field.size() == 80); + + // Most cells should be inside + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + REQUIRE(inside_count <= 80); + } + + SECTION("Binary field with equal divisions convenience function") + { + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); + + auto field = CreateUniformGridBinaryField<2>(mesh, 8); + + REQUIRE(field.size() == 64); + + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + } + + SECTION("Fine grid over coarse mesh") + { + // Create a simple mesh + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 2.0, 0.0, 4, 4, 0, false); + + // Create a fine grid (20x20) + auto grid = CreateUniformGridFromMesh<2>(mesh, {20, 20}); + auto field = CreateUniformGridBinaryField<2>(mesh, {20, 20}); + + REQUIRE(field.size() == 400); + + // Verify consistency: check some specific cells + // Center cells should be inside + auto center_idx = grid.GetCellIndex({10, 10}); + REQUIRE(field[center_idx] == 1); + + // Corner cells should be inside + auto corner_idx = grid.GetCellIndex({0, 0}); + REQUIRE(field[corner_idx] == 1); + } + + SECTION("Verify field values are binary") + { + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); + + auto field = CreateUniformGridBinaryField<2>(mesh, 10); + + // All values should be 0 or 1 + for (const auto& val : field) { + REQUIRE((val == 0 || val == 1)); + } + } + + SECTION("Grid extends beyond mesh - cells outside should be marked 0") + { + // Create a small mesh in the center of a domain + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.5, 0.5, 0.0, 5, 5, 0, false); + + // The mesh occupies [0, 0.5] x [0, 0.5] + // Create a grid that would cover this + auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 10}); + auto field = CreateUniformGridBinaryField<2>(mesh, {10, 10}); + + REQUIRE(field.size() == 100); + + // All cells should be inside since grid is exactly on mesh bbox + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + } + + SECTION("Test with different aspect ratio") + { + // Create a rectangular mesh + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 3.0, 1.0, 0.0, 12, 4, + 0, false); + + auto field = CreateUniformGridBinaryField<2>(mesh, {30, 10}); + + REQUIRE(field.size() == 300); + + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + + // Calculate percentage inside + double inside_percent = 100.0 * inside_count / field.size(); + // Most cells should be inside + REQUIRE(inside_percent > 50.0); + } +} + +TEST_CASE("Binary field integration with grid methods") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, 10, 0, false); + + SECTION("Query field value at specific grid cell") + { + auto grid = CreateUniformGridFromMesh<2>(mesh, {8, 8}); + auto field = CreateUniformGridBinaryField<2>(mesh, {8, 8}); + + // Get field value for a specific cell + pcms::LO cell_id = grid.GetCellIndex({4, 4}); // Middle cell + REQUIRE(field[cell_id] == 1); // Should be inside + + // Get bbox of that cell + auto bbox = grid.GetCellBBOX(cell_id); + // Verify it's roughly in the middle + REQUIRE(bbox.center[0] == Catch::Approx(0.5).margin(0.1)); + REQUIRE(bbox.center[1] == Catch::Approx(0.5).margin(0.1)); + } + + SECTION("Count cells by region") + { + auto grid = CreateUniformGridFromMesh<2>(mesh, 10); + auto field = CreateUniformGridBinaryField<2>(mesh, 10); + + // Count cells in different quadrants + int q1 = 0, q2 = 0, q3 = 0, q4 = 0; // quadrants + + for (pcms::LO i = 0; i < grid.GetNumCells(); ++i) { + if (field[i] == 1) { + auto bbox = grid.GetCellBBOX(i); + if (bbox.center[0] < 0.5 && bbox.center[1] < 0.5) + q1++; + else if (bbox.center[0] >= 0.5 && bbox.center[1] < 0.5) + q2++; + else if (bbox.center[0] < 0.5 && bbox.center[1] >= 0.5) + q3++; + else + q4++; + } + } + + // All quadrants should have some inside cells + REQUIRE(q1 > 0); + REQUIRE(q2 > 0); + REQUIRE(q3 > 0); + REQUIRE(q4 > 0); + + // Distribution should be roughly equal + int total = q1 + q2 + q3 + q4; + REQUIRE(q1 == Catch::Approx(total / 4.0).margin(5)); + REQUIRE(q2 == Catch::Approx(total / 4.0).margin(5)); + REQUIRE(q3 == Catch::Approx(total / 4.0).margin(5)); + REQUIRE(q4 == Catch::Approx(total / 4.0).margin(5)); + } +} + +TEST_CASE("Performance and edge cases") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + SECTION("Very fine grid") + { + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); + + // Create a very fine grid + auto field = CreateUniformGridBinaryField<2>(mesh, 50); + + REQUIRE(field.size() == 2500); + + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + } + + SECTION("Coarse grid") + { + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, + 10, 0, false); + + // Very coarse grid + auto field = CreateUniformGridBinaryField<2>(mesh, {2, 2}); + + REQUIRE(field.size() == 4); + + // All 4 cells should be inside for this configuration + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + } + + SECTION("Non-square domain") + { + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 5.0, 2.0, 0.0, 20, 8, + 0, false); + + auto field = CreateUniformGridBinaryField<2>(mesh, {25, 10}); + + REQUIRE(field.size() == 250); + + int inside_count = std::accumulate(field.begin(), field.end(), 0); + REQUIRE(inside_count > 0); + } + + SECTION("Grid larger than mesh - cells outside marked as 0") + { + // Create a mesh covering [0, 0.5] x [0, 0.5] + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.5, 0.5, 0.0, 5, 5, 0, false); + + // Manually create a larger grid covering [0, 1] x [0, 1] + pcms::UniformGrid<2> grid; + grid.edge_length = {1.0, 1.0}; + grid.bot_left = {0.0, 0.0}; + grid.divisions = {10, 10}; + + // Create binary field on the larger grid + auto field = pcms::CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid); + + REQUIRE(field.size() == 100); + + // Count cells inside and outside + int inside_count = std::accumulate(field.begin(), field.end(), 0); + int outside_count = field.size() - inside_count; + + // Should have both inside (1) and outside (0) cells + REQUIRE(inside_count > 0); + REQUIRE(outside_count > 0); + + // Check specific cells + // Bottom-left corner [0, 0.1] x [0, 0.1] should be inside (mesh covers [0, + // 0.5]) + auto corner_id = grid.GetCellIndex({0, 0}); + auto bbox = grid.GetCellBBOX(corner_id); + REQUIRE(field[corner_id] == 1); + + // Top-right corner [0.9, 1.0] x [0.9, 1.0] should be outside (mesh ends at + // 0.5) + corner_id = grid.GetCellIndex({9, 9}); + bbox = grid.GetCellBBOX(corner_id); + REQUIRE(field[corner_id] == 0); + + // Cell at [0.5, 0.6] x [0.5, 0.6] should be outside + corner_id = grid.GetCellIndex({6, 6}); + bbox = grid.GetCellBBOX(corner_id); + REQUIRE(field[corner_id] == 0); + + // Cell at [0.2, 0.3] x [0.2, 0.3] should be inside + auto center_id = grid.GetCellIndex({2, 2}); + bbox = grid.GetCellBBOX(center_id); + REQUIRE(field[center_id] == 1); + } +} + +TEST_CASE("UniformGrid workflow") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + + // Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, false); + auto grid = pcms::CreateUniformGridFromMesh<2>(mesh, {4, 4}); + auto mask_field = pcms::CreateUniformGridBinaryField<2>(mesh, {4, 4}); + + // Create OmegaH field layout with linear elements + auto omega_h_layout = + pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto omega_h_field = omega_h_layout->CreateField(); + + // Initialize omega_h field with a simple function f(x,y) = x + 2*y + // Get coordinates from the layout + auto coords = omega_h_layout->GetDOFHolderCoordinates(); + auto coords_data = coords.GetCoordinates(); + int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); + + std::vector omega_h_data(num_nodes); + for (int i = 0; i < num_nodes; ++i) { + pcms::Real x = coords_data(i, 0); + pcms::Real y = coords_data(i, 1); + omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y + } + + auto omega_h_data_view = + pcms::Rank1View( + omega_h_data.data(), omega_h_data.size()); + omega_h_field->SetDOFHolderData(omega_h_data_view); + + // Create uniform grid field layout + pcms::UniformGridFieldLayout<2> ug_layout(grid, 1, + pcms::CoordinateSystem::Cartesian); + auto ug_field = ug_layout.CreateField(); + + auto coords_interpolation = ug_layout.GetDOFHolderCoordinates(); + std::vector evaluation( + coords_interpolation.GetCoordinates().size() / 2); + auto evaluation_view = pcms::Rank1View( + evaluation.data(), evaluation.size()); + pcms::FieldDataView data_view{ + evaluation_view, omega_h_field->GetCoordinateSystem()}; + auto locale = omega_h_field->GetLocalizationHint(coords_interpolation); + omega_h_field->Evaluate(locale, data_view); + auto evaluation_view_const = + pcms::Rank1View(evaluation.data(), + evaluation.size()); + ug_field->SetDOFHolderData(evaluation_view_const); + + // Verify transferred data at grid vertices using evaluation array + auto ug_coords = ug_layout.GetDOFHolderCoordinates(); + for (int j = 0; j < grid.divisions[1]; ++j) { + for (int i = 0; i < grid.divisions[0]; ++i) { + int vertex_id = j * (grid.divisions[0] + 1) + i; + pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0); + pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1); + pcms::Real expected = x + 2.0 * y; + pcms::Real actual = evaluation[vertex_id]; + fprintf(stderr, + "Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n", i, + j, x, y, expected, actual); + if (std::abs(expected - actual) > 1e-10) { + throw std::runtime_error("Data mismatch at vertex (" + + std::to_string(i) + ", " + std::to_string(j) + + "): expected " + std::to_string(expected) + + ", got " + std::to_string(actual)); + } + } + } + + // Verify ug_field values directly from the field object + fprintf(stderr, "\nVerifying ug_field values:\n"); + auto ug_field_data = ug_field->GetDOFHolderData(); + for (int j = 0; j <= grid.divisions[1]; ++j) { + for (int i = 0; i <= grid.divisions[0]; ++i) { + int vertex_id = j * (grid.divisions[0] + 1) + i; + pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0); + pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1); + pcms::Real expected = x + 2.0 * y; + pcms::Real actual = ug_field_data(vertex_id); + fprintf( + stderr, + "ug_field Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n", + i, j, x, y, expected, actual); + REQUIRE(std::abs(expected - actual) <= 1e-10); + } + } + // Verify mask field values + fprintf(stderr, "\nVerifying mask_field values:\n"); + auto mask_data = mask_field.data(); + for (int j = 0; j < grid.divisions[1]; ++j) { + for (int i = 0; i < grid.divisions[0]; ++i) { + int cell_id = j * grid.divisions[0] + i; + int mask_value = mask_data[cell_id]; + fprintf(stderr, "Cell (%d, %d): mask = %d\n", i, j, mask_value); + REQUIRE(mask_value == 1); + } + } +} From 26f4ad427aa4bdeb89eadf1ba61f4e77664d7481 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Wed, 21 Jan 2026 19:09:42 -0500 Subject: [PATCH 2/3] clean uniform tests --- .../mesh_to_uniform_grid_workflow.cpp | 337 ------------------ src/pcms/transfer_field2.h | 4 +- test/test_uniform_grid_field.cpp | 134 +++---- 3 files changed, 60 insertions(+), 415 deletions(-) delete mode 100644 examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp diff --git a/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp b/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp deleted file mode 100644 index 4f944f5c..00000000 --- a/examples/mesh-to-uniform-grid-example/mesh_to_uniform_grid_workflow.cpp +++ /dev/null @@ -1,337 +0,0 @@ -/** - * @file mesh_to_uniform_grid_workflow.cpp - * @brief Complete workflow: Transfer field data from Omega_h mesh to uniform - * grid with mask - * - * This example demonstrates a complete workflow starting from: - * 1. An Omega_h mesh with field data - * 2. Creating a uniform grid that covers the mesh domain - * 3. Creating a binary mask field (inside/outside mesh) - * 4. Transferring field data from mesh to uniform grid via interpolation - * 5. Using the mask to identify valid grid regions - * 6. Evaluating the grid field at arbitrary points - */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// Helper function to print field statistics -void print_field_stats(const std::string& name, - const std::vector& data) -{ - auto min = *std::min_element(data.begin(), data.end()); - auto max = *std::max_element(data.begin(), data.end()); - auto sum = std::accumulate(data.begin(), data.end(), 0.0); - - std::cout << name << " statistics:\n"; - std::cout << " Size: " << data.size() << "\n"; - std::cout << " Min: " << min << "\n"; - std::cout << " Max: " << max << "\n"; -} - -// Helper function to visualize mask field (2D only) -void print_mask_visualization(const pcms::UniformGrid<2>& grid, - const std::vector& mask, - int sample_size = 20) -{ - // Only visualize if grid is reasonable size - if (grid.divisions[0] > sample_size || grid.divisions[1] > sample_size) { - std::cout << "(Grid too large for visualization, skipping...)\n"; - return; - } - - std::cout << "\nMask field visualization (1=inside, 0=outside):\n"; - std::cout << std::string(grid.divisions[0] * 2 + 1, '-') << "\n"; - - // Print from top to bottom - for (int j = grid.divisions[1] - 1; j >= 0; --j) { - std::cout << "|"; - for (int i = 0; i < grid.divisions[0]; ++i) { - auto cell_id = grid.GetCellIndex({i, j}); - std::cout << mask[cell_id] << "|"; - } - std::cout << "\n"; - } - std::cout << std::string(grid.divisions[0] * 2 + 1, '-') << "\n"; -} - -int main(int argc, char** argv) -{ - // Initialize Omega_h - auto lib = Omega_h::Library(&argc, &argv); - auto world = lib.world(); - - std::cout << "========================================================\n"; - std::cout << "Mesh-to-UniformGrid Field Transfer Workflow\n"; - std::cout << "========================================================\n\n"; - - // ============================================================================ - // STEP 1: Create an Omega_h mesh - // ============================================================================ - std::cout << "STEP 1: Creating Omega_h mesh\n"; - std::cout << "------------------------------\n"; - - // Create a 2D mesh: domain [0, 2] x [0, 2] with triangular elements - auto mesh = - Omega_h::build_box(world, OMEGA_H_SIMPLEX, 2.0, 2.0, 0.0, 8, 8, 0, false); - - std::cout << "Created 2D mesh:\n"; - std::cout << " Domain: [0, 2] x [0, 2]\n"; - std::cout << " Elements: " << mesh.nelems() << "\n"; - std::cout << " Vertices: " << mesh.nverts() << "\n\n"; - - // ============================================================================ - // STEP 2: Create an Omega_h field with mathematical data - // ============================================================================ - std::cout << "STEP 2: Creating Omega_h field with data\n"; - std::cout << "-----------------------------------------\n"; - - // Create a field layout on the mesh (P1 linear elements) - auto omega_h_layout = - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); - auto omega_h_field = omega_h_layout->CreateField(); - - // Get mesh vertex coordinates - auto mesh_coords = omega_h_layout->GetDOFHolderCoordinates(); - auto mesh_coords_data = mesh_coords.GetCoordinates(); - int num_mesh_nodes = omega_h_layout->GetNumOwnedDofHolder(); - - // Initialize field with a mathematical function: f(x,y) = sin(πx) * cos(πy) - std::cout << "Setting field values: f(x,y) = sin(πx) * cos(πy)\n"; - std::vector mesh_field_data(num_mesh_nodes); - for (int i = 0; i < num_mesh_nodes; ++i) { - pcms::Real x = mesh_coords_data(i, 0); - pcms::Real y = mesh_coords_data(i, 1); - mesh_field_data[i] = std::sin(M_PI * x) * std::cos(M_PI * y); - } - - auto mesh_data_view = - pcms::Rank1View( - mesh_field_data.data(), mesh_field_data.size()); - omega_h_field->SetDOFHolderData(mesh_data_view); - - print_field_stats("Omega_h field", mesh_field_data); - std::cout << "\n"; - - // ============================================================================ - // STEP 3: Create uniform grid from mesh - // ============================================================================ - std::cout << "STEP 3: Creating uniform grid from mesh\n"; - std::cout << "----------------------------------------\n"; - - // Create a 10x10 uniform grid covering the mesh bounding box - auto grid = pcms::CreateUniformGridFromMesh<2>(mesh, {10, 10}); - - std::cout << "Created uniform grid:\n"; - std::cout << " Number of cells: " << grid.GetNumCells() << "\n"; - std::cout << " Number of vertices: " - << (grid.divisions[0] + 1) * (grid.divisions[1] + 1) << "\n"; - std::cout << " Domain: [" << grid.bot_left[0] << ", " - << grid.bot_left[0] + grid.edge_length[0] << "] x [" - << grid.bot_left[1] << ", " - << grid.bot_left[1] + grid.edge_length[1] << "]\n"; - std::cout << " Cell size: " - << grid.edge_length[0] / grid.divisions[0] << " x " - << grid.edge_length[1] / grid.divisions[1] << "\n\n"; - - // ============================================================================ - // STEP 4: Create binary mask field (inside/outside mesh) - // ============================================================================ - std::cout << "STEP 4: Creating binary mask field\n"; - std::cout << "-----------------------------------\n"; - - auto mask_field = pcms::CreateUniformGridBinaryField<2>(mesh, {10, 10}); - - int inside_count = std::accumulate(mask_field.begin(), mask_field.end(), 0); - int outside_count = mask_field.size() - inside_count; - - std::cout << "Mask field created:\n"; - std::cout << " Total cells: " << mask_field.size() << "\n"; - std::cout << " Inside mesh: " << inside_count << " cells (" - << 100.0 * inside_count / mask_field.size() << "%)\n"; - std::cout << " Outside mesh: " << outside_count << " cells (" - << 100.0 * outside_count / mask_field.size() << "%)\n"; - - print_mask_visualization(grid, mask_field); - std::cout << "\n"; - - // ============================================================================ - // STEP 5: Create uniform grid field and transfer data - // ============================================================================ - std::cout << "STEP 5: Transferring field data to uniform grid\n"; - std::cout << "------------------------------------------------\n"; - - // Create uniform grid field layout - pcms::UniformGridFieldLayout<2> ug_layout(grid, 1, - pcms::CoordinateSystem::Cartesian); - auto ug_field = ug_layout.CreateField(); - - std::cout << "Grid field layout:\n"; - std::cout << " DOF holders (vertices): " << ug_layout.GetNumOwnedDofHolder() - << "\n"; - std::cout << " Components: " << ug_layout.GetNumComponents() - << "\n\n"; - - // Get uniform grid vertex coordinates for evaluation - auto ug_coords = ug_layout.GetDOFHolderCoordinates(); - - // Evaluate omega_h field at uniform grid vertices - std::cout << "Interpolating mesh field to grid vertices...\n"; - std::vector ug_field_data(ug_layout.GetNumOwnedDofHolder()); - auto ug_data_view = pcms::Rank1View( - ug_field_data.data(), ug_field_data.size()); - pcms::FieldDataView field_data_view{ - ug_data_view, omega_h_field->GetCoordinateSystem()}; - - auto localization_hint = omega_h_field->GetLocalizationHint(ug_coords); - omega_h_field->Evaluate(localization_hint, field_data_view); - - // Set the evaluated data on the uniform grid field - auto ug_data_const_view = - pcms::Rank1View( - ug_field_data.data(), ug_field_data.size()); - ug_field->SetDOFHolderData(ug_data_const_view); - - print_field_stats("Uniform grid field", ug_field_data); - std::cout << "\n"; - - // ============================================================================ - // STEP 6: Verify the transferred data - // ============================================================================ - std::cout << "STEP 6: Verifying data transfer accuracy\n"; - std::cout << "-----------------------------------------\n"; - - auto transferred_data = ug_field->GetDOFHolderData(); - auto ug_coords_data = ug_coords.GetCoordinates(); - - // Check several sample points - std::vector sample_indices = {0, 5, 10, 50, 100, 120}; // Sample vertices - double max_error = 0.0; - - std::cout << "Sample verification points:\n"; - std::cout << std::fixed << std::setprecision(6); - for (int idx : sample_indices) { - if (idx >= ug_layout.GetNumOwnedDofHolder()) - continue; - - pcms::Real x = ug_coords_data(idx, 0); - pcms::Real y = ug_coords_data(idx, 1); - pcms::Real expected = std::sin(M_PI * x) * std::cos(M_PI * y); - pcms::Real actual = transferred_data[idx]; - double error = std::abs(actual - expected); - max_error = std::max(max_error, error); - - std::cout << " Vertex " << std::setw(3) << idx << " at (" << std::setw(8) - << x << ", " << std::setw(8) << y - << "): " << "expected=" << std::setw(10) << expected - << ", actual=" << std::setw(10) << actual - << ", error=" << std::scientific << error << std::fixed << "\n"; - } - - std::cout << "\nMaximum interpolation error: " << std::scientific << max_error - << std::fixed << "\n\n"; - - // ============================================================================ - // STEP 7: Use mask field to process only valid grid regions - // ============================================================================ - std::cout << "STEP 7: Processing valid grid regions using mask\n"; - std::cout << "-------------------------------------------------\n"; - - double sum_inside = 0.0; - int count_inside = 0; - - for (int j = 0; j < grid.divisions[1]; ++j) { - for (int i = 0; i < grid.divisions[0]; ++i) { - pcms::LO cell_id = grid.GetCellIndex({i, j}); - - if (mask_field[cell_id] == 1) { - // This cell is inside the mesh - we can safely work with it - auto cell_bbox = grid.GetCellBBOX(cell_id); - - // For demonstration, accumulate field values at cell centers - // In practice, you would evaluate the field or perform other operations - count_inside++; - - // Get the four vertices of this cell (bottom-left, bottom-right, - // top-left, top-right) - int v_bl = j * (grid.divisions[0] + 1) + i; - int v_br = j * (grid.divisions[0] + 1) + (i + 1); - int v_tl = (j + 1) * (grid.divisions[0] + 1) + i; - int v_tr = (j + 1) * (grid.divisions[0] + 1) + (i + 1); - - // Average the field values at the four corners - double cell_avg = (transferred_data[v_bl] + transferred_data[v_br] + - transferred_data[v_tl] + transferred_data[v_tr]) / - 4.0; - sum_inside += cell_avg; - } - } - } - - std::cout << "Processed " << count_inside << " valid cells (inside mesh)\n"; - std::cout << "Average field value in valid region: " - << sum_inside / count_inside << "\n\n"; - - // ============================================================================ - // STEP 8: Evaluate grid field at arbitrary points - // ============================================================================ - std::cout << "STEP 8: Evaluating grid field at arbitrary points\n"; - std::cout << "--------------------------------------------------\n"; - - // Define some test points - std::vector eval_points = { - 0.5, 0.5, // Point 1 - 1.0, 1.0, // Point 2 (center) - 1.5, 1.5, // Point 3 - 0.25, 1.75 // Point 4 - }; - - auto eval_coords_view = - pcms::Rank2View(eval_points.data(), - 4, 2); - auto eval_coord_view = pcms::CoordinateView( - pcms::CoordinateSystem::Cartesian, eval_coords_view); - - auto eval_hint = ug_field->GetLocalizationHint(eval_coord_view); - - std::vector eval_results(4); - auto eval_results_view = - pcms::Rank1View(eval_results.data(), 4); - auto eval_field_view = pcms::FieldDataView( - eval_results_view, pcms::CoordinateSystem::Cartesian); - - ug_field->Evaluate(eval_hint, eval_field_view); - - std::cout << "Evaluation results:\n"; - for (int i = 0; i < 4; ++i) { - pcms::Real x = eval_points[2 * i]; - pcms::Real y = eval_points[2 * i + 1]; - pcms::Real expected = std::sin(M_PI * x) * std::cos(M_PI * y); - pcms::Real actual = eval_results[i]; - - std::cout << " Point " << i + 1 << " (" << std::setw(8) << x << ", " - << std::setw(8) << y << "): " << "interpolated=" << std::setw(10) - << actual << ", exact=" << std::setw(10) << expected - << ", error=" << std::scientific << std::abs(actual - expected) - << std::fixed << "\n"; - } - - std::cout << "\n"; - - // ============================================================================ - // End - // ============================================================================ - std::cout << "Workflow Complete!\n"; - - return 0; -} diff --git a/src/pcms/transfer_field2.h b/src/pcms/transfer_field2.h index 60a0266f..df5712ef 100644 --- a/src/pcms/transfer_field2.h +++ b/src/pcms/transfer_field2.h @@ -34,8 +34,8 @@ void interpolate_field2(const FieldT& source, FieldT& target) } auto coords = target.GetLayout().GetDOFHolderCoordinates(); - std::vector evaluation(coords.GetCoordinates().size() / 2); - FieldDataView data_view{make_array_view(evaluation), + std::vector evaluation(coords.GetCoordinates().size() / 2); + FieldDataView data_view{make_array_view(evaluation), source.GetCoordinateSystem()}; auto locale = source.GetLocalizationHint(coords); source.Evaluate(locale, data_view); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index c1d11f05..0ce22eba 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -16,6 +16,57 @@ using pcms::CreateUniformGridBinaryField; using pcms::CreateUniformGridFromMesh; +// Helper function to initialize omega_h field data with f(x,y) = x + 2*y +std::vector CreateOmegaHFieldData( + const pcms::CoordinateView& coords, int num_nodes) { + auto coords_data = coords.GetCoordinates(); + std::vector omega_h_data(num_nodes); + for (int i = 0; i < num_nodes; ++i) { + pcms::Real x = coords_data(i, 0); + pcms::Real y = coords_data(i, 1); + omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y + } + return omega_h_data; +} + +// Helper function to verify ug_field values +void VerifyUniformGridFieldValues( + const pcms::UniformGrid<2>& grid, + const pcms::CoordinateView& ug_coords, + const pcms::Rank1View& ug_field_data) { + fprintf(stderr, "\nVerifying ug_field values:\n"); + for (int j = 0; j <= grid.divisions[1]; ++j) { + for (int i = 0; i <= grid.divisions[0]; ++i) { + int vertex_id = j * (grid.divisions[0] + 1) + i; + pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0); + pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1); + pcms::Real expected = x + 2.0 * y; + pcms::Real actual = ug_field_data(vertex_id); + fprintf( + stderr, + "ug_field Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n", + i, j, x, y, expected, actual); + REQUIRE(std::abs(expected - actual) <= 1e-10); + } + } +} + +// Helper function to verify mask field values +void VerifyMaskFieldValues( + const pcms::UniformGrid<2>& grid, + const std::vector& mask_field) { + fprintf(stderr, "\nVerifying mask_field values:\n"); + auto mask_data = mask_field.data(); + for (int j = 0; j < grid.divisions[1]; ++j) { + for (int i = 0; i < grid.divisions[0]; ++i) { + int cell_id = j * grid.divisions[0] + i; + int mask_value = mask_data[cell_id]; + fprintf(stderr, "Cell (%d, %d): mask = %d\n", i, j, mask_value); + REQUIRE(mask_value == 1); + } + } +} + TEST_CASE("UniformGrid field creation") { // Create a simple 2D uniform grid @@ -263,17 +314,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") auto omega_h_field = omega_h_layout->CreateField(); // Initialize omega_h field with a simple function f(x,y) = x + 2*y - // Get coordinates from the layout auto coords = omega_h_layout->GetDOFHolderCoordinates(); - auto coords_data = coords.GetCoordinates(); int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); - - std::vector omega_h_data(num_nodes); - for (int i = 0; i < num_nodes; ++i) { - pcms::Real x = coords_data(i, 0); - pcms::Real y = coords_data(i, 1); - omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y - } + std::vector omega_h_data = CreateOmegaHFieldData(coords, num_nodes); auto omega_h_data_view = pcms::Rank1View( @@ -620,17 +663,9 @@ TEST_CASE("UniformGrid workflow") auto omega_h_field = omega_h_layout->CreateField(); // Initialize omega_h field with a simple function f(x,y) = x + 2*y - // Get coordinates from the layout auto coords = omega_h_layout->GetDOFHolderCoordinates(); - auto coords_data = coords.GetCoordinates(); int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); - - std::vector omega_h_data(num_nodes); - for (int i = 0; i < num_nodes; ++i) { - pcms::Real x = coords_data(i, 0); - pcms::Real y = coords_data(i, 1); - omega_h_data[i] = x + 2.0 * y; // f(x,y) = x + 2y - } + std::vector omega_h_data = CreateOmegaHFieldData(coords, num_nodes); auto omega_h_data_view = pcms::Rank1View( @@ -642,67 +677,14 @@ TEST_CASE("UniformGrid workflow") pcms::CoordinateSystem::Cartesian); auto ug_field = ug_layout.CreateField(); - auto coords_interpolation = ug_layout.GetDOFHolderCoordinates(); - std::vector evaluation( - coords_interpolation.GetCoordinates().size() / 2); - auto evaluation_view = pcms::Rank1View( - evaluation.data(), evaluation.size()); - pcms::FieldDataView data_view{ - evaluation_view, omega_h_field->GetCoordinateSystem()}; - auto locale = omega_h_field->GetLocalizationHint(coords_interpolation); - omega_h_field->Evaluate(locale, data_view); - auto evaluation_view_const = - pcms::Rank1View(evaluation.data(), - evaluation.size()); - ug_field->SetDOFHolderData(evaluation_view_const); - - // Verify transferred data at grid vertices using evaluation array + // Transfer from omega_h field to uniform grid field using interpolation + pcms::interpolate_field2(*omega_h_field, *ug_field); auto ug_coords = ug_layout.GetDOFHolderCoordinates(); - for (int j = 0; j < grid.divisions[1]; ++j) { - for (int i = 0; i < grid.divisions[0]; ++i) { - int vertex_id = j * (grid.divisions[0] + 1) + i; - pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0); - pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1); - pcms::Real expected = x + 2.0 * y; - pcms::Real actual = evaluation[vertex_id]; - fprintf(stderr, - "Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n", i, - j, x, y, expected, actual); - if (std::abs(expected - actual) > 1e-10) { - throw std::runtime_error("Data mismatch at vertex (" + - std::to_string(i) + ", " + std::to_string(j) + - "): expected " + std::to_string(expected) + - ", got " + std::to_string(actual)); - } - } - } // Verify ug_field values directly from the field object - fprintf(stderr, "\nVerifying ug_field values:\n"); auto ug_field_data = ug_field->GetDOFHolderData(); - for (int j = 0; j <= grid.divisions[1]; ++j) { - for (int i = 0; i <= grid.divisions[0]; ++i) { - int vertex_id = j * (grid.divisions[0] + 1) + i; - pcms::Real x = ug_coords.GetCoordinates()(vertex_id, 0); - pcms::Real y = ug_coords.GetCoordinates()(vertex_id, 1); - pcms::Real expected = x + 2.0 * y; - pcms::Real actual = ug_field_data(vertex_id); - fprintf( - stderr, - "ug_field Vertex (%d, %d) at (%.2f, %.2f): expected %.4f, got %.4f\n", - i, j, x, y, expected, actual); - REQUIRE(std::abs(expected - actual) <= 1e-10); - } - } + VerifyUniformGridFieldValues(grid, ug_coords, ug_field_data); + // Verify mask field values - fprintf(stderr, "\nVerifying mask_field values:\n"); - auto mask_data = mask_field.data(); - for (int j = 0; j < grid.divisions[1]; ++j) { - for (int i = 0; i < grid.divisions[0]; ++i) { - int cell_id = j * grid.divisions[0] + i; - int mask_value = mask_data[cell_id]; - fprintf(stderr, "Cell (%d, %d): mask = %d\n", i, j, mask_value); - REQUIRE(mask_value == 1); - } - } + VerifyMaskFieldValues(grid, mask_field); } From 0e0919c23c8aaea2d5f30f9817fd0372a6df0228 Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 6 Feb 2026 18:06:10 -0500 Subject: [PATCH 3/3] create pythonapi --- CMakeLists.txt | 4 + src/CMakeLists.txt | 4 + .../uniform_grid/uniform_grid_field.cpp | 38 ++ .../adapter/uniform_grid/uniform_grid_field.h | 3 + src/pcms/create_field.cpp | 76 ++- src/pcms/create_field.h | 42 +- src/pcms/pythonapi/CMakeLists.txt | 71 +++ src/pcms/pythonapi/FindPythonModule.cmake | 82 +++ src/pcms/pythonapi/bind_field_base.cpp | 282 ++++++++++ src/pcms/pythonapi/bind_field_layout.cpp | 88 +++ src/pcms/pythonapi/bind_omega_h.cpp | 526 ++++++++++++++++++ src/pcms/pythonapi/bind_omega_h_field2.cpp | 127 +++++ .../pythonapi/bind_omega_h_field_layout.cpp | 119 ++++ src/pcms/pythonapi/bind_transfer_field2.cpp | 62 +++ .../pythonapi/bind_uniform_grid_field.cpp | 230 ++++++++ .../bind_uniform_grid_field_layout.cpp | 222 ++++++++ src/pcms/pythonapi/numpy_array_transform.h | 232 ++++++++ src/pcms/pythonapi/pythonapi.cpp | 53 ++ src/pcms/pythonapi/test_field_copy.py | 84 +++ src/pcms/pythonapi/test_file_io.py | 452 +++++++++++++++ src/pcms/pythonapi/test_omega_h_field.py | 276 +++++++++ src/pcms/pythonapi/test_uniform_grid_field.py | 476 ++++++++++++++++ .../pythonapi/uniform_grid_workflow_guide.md | 264 +++++++++ test/test_uniform_grid_field.cpp | 308 +++++----- 24 files changed, 3949 insertions(+), 172 deletions(-) create mode 100644 src/pcms/pythonapi/CMakeLists.txt create mode 100644 src/pcms/pythonapi/FindPythonModule.cmake create mode 100644 src/pcms/pythonapi/bind_field_base.cpp create mode 100644 src/pcms/pythonapi/bind_field_layout.cpp create mode 100644 src/pcms/pythonapi/bind_omega_h.cpp create mode 100644 src/pcms/pythonapi/bind_omega_h_field2.cpp create mode 100644 src/pcms/pythonapi/bind_omega_h_field_layout.cpp create mode 100644 src/pcms/pythonapi/bind_transfer_field2.cpp create mode 100644 src/pcms/pythonapi/bind_uniform_grid_field.cpp create mode 100644 src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp create mode 100644 src/pcms/pythonapi/numpy_array_transform.h create mode 100644 src/pcms/pythonapi/pythonapi.cpp create mode 100644 src/pcms/pythonapi/test_field_copy.py create mode 100644 src/pcms/pythonapi/test_file_io.py create mode 100644 src/pcms/pythonapi/test_omega_h_field.py create mode 100644 src/pcms/pythonapi/test_uniform_grid_field.py create mode 100644 src/pcms/pythonapi/uniform_grid_workflow_guide.md diff --git a/CMakeLists.txt b/CMakeLists.txt index d9c519ca..72cfaf09 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,6 +24,7 @@ option(PCMS_ENABLE_CLIENT "enable the coupling client implementation" ON) option(PCMS_ENABLE_XGC "enable xgc field adapter" ON) option(PCMS_ENABLE_OMEGA_H "enable Omega_h field adapter" OFF) option(PCMS_ENABLE_C "Enable pcms C api" ON) +option(PCMS_ENABLE_Python "Enable pcms Python api" OFF) option(PCMS_ENABLE_PRINT "PCMS print statements enabled" ON) option(PCMS_ENABLE_SPDLOG "use spdlog for logging" ON) @@ -37,6 +38,9 @@ endif() # broken find_package(redev 4.3.0 REQUIRED) +if (PCMS_ENABLE_Python) + set(CMAKE_POSITION_INDEPENDENT_CODE TRUE) +endif () if(PCMS_ENABLE_C) enable_language(C) option(PCMS_ENABLE_Fortran "Enable pcms fortran api" ON) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3c6addad..18b6bc61 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -144,6 +144,10 @@ install( add_library(pcms_pcms INTERFACE) target_link_libraries(pcms_pcms INTERFACE pcms::core) set_target_properties(pcms_pcms PROPERTIES EXPORT_NAME pcms) +if (PCMS_ENABLE_Python) + add_subdirectory(pcms/pythonapi) + #target_link_libraries(pcms_pcms INTERFACE pcms::pythonapi) +endif() if(PCMS_ENABLE_C) add_subdirectory(pcms/capi) target_link_libraries(pcms_pcms INTERFACE pcms::capi) diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp index 3ae3c3a4..dc584b74 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp @@ -55,6 +55,44 @@ void UniformGridField::SetDOFHolderData( } } +template +View UniformGridField::to_mdspan() +{ + PCMS_FUNCTION_TIMER; + + if constexpr (Dim == 2) { + return View( + dof_holder_data_.data(), grid_.divisions[0] + 1, + grid_.divisions[1] + 1); + } else if constexpr (Dim == 3) { + return View( + dof_holder_data_.data(), grid_.divisions[0] + 1, + grid_.divisions[1] + 1, grid_.divisions[2] + 1); + } else { + static_assert(Dim == 2 || Dim == 3, + "to_mdspan only supports 2D or 3D uniform grids"); + } +} + +template +View UniformGridField::to_mdspan() const +{ + PCMS_FUNCTION_TIMER; + + if constexpr (Dim == 2) { + return View( + dof_holder_data_.data(), grid_.divisions[0] + 1, + grid_.divisions[1] + 1); + } else if constexpr (Dim == 3) { + return View( + dof_holder_data_.data(), grid_.divisions[0] + 1, + grid_.divisions[1] + 1, grid_.divisions[2] + 1); + } else { + static_assert(Dim == 2 || Dim == 3, + "to_mdspan only supports 2D or 3D uniform grids"); + } +} + template LocalizationHint UniformGridField::GetLocalizationHint( CoordinateView coordinate_view) const diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.h b/src/pcms/adapter/uniform_grid/uniform_grid_field.h index 49006b0d..eba81048 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field.h +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.h @@ -39,6 +39,9 @@ class UniformGridField : public FieldT Rank1View GetDOFHolderData() const override; void SetDOFHolderData(Rank1View data) override; + View to_mdspan(); + View to_mdspan() const; + ~UniformGridField() noexcept = default; private: diff --git a/src/pcms/create_field.cpp b/src/pcms/create_field.cpp index f9fb3809..c864ddcd 100644 --- a/src/pcms/create_field.cpp +++ b/src/pcms/create_field.cpp @@ -1,6 +1,8 @@ #include "create_field.h" #include "adapter/omega_h/omega_h_field2.h" #include "adapter/omega_h/omega_h_field_layout.h" +#include "adapter/uniform_grid/uniform_grid_field.h" +#include "adapter/uniform_grid/uniform_grid_field_layout.h" #include "uniform_grid.h" #include "point_search.h" @@ -28,53 +30,74 @@ std::unique_ptr CreateLagrangeLayout( } template <> -std::vector CreateUniformGridBinaryFieldFromGrid<2>( - Omega_h::Mesh& mesh, const UniformGrid<2>& grid) +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryFieldFromGrid<2>(Omega_h::Mesh& mesh, + UniformGrid<2>& grid) { constexpr unsigned dim = 2; - // Get total number of cells - const LO num_cells = grid.GetNumCells(); - - // Initialize result vector - std::vector binary_field(num_cells, 0); + // Get total number of vertices + const LO num_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1); // Create GridPointSearch for point-in-mesh queries GridPointSearch point_search(mesh, grid.divisions[0], grid.divisions[1]); - // Create array of grid cell center points - Kokkos::View cell_centers("cell centers", num_cells); - auto cell_centers_h = Kokkos::create_mirror_view(cell_centers); - - // Fill cell center coordinates - for (LO i = 0; i < num_cells; ++i) { - auto bbox = grid.GetCellBBOX(i); - for (unsigned d = 0; d < dim; ++d) { - cell_centers_h(i, d) = bbox.center[d]; + // Create array of grid vertex points + Kokkos::View vertices("vertices", num_vertices); + auto vertices_h = Kokkos::create_mirror_view(vertices); + + // Fill vertex coordinates + const Real dx = grid.edge_length[0] / grid.divisions[0]; + const Real dy = grid.edge_length[1] / grid.divisions[1]; + + for (LO j = 0; j <= grid.divisions[1]; ++j) { + for (LO i = 0; i <= grid.divisions[0]; ++i) { + const LO vertex_id = j * (grid.divisions[0] + 1) + i; + vertices_h(vertex_id, 0) = grid.bot_left[0] + i * dx; + vertices_h(vertex_id, 1) = grid.bot_left[1] + j * dy; } } // Copy to device - Kokkos::deep_copy(cell_centers, cell_centers_h); + Kokkos::deep_copy(vertices, vertices_h); // Perform point search - auto results = point_search(cell_centers); + auto results = point_search(vertices); + fprintf(stderr, "Completed point-in-mesh search for %d vertices.\n", + num_vertices); // Copy results back to host auto results_h = Kokkos::create_mirror_view(results); Kokkos::deep_copy(results_h, results); + fprintf(stderr, "Copied point search results back to host.\n"); - // Process results: set field to 1 if point is inside mesh (tri_id >= 0) - for (LO i = 0; i < num_cells; ++i) { - binary_field[i] = (results_h(i).tri_id >= 0) ? 1 : 0; + // Create UniformGridFieldLayout and field + auto layout = std::make_unique>( + grid, 1, CoordinateSystem::Cartesian); + auto field = std::make_unique>(*layout); + + // Create binary data as Real values (0.0 or 1.0) + Kokkos::View binary_data("binary_data", num_vertices); + for (LO i = 0; i < num_vertices; ++i) { + binary_data(i) = (results_h(i).tri_id >= 0) ? 1.0 : 0.0; } + fprintf(stderr, "Generated binary inside/outside data for grid vertices.\n"); + + // Set the DOF holder data + Rank1View data_view( + binary_data.data(), binary_data.extent(0)); + field->SetDOFHolderData(data_view); + fprintf(stderr, "Set binary data on UniformGridField.\n"); - return binary_field; + return {std::move(layout), std::move(field)}; } template <> -std::vector CreateUniformGridBinaryField<2>( - Omega_h::Mesh& mesh, const std::array& divisions) +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh, + const std::array& divisions) { constexpr unsigned dim = 2; @@ -86,8 +109,9 @@ std::vector CreateUniformGridBinaryField<2>( } template <> -std::vector CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh, - LO cells_per_dim) +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField<2>(Omega_h::Mesh& mesh, LO cells_per_dim) { std::array divisions; divisions.fill(cells_per_dim); diff --git a/src/pcms/create_field.h b/src/pcms/create_field.h index b595d1bf..0cc8b993 100644 --- a/src/pcms/create_field.h +++ b/src/pcms/create_field.h @@ -2,6 +2,8 @@ #define CREATE_FIELD_H_ #include "adapter/omega_h/omega_h_field_layout.h" +#include "adapter/uniform_grid/uniform_grid_field.h" +#include "adapter/uniform_grid/uniform_grid_field_layout.h" #include "field_layout.h" #include "field.h" #include "coordinate_system.h" @@ -23,9 +25,9 @@ std::unique_ptr CreateLagrangeLayout( * \brief Create a binary field on a uniform grid indicating inside/outside mesh * * This function creates a uniform grid from the mesh and generates a binary - * field where each grid cell is assigned: - * - 1 if the cell center is inside the original mesh - * - 0 if the cell center is outside the original mesh + * field where each grid vertex is assigned: + * - 1 if the vertex is inside the original mesh + * - 0 if the vertex is outside the original mesh * * Uses GridPointSearch to determine if a point lies within the mesh domain. * Currently only supports 2D meshes. @@ -33,14 +35,17 @@ std::unique_ptr CreateLagrangeLayout( * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported) * \param mesh The Omega_h mesh to create the grid from * \param divisions Array specifying the number of cells in each dimension - * \return std::vector Binary field values (0 or 1) for each grid cell + * \return std::pair containing the layout and field (layout must outlive field) * - * \note The returned vector has length equal to the number of grid cells. - * Cell values are ordered according to the grid's internal indexing. + * \note The returned field has vertex-centered data with values 0.0 or 1.0. + * Vertex values are ordered according to the grid's internal indexing. + * The layout must be kept alive as long as the field is used. */ template -std::vector CreateUniformGridBinaryField( - Omega_h::Mesh& mesh, const std::array& divisions); +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField(Omega_h::Mesh& mesh, + const std::array& divisions); /** * \brief Create a binary field with equal divisions in all dimensions @@ -51,19 +56,20 @@ std::vector CreateUniformGridBinaryField( * \tparam dim Spatial dimension of the mesh (currently only dim=2 is supported) * \param mesh The Omega_h mesh to create the grid from * \param cells_per_dim Number of cells per dimension (same for all dimensions) - * \return std::vector Binary field values (0 or 1) for each grid cell + * \return std::pair containing the layout and field (layout must outlive field) */ template -std::vector CreateUniformGridBinaryField(Omega_h::Mesh& mesh, - LO cells_per_dim); +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField(Omega_h::Mesh& mesh, LO cells_per_dim); /** * \brief Create a binary field on a given uniform grid * * This function takes a pre-defined uniform grid and generates a binary field - * where each grid cell is assigned: - * - 1 if the cell center is inside the mesh - * - 0 if the cell center is outside the mesh + * where each grid vertex is assigned: + * - 1 if the vertex is inside the mesh + * - 0 if the vertex is outside the mesh * * This allows testing with custom grids that may extend beyond the mesh * boundaries. @@ -71,11 +77,13 @@ std::vector CreateUniformGridBinaryField(Omega_h::Mesh& mesh, * \tparam dim Spatial dimension (currently only dim=2 is supported) * \param mesh The Omega_h mesh to test against * \param grid The uniform grid to evaluate - * \return std::vector Binary field values (0 or 1) for each grid cell + * \return std::pair containing the layout and field (layout must outlive field) */ template -std::vector CreateUniformGridBinaryFieldFromGrid( - Omega_h::Mesh& mesh, const UniformGrid& grid); +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryFieldFromGrid(Omega_h::Mesh& mesh, + UniformGrid& grid); } // namespace pcms diff --git a/src/pcms/pythonapi/CMakeLists.txt b/src/pcms/pythonapi/CMakeLists.txt new file mode 100644 index 00000000..c33b2c70 --- /dev/null +++ b/src/pcms/pythonapi/CMakeLists.txt @@ -0,0 +1,71 @@ +find_package(Python COMPONENTS Interpreter Development REQUIRED) +find_package(pybind11 REQUIRED) +include(FindPythonModule.cmake) +message(STATUS "Python_EXECUTABLE ${Python_EXECUTABLE}") +find_python_module(mpi4py REQUIRED) + +if(mpi4py_FOUND) + execute_process( + COMMAND + "${Python_EXECUTABLE}" "-c" + "import mpi4py as m; print(m.__version__); print(m.get_include());" + RESULT_VARIABLE + _mpi4py_SEARCH_SUCCESS + OUTPUT_VARIABLE + _mpi4py_VALUES + ERROR_VARIABLE + _mpi4py_ERROR_VALUE + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + + # Convert the process output into a list + string(REGEX REPLACE ";" "\\\\;" _mpi4py_VALUES ${_mpi4py_VALUES}) + string(REGEX REPLACE "\n" ";" _mpi4py_VALUES ${_mpi4py_VALUES}) + list(GET _mpi4py_VALUES 0 mpi4py_VERSION) + list(GET _mpi4py_VALUES 1 mpi4py_INCLUDE_DIRS) + + # Make sure all directory separators are '/' + string(REGEX REPLACE "\\\\" "/" mpi4py_INCLUDE_DIRS ${mpi4py_INCLUDE_DIRS}) + + # Get the major and minor version numbers + string(REGEX REPLACE "\\." ";" _mpi4py_VERSION_LIST ${mpi4py_VERSION}) + list(GET _mpi4py_VERSION_LIST 0 mpi4py_VERSION_MAJOR) + list(GET _mpi4py_VERSION_LIST 1 mpi4py_VERSION_MINOR) + list(GET _mpi4py_VERSION_LIST 2 mpi4py_VERSION_PATCH) + string(REGEX MATCH "[0-9]*" mpi4py_VERSION_PATCH ${mpi4py_VERSION_PATCH}) + math(EXPR mpi4py_VERSION_DECIMAL + "(${mpi4py_VERSION_MAJOR} * 10000) + (${mpi4py_VERSION_MINOR} * 100) + ${mpi4py_VERSION_PATCH}") +endif() + +pybind11_add_module(py_pcms + pythonapi.cpp + bind_field_base.cpp + bind_field_layout.cpp + bind_omega_h_field2.cpp + bind_omega_h_field_layout.cpp + bind_uniform_grid_field.cpp + bind_uniform_grid_field_layout.cpp + bind_omega_h.cpp + bind_transfer_field2.cpp +) +target_include_directories(py_pcms + SYSTEM PRIVATE + ${mpi4py_INCLUDE_DIRS} +) +set_target_properties(py_pcms PROPERTIES CXX_STANDARD 17) +target_link_libraries(py_pcms PRIVATE Kokkos::kokkos pybind11::module MPI::MPI_C pcms::core pcms::interpolator) + + + +install( + TARGETS py_pcms + EXPORT pcms_python-targets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR} + PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pcms/pythonapi/) +install( + EXPORT pcms_python-targets + NAMESPACE pcms:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms) \ No newline at end of file diff --git a/src/pcms/pythonapi/FindPythonModule.cmake b/src/pcms/pythonapi/FindPythonModule.cmake new file mode 100644 index 00000000..5583442b --- /dev/null +++ b/src/pcms/pythonapi/FindPythonModule.cmake @@ -0,0 +1,82 @@ +# Find if a Python module is installed. +# Usage: find_python_module( [[ATLEAST | EXACT] version] [QUIET] [REQUIRED]) + +macro(find_python_module module) + cmake_parse_arguments(ARG + "QUIET;REQUIRED" # options + "ATLEAST;EXACT" # one-value arguments + "" # multi-value arguments + ${ARGN} # everything else + ) + + if(ARG_QUIET) + set(${module}_FIND_QUIETLY TRUE) + endif() + + if(ARG_REQUIRED) + set(${module}_FIND_REQUIRED TRUE) + endif() + + if(ARG_ATLEAST AND ARG_EXACT) + message(FATAL_ERROR "Can't be both ATLEAST and EXACT") + endif() + if(ARG_ATLEAST) + set(_op ">=") + set(${module}_tgtver ${ARG_ATLEAST}) + elseif(ARG_EXACT) + set(_op "==") + set(${module}_tgtver ${ARG_EXACT}) + else() + # deceive handle_standard_arguments into not caring about version + set(_${module}_requested_version_found "${Python_EXECUTABLE}") + endif() + + unset(PY_${module} CACHE) + unset(${module}_VERSION CACHE) + + execute_process( + COMMAND "${Python_EXECUTABLE}" "-c" + "import re; \ + import ${module}; \ + print(re.compile('/__init__.py.*').sub('', ${module}.__file__))" + RESULT_VARIABLE _${module}_status + OUTPUT_VARIABLE _${module}_location + ERROR_QUIET + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + if(NOT ${_${module}_status}) + set(PY_${module} ${_${module}_location} CACHE STRING "Location of Python module ${module}" FORCE) + execute_process( + COMMAND "${Python_EXECUTABLE}" "-c" + "import sys; \ + import ${module}; \ + print(${module}.__version__)" + RESULT_VARIABLE _${module}_ver_status + OUTPUT_VARIABLE _${module}_version + ERROR_QUIET + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + if(NOT ${_${module}_ver_status}) + set(${module}_VERSION ${_${module}_version} CACHE STRING + "Version of Python module ${module}" FORCE) + + if(${module}_tgtver) + execute_process( + COMMAND "${Python_EXECUTABLE}" "-c" + "from pkg_resources import parse_version; \ + print(parse_version('${${module}_VERSION}') ${_op} parse_version('${${module}_tgtver}'))" + RESULT_VARIABLE _${module}_verenuf_status + OUTPUT_VARIABLE _${module}_verenuf + ERROR_QUIET + OUTPUT_STRIP_TRAILING_WHITESPACE + ) + if(NOT ${_${module}_verenuf_status}) + if(${_${module}_verenuf} STREQUAL "True") + set(_${module}_requested_version_found "${Python_EXECUTABLE}") + endif() + endif() + endif() + endif() + endif() + find_package_handle_standard_args(${module} DEFAULT_MSG PY_${module} _${module}_requested_version_found) +endmacro() \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp new file mode 100644 index 00000000..d98e0fb3 --- /dev/null +++ b/src/pcms/pythonapi/bind_field_base.cpp @@ -0,0 +1,282 @@ +#include +#include +#include +#include "pcms/coordinate_system.h" +#include "pcms/coordinate.h" +#include "pcms/create_field.h" +#include "pcms/uniform_grid.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_coordinate_system_module(py::module& m) { + // Bind CoordinateSystem enum + py::enum_(m, "CoordinateSystem") + .value("Cartesian", CoordinateSystem::Cartesian) + .value("Cylindrical", CoordinateSystem::Cylindrical) + .value("XGC", CoordinateSystem::XGC) + .value("GNET", CoordinateSystem::GNET) + .value("BEAMS3D", CoordinateSystem::BEAMS3D) + .export_values(); + + // Bind CoordinateView for HostMemorySpace + py::class_>(m, "CoordinateView") + .def(py::init([](CoordinateSystem cs, py::array_t coords) { + // Convert 2D numpy array to Rank2View + auto buf = coords.request(); + if (buf.ndim != 2) { + throw std::runtime_error("Coordinates must be a 2D array"); + } + // Create a view from the numpy array + Rank2View coords_view( + static_cast(buf.ptr), buf.shape[0], buf.shape[1]); + return CoordinateView(cs, coords_view); + }), + py::arg("coordinate_system"), + py::arg("coordinates"), + "Constructor for CoordinateView") + + .def("get_coordinate_system", &CoordinateView::GetCoordinateSystem, + "Get the coordinate system") + + .def("set_coordinate_system", &CoordinateView::SetCoordinateSystem, + py::arg("cs"), + "Set the coordinate system") + + .def("get_coordinates", [](const CoordinateView& self) { + auto coords = self.GetCoordinates(); + // Convert to numpy array + py::array_t result({static_cast(coords.extent(0)), + static_cast(coords.extent(1))}); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < coords.extent(0); ++i) { + for (size_t j = 0; j < coords.extent(1); ++j) { + ptr[i * coords.extent(1) + j] = coords(i, j); + } + } + return result; + }, + "Get the coordinates as numpy array"); + + // Bind CoordinateTransformation abstract base class + py::class_>( + m, "CoordinateTransformation") + .def("evaluate", &CoordinateTransformation::Evaluate, + py::arg("from"), + py::arg("to"), + "Evaluate the coordinate transformation"); +} + +void bind_coordinate_module(py::module& m) { + // Bind Coordinate template for common cases + // 3D Cartesian coordinates with Real type + py::class_>(m, "Coordinate3D") + .def(py::init(), + py::arg("x"), + py::arg("y"), + py::arg("z"), + "Constructor for 3D coordinate") + + .def("values", [](const Coordinate& self) { + auto vals = self.Values(); + py::array_t result(3); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + ptr[0] = vals[0]; + ptr[1] = vals[1]; + ptr[2] = vals[2]; + return result; + }, + "Get coordinate values as numpy array") + + .def("__getitem__", &Coordinate::operator[], + py::arg("index"), + "Get coordinate value by index"); + + // 2D coordinates + py::class_>(m, "Coordinate2D") + .def(py::init(), + py::arg("x"), + py::arg("y"), + "Constructor for 2D coordinate") + + .def("values", [](const Coordinate& self) { + auto vals = self.Values(); + py::array_t result(2); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + ptr[0] = vals[0]; + ptr[1] = vals[1]; + return result; + }, + "Get coordinate values as numpy array") + + .def("__getitem__", &Coordinate::operator[], + py::arg("index"), + "Get coordinate value by index"); + + // Bind CoordinateElement for common types + py::class_>(m, "CoordinateElement") + .def(py::init(), + py::arg("data"), + "Constructor for CoordinateElement") + + .def("underlying", + py::overload_cast<>(&CoordinateElement::underlying, py::const_), + "Get the underlying value"); +} + +void bind_create_field_module(py::module& m) { + // Bind CreateLagrangeLayout function with shared_ptr wrapper + // pybind11 handles shared_ptr better than unique_ptr for Python ownership + m.def("create_lagrange_layout", + [](Omega_h::Mesh& mesh, int order, int num_components, + CoordinateSystem coordinate_system) { + return std::shared_ptr( + CreateLagrangeLayout(mesh, order, num_components, coordinate_system)); + }, + py::arg("mesh"), + py::arg("order"), + py::arg("num_components"), + py::arg("coordinate_system")); + + // Bind CreateUniformGridFromMesh for 2D + m.def("create_uniform_grid_from_mesh", + [](Omega_h::Mesh& mesh, const std::array& divisions) { + return CreateUniformGridFromMesh<2>(mesh, divisions); + }, + py::arg("mesh"), + py::arg("divisions"), + "Create a 2D uniform grid from an Omega_h mesh"); + + // Bind CreateUniformGridFromMesh for 3D + m.def("create_uniform_grid_from_mesh", + [](Omega_h::Mesh& mesh, const std::array& divisions) { + return CreateUniformGridFromMesh<3>(mesh, divisions); + }, + py::arg("mesh"), + py::arg("divisions"), + "Create a 3D uniform grid from an Omega_h mesh"); + + // Bind CreateUniformGridBinaryField for 2D + m.def("create_uniform_grid_binary_field", + [](Omega_h::Mesh& mesh, const std::array& divisions) { + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, divisions); + // Wrap in shared_ptr for proper Python ownership and lifetime management + return py::make_tuple(std::shared_ptr>(std::move(layout)), + std::shared_ptr>(std::move(field))); + }, + py::arg("mesh"), + py::arg("divisions"), + "Create a 2D binary field on a uniform grid indicating inside/outside mesh. " + "Returns tuple of (layout, field). Layout lifetime is properly managed via shared_ptr."); +} + +template +void bind_field_t(py::module& m, const std::string& type_suffix) { + std::string class_name = "FieldT_" + type_suffix; + + py::class_, std::shared_ptr>>(m, class_name.c_str()) + .def("get_coordinate_system", &FieldT::GetCoordinateSystem, + "Get the coordinate system of the field") + + .def("get_localization_hint", &FieldT::GetLocalizationHint, + py::arg("coordinates"), + "Get a localization hint for a set of coordinates") + + .def("evaluate", [](const FieldT& self, + LocalizationHint hint, + py::array_t results_array, + CoordinateSystem coord_sys) { + auto results_view = numpy_to_view(results_array); + FieldDataView results(results_view, coord_sys); + self.Evaluate(hint, results); + }, + py::arg("hint"), + py::arg("results"), + py::arg("coordinate_system"), + "Evaluate the field at given locations") + + .def("evaluate_gradient", [](FieldT& self, + py::array_t results_array, + CoordinateSystem coord_sys) { + auto results_view = numpy_to_view(results_array); + FieldDataView results(results_view, coord_sys); + self.EvaluateGradient(results); + }, + py::arg("results"), + py::arg("coordinate_system"), + "Evaluate the gradient of the field") + + .def("get_dof_holder_data", [](const FieldT& self) { + auto data = self.GetDOFHolderData(); + return view_to_numpy(data); + }, + "Get the DOF holder data") + + .def("set_dof_holder_data", [](FieldT& self, py::array_t data) { + auto data_view = numpy_to_view(data); + self.SetDOFHolderData(data_view); + }, + py::arg("data"), + "Set the DOF holder data") + + .def("get_layout", &FieldT::GetLayout, + py::return_value_policy::reference, + "Get the field layout") + + .def("can_evaluate_gradient", &FieldT::CanEvaluateGradient, + "Check if the field can evaluate gradients") + + .def("serialize", [](const FieldT& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + return self.Serialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Serialize the field data") + + .def("deserialize", [](FieldT& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + self.Deserialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Deserialize field data from buffer"); +} + +void bind_field_module(py::module& m) { + // Bind LocalizationHint (opaque type - data member is internal only) + py::class_(m, "LocalizationHint") + .def(py::init<>()); + + // Bind FieldDataView for common types + py::class_>(m, "FieldDataView_Real") + .def(py::init, CoordinateSystem>(), + py::arg("values"), + py::arg("coordinate_system")) + .def("size", &FieldDataView::Size) + .def("get_coordinate_system", &FieldDataView::GetCoordinateSystem) + // TODO: const GetValues? + .def("get_values", [](FieldDataView& self) { + return view_to_numpy(self.GetValues()); + }); + + // Bind FieldT for common types + // bind_field_t(m, "Real"); + // bind_field_t(m, "Float"); + + // Only bind double separately if Real is not already double + bind_field_t(m, "Double"); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_field_layout.cpp b/src/pcms/pythonapi/bind_field_layout.cpp new file mode 100644 index 00000000..0f7c925b --- /dev/null +++ b/src/pcms/pythonapi/bind_field_layout.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include "pcms/field_layout.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_field_layout_module(py::module& m) { + // Bind the base FieldLayout class as an abstract base + py::class_>(m, "FieldLayout") + .def("create_field", &FieldLayout::CreateField, + "Create a field with this layout") + + .def("get_num_components", &FieldLayout::GetNumComponents, + "Get the number of components in the field") + + .def("get_num_owned_dof_holder", &FieldLayout::GetNumOwnedDofHolder, + "Get the number of owned DOF holders") + + .def("get_num_global_dof_holder", &FieldLayout::GetNumGlobalDofHolder, + "Get the number of global DOF holders") + + .def("owned_size", &FieldLayout::OwnedSize, + "Get the owned size (num_components * num_owned_dof_holder)") + + .def("global_size", &FieldLayout::GlobalSize, + "Get the global size (num_components * num_global_dof_holder)") + + .def("get_owned", [](const FieldLayout& self) { + auto owned = self.GetOwned(); + // Convert to numpy array + py::array_t result(owned.extent(0)); + auto buf = result.request(); + bool* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < owned.extent(0); ++i) { + ptr[i] = owned[i]; + } + return result; + }, + "Get the owned mask array") + + .def("get_gids", [](const FieldLayout& self) { + auto gids = self.GetGids(); + // Convert to numpy array + py::array_t result(gids.extent(0)); + auto buf = result.request(); + GO* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < gids.extent(0); ++i) { + ptr[i] = gids[i]; + } + return result; + }, + "Get the global IDs array") + + .def("is_distributed", &FieldLayout::IsDistributed, + "Check if the field layout is distributed") + + .def("get_ent_offsets", [](const FieldLayout& self) { + auto offsets = self.GetEntOffsets(); + // Convert std::array to list + py::list result; + for (size_t i = 0; i < offsets.size(); ++i) { + result.append(offsets[i]); + } + return result; + }, + "Get the entity offsets array") + + .def("get_dof_holder_coordinates", [](const FieldLayout& self) { + auto coords = self.GetDOFHolderCoordinates().GetCoordinates(); + // Convert to numpy array (2D) + py::array_t result({coords.extent(0), coords.extent(1)}); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < coords.extent(0); ++i) { + for (size_t j = 0; j < coords.extent(1); ++j) { + ptr[i * coords.extent(1) + j] = coords(i, j); + } + } + return result; + }, + "Get the DOF holder coordinates"); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_omega_h.cpp b/src/pcms/pythonapi/bind_omega_h.cpp new file mode 100644 index 00000000..eaed81c8 --- /dev/null +++ b/src/pcms/pythonapi/bind_omega_h.cpp @@ -0,0 +1,526 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef OMEGA_H_USE_ADIOS2 +#include +#endif +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_omega_h_mesh_module(py::module& m) { + // Bind Omega_h::Library + py::class_>(m, "OmegaHLibrary") + .def(py::init<>(), "Default constructor") + .def("world", &Omega_h::Library::world, + py::return_value_policy::reference, + "Get the world communicator"); + + // Bind the Comm class + py::class_>(m, "Comm") + // Constructors +#ifdef OMEGA_H_USE_MPI + .def(py::init(), + py::arg("library"), + py::arg("impl")) + .def(py::init([](Omega_h::Library* library, MPI_Comm impl, py::array_t srcs, + py::array_t dsts) { + auto srcs_view = numpy_to_omega_h_read(srcs); + auto dsts_view = numpy_to_omega_h_read(dsts); + return new Omega_h::Comm(library, impl, srcs_view, dsts_view); + })) + .def("get_impl", &Omega_h::Comm::get_impl, + "Get the underlying MPI communicator") +#else + .def(py::init(), + py::arg("library"), + py::arg("is_graph"), + py::arg("sends_to_self")) +#endif + // Methods + .def("library", &Omega_h::Comm::library, + py::return_value_policy::reference, + "Get the library pointer") + .def("rank", &Omega_h::Comm::rank, + "Get the rank of this process") + .def("size", &Omega_h::Comm::size, + "Get the total number of processes") + .def("dup", &Omega_h::Comm::dup, + "Duplicate the communicator") + .def("split", &Omega_h::Comm::split, + py::arg("color"), + py::arg("key"), + "Split the communicator") + .def("graph", &Omega_h::Comm::graph, + py::arg("dsts"), + "Create a graph communicator") + .def("graph_adjacent", &Omega_h::Comm::graph_adjacent, + py::arg("srcs"), + py::arg("dsts"), + "Create an adjacent graph communicator") + .def("graph_inverse", &Omega_h::Comm::graph_inverse, + "Get the inverse graph communicator") + .def("sources", &Omega_h::Comm::sources, + "Get source ranks") + .def("destinations", &Omega_h::Comm::destinations, + "Get destination ranks") + .def("reduce_or", &Omega_h::Comm::reduce_or, + py::arg("x"), + "Reduce using logical OR") + .def("reduce_and", &Omega_h::Comm::reduce_and, + py::arg("x"), + "Reduce using logical AND") + .def("add_int128", &Omega_h::Comm::add_int128, + py::arg("x"), + "Add Int128 values across processes") + .def("bcast_string", &Omega_h::Comm::bcast_string, + py::arg("s"), + py::arg("root_rank") = 0, + "Broadcast a string") + .def("barrier", &Omega_h::Comm::barrier, + "Synchronize all processes"); + m.def("build_box", &Omega_h::build_box, + py::arg("comm"), + py::arg("family"), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("nx"), + py::arg("ny"), + py::arg("nz"), + py::arg("symmetric") = false, + "Build a box mesh with specified dimensions and discretization", + py::return_value_policy::move); + + // Bind Omega_h::Mesh + py::class_>(m, "OmegaHMesh") + .def(py::init<>(), "Default constructor") + .def(py::init(), py::arg("library"), + "Constructor with library") + + .def("set_library", &Omega_h::Mesh::set_library, + py::arg("library"), + "Set the library") + + .def("library", &Omega_h::Mesh::library, + py::return_value_policy::reference, + "Get the library") + + .def("set_comm", &Omega_h::Mesh::set_comm, + py::arg("comm"), + "Set the communicator") + + .def("comm", &Omega_h::Mesh::comm, + "Get the communicator") + + .def("set_dim", &Omega_h::Mesh::set_dim, + py::arg("dim"), + "Set mesh dimension") + + .def("dim", &Omega_h::Mesh::dim, + "Get mesh dimension") + + .def("set_family", &Omega_h::Mesh::set_family, + py::arg("family"), + "Set mesh family (simplex, hypercube, etc.)") + + .def("family", &Omega_h::Mesh::family, + "Get mesh family") + + .def("nverts", &Omega_h::Mesh::nverts, + "Get number of vertices") + + .def("nedges", &Omega_h::Mesh::nedges, + "Get number of edges") + + .def("nfaces", &Omega_h::Mesh::nfaces, + "Get number of faces") + + .def("nregions", &Omega_h::Mesh::nregions, + "Get number of regions") + + .def("nelems", &Omega_h::Mesh::nelems, + "Get number of elements") + + .def("nents", &Omega_h::Mesh::nents, + py::arg("ent_dim"), + "Get number of entities of given dimension") + + .def("nglobal_ents", &Omega_h::Mesh::nglobal_ents, + py::arg("ent_dim"), + "Get global number of entities") + + .def("coords", [](const Omega_h::Mesh& mesh) { + auto coords = mesh.coords(); + return omega_h_read_to_numpy(coords); + }, + "Get mesh coordinates as numpy array") + + .def("set_coords", [](Omega_h::Mesh& mesh, py::array_t coords) { + auto coords_write = numpy_to_omega_h_write(coords); + mesh.set_coords(Omega_h::Reals(coords_write)); + }, + py::arg("coords"), + "Set mesh coordinates from numpy array") + + .def("add_coords", [](Omega_h::Mesh& mesh, py::array_t coords) { + auto coords_write = numpy_to_omega_h_write(coords); + mesh.add_coords(Omega_h::Reals(coords_write)); + }, + py::arg("coords"), + "Add mesh coordinates from numpy array") + + .def("has_tag", &Omega_h::Mesh::has_tag, + py::arg("ent_dim"), + py::arg("name"), + "Check if mesh has a tag") + + .def("ntags", &Omega_h::Mesh::ntags, + py::arg("ent_dim"), + "Get number of tags for entity dimension") + + .def("remove_tag", &Omega_h::Mesh::remove_tag, + py::arg("ent_dim"), + py::arg("name"), + "Remove a tag") + + .def("has_ents", &Omega_h::Mesh::has_ents, + py::arg("ent_dim"), + "Check if mesh has entities of given dimension") + + .def("ask_lengths", [](Omega_h::Mesh& mesh) { + auto lengths = mesh.ask_lengths(); + return omega_h_read_to_numpy(lengths); + }, + "Get edge lengths") + + .def("ask_qualities", [](Omega_h::Mesh& mesh) { + auto qualities = mesh.ask_qualities(); + return omega_h_read_to_numpy(qualities); + }, + "Get element qualities") + + .def("min_quality", &Omega_h::Mesh::min_quality, + "Get minimum element quality") + + .def("max_length", &Omega_h::Mesh::max_length, + "Get maximum edge length") + + .def("balance", py::overload_cast(&Omega_h::Mesh::balance), + py::arg("predictive") = false, + "Balance the mesh across processors") + + .def("set_parting", + py::overload_cast(&Omega_h::Mesh::set_parting), + py::arg("parting"), + py::arg("verbose") = false, + "Set mesh partitioning") + + .def("parting", &Omega_h::Mesh::parting, + "Get mesh partitioning type") + + .def("nghost_layers", &Omega_h::Mesh::nghost_layers, + "Get number of ghost layers") + + .def("owned", [](Omega_h::Mesh& mesh, Omega_h::Int ent_dim) { + auto owned = mesh.owned(ent_dim); + return omega_h_read_to_numpy(owned); + }, + py::arg("ent_dim"), + "Get ownership flags for entities"); + + // Bind mesh I/O functions + + // General mesh file reader (detects format) + m.def("read_mesh_file", + [](const std::string& filepath, std::shared_ptr comm) { + return Omega_h::read_mesh_file(filepath, comm); + }, + py::arg("filepath"), + py::arg("comm"), + "Read mesh from file (auto-detects format)", + py::return_value_policy::move); + + // Binary format I/O + m.def("read_mesh_binary", + [](const std::string& filepath, Omega_h::Library* lib) { + Omega_h::Mesh mesh(lib); + Omega_h::binary::read(filepath, lib->world(), &mesh); + return mesh; + }, + py::arg("filepath"), + py::arg("library"), + "Read mesh from binary file", + py::return_value_policy::move); + + m.def("read_mesh_binary", + [](const std::string& filepath, std::shared_ptr comm) { + return Omega_h::binary::read(filepath, comm); + }, + py::arg("filepath"), + py::arg("comm"), + "Read mesh from binary file", + py::return_value_policy::move); + + m.def("write_mesh_binary", + [](const std::string& filepath, Omega_h::Mesh& mesh) { + Omega_h::binary::write(filepath, &mesh); + }, + py::arg("filepath"), + py::arg("mesh"), + "Write mesh to binary file"); + + // Gmsh format I/O + m.def("read_mesh_gmsh", + [](const std::string& filepath, std::shared_ptr comm) { + return Omega_h::gmsh::read(filepath, comm); + }, + py::arg("filepath"), + py::arg("comm"), + "Read mesh from Gmsh file", + py::return_value_policy::move); + + m.def("write_mesh_gmsh", + [](const std::string& filepath, Omega_h::Mesh& mesh) { + Omega_h::gmsh::write(filepath, &mesh); + }, + py::arg("filepath"), + py::arg("mesh"), + "Write mesh to Gmsh file"); + +#ifdef OMEGA_H_USE_GMSH + m.def("read_mesh_gmsh_parallel", + [](const std::string& filepath, std::shared_ptr comm) { + return Omega_h::gmsh::read_parallel(filepath, comm); + }, + py::arg("filepath"), + py::arg("comm"), + "Read parallel Gmsh mesh (MSH 4.1+)", + py::return_value_policy::move); + + m.def("write_mesh_gmsh_parallel", + [](const std::string& filepath, Omega_h::Mesh& mesh) { + Omega_h::gmsh::write_parallel(filepath, mesh); + }, + py::arg("filepath"), + py::arg("mesh"), + "Write parallel Gmsh mesh (MSH 4.1)"); +#endif + + // VTK format I/O + m.def("write_mesh_vtu", + [](const std::string& filepath, Omega_h::Mesh& mesh, bool compress) { + Omega_h::vtk::write_vtu(filepath, &mesh, compress); + }, + py::arg("filepath"), + py::arg("mesh"), + py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS, + "Write mesh to VTU file"); + + m.def("write_mesh_vtu", + [](const std::string& filepath, Omega_h::Mesh& mesh, Omega_h::Int cell_dim, bool compress) { + Omega_h::vtk::write_vtu(filepath, &mesh, cell_dim, compress); + }, + py::arg("filepath"), + py::arg("mesh"), + py::arg("cell_dim"), + py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS, + "Write mesh to VTU file with specified cell dimension"); + + m.def("write_mesh_parallel_vtk", + [](const std::string& filepath, Omega_h::Mesh& mesh, bool compress) { + Omega_h::vtk::write_parallel(filepath, &mesh, compress); + }, + py::arg("filepath"), + py::arg("mesh"), + py::arg("compress") = OMEGA_H_DEFAULT_COMPRESS, + "Write mesh to parallel VTK files"); + + m.def("read_mesh_parallel_vtk", + [](const std::string& pvtupath, std::shared_ptr comm) { + Omega_h::Mesh mesh; + Omega_h::vtk::read_parallel(pvtupath, comm, &mesh); + return mesh; + }, + py::arg("pvtupath"), + py::arg("comm"), + "Read parallel VTK mesh", + py::return_value_policy::move); + +#ifdef OMEGA_H_USE_SEACASEXODUS + // Exodus ClassifyWith enum + py::enum_(m, "ExodusClassifyWith") + .value("NODE_SETS", Omega_h::exodus::NODE_SETS) + .value("SIDE_SETS", Omega_h::exodus::SIDE_SETS) + .export_values(); + + // Exodus file operations + m.def("exodus_open", + [](const std::string& filepath, bool verbose) { + return Omega_h::exodus::open(filepath, verbose); + }, + py::arg("filepath"), + py::arg("verbose") = false, + "Open an Exodus file and return file handle"); + + m.def("exodus_close", + [](int exodus_file) { + Omega_h::exodus::close(exodus_file); + }, + py::arg("exodus_file"), + "Close an Exodus file"); + + m.def("exodus_get_num_time_steps", + [](int exodus_file) { + return Omega_h::exodus::get_num_time_steps(exodus_file); + }, + py::arg("exodus_file"), + "Get the number of time steps in an Exodus file"); + + m.def("read_mesh_exodus", + [](int exodus_file, Omega_h::Mesh& mesh, bool verbose, int classify_with) { + Omega_h::exodus::read_mesh(exodus_file, &mesh, verbose, classify_with); + }, + py::arg("exodus_file"), + py::arg("mesh"), + py::arg("verbose") = false, + py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS, + "Read mesh from an open Exodus file"); + + m.def("read_exodus_nodal_fields", + [](int exodus_file, Omega_h::Mesh& mesh, int time_step, + const std::string& prefix, const std::string& postfix, bool verbose) { + Omega_h::exodus::read_nodal_fields(exodus_file, &mesh, time_step, + prefix, postfix, verbose); + }, + py::arg("exodus_file"), + py::arg("mesh"), + py::arg("time_step"), + py::arg("prefix") = "", + py::arg("postfix") = "", + py::arg("verbose") = false, + "Read nodal fields from an open Exodus file at a specific time step"); + + m.def("read_exodus_element_fields", + [](int exodus_file, Omega_h::Mesh& mesh, int time_step, + const std::string& prefix, const std::string& postfix, bool verbose) { + Omega_h::exodus::read_element_fields(exodus_file, &mesh, time_step, + prefix, postfix, verbose); + }, + py::arg("exodus_file"), + py::arg("mesh"), + py::arg("time_step"), + py::arg("prefix") = "", + py::arg("postfix") = "", + py::arg("verbose") = false, + "Read element fields from an open Exodus file at a specific time step"); + + m.def("read_mesh_exodus_sliced", + [](const std::string& filepath, std::shared_ptr comm, + bool verbose, int classify_with, int time_step) { + return Omega_h::exodus::read_sliced(filepath, comm, verbose, + classify_with, time_step); + }, + py::arg("filepath"), + py::arg("comm"), + py::arg("verbose") = false, + py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS, + py::arg("time_step") = -1, + "Read sliced Exodus mesh in parallel", + py::return_value_policy::move); + + m.def("write_mesh_exodus", + [](const std::string& filepath, Omega_h::Mesh& mesh, bool verbose, int classify_with) { + Omega_h::exodus::write(filepath, &mesh, verbose, classify_with); + }, + py::arg("filepath"), + py::arg("mesh"), + py::arg("verbose") = false, + py::arg("classify_with") = Omega_h::exodus::NODE_SETS | Omega_h::exodus::SIDE_SETS, + "Write mesh to Exodus file"); +#endif + +#ifdef OMEGA_H_USE_LIBMESHB + // MESHB format I/O + m.def("read_mesh_meshb", + [](Omega_h::Mesh& mesh, const std::string& filepath) { + Omega_h::meshb::read(&mesh, filepath); + }, + py::arg("mesh"), + py::arg("filepath"), + "Read mesh from MESHB file"); + + m.def("write_mesh_meshb", + [](Omega_h::Mesh& mesh, const std::string& filepath, int version) { + Omega_h::meshb::write(&mesh, filepath, version); + }, + py::arg("mesh"), + py::arg("filepath"), + py::arg("version") = 2, + "Write mesh to MESHB file"); + + m.def("read_meshb_sol", + [](Omega_h::Mesh& mesh, const std::string& filepath, const std::string& sol_name) { + Omega_h::meshb::read_sol(&mesh, filepath, sol_name); + }, + py::arg("mesh"), + py::arg("filepath"), + py::arg("sol_name"), + "Read solution/resolution data from MESHB .sol file"); + + m.def("write_meshb_sol", + [](Omega_h::Mesh& mesh, const std::string& filepath, const std::string& sol_name, int version) { + Omega_h::meshb::write_sol(&mesh, filepath, sol_name, version); + }, + py::arg("mesh"), + py::arg("filepath"), + py::arg("sol_name"), + py::arg("version") = 2, + "Write solution/resolution data to MESHB .sol file"); +#endif + +#ifdef OMEGA_H_USE_ADIOS2 + // ADIOS2 format I/O + m.def("read_mesh_adios2", + [](const std::string& filepath, Omega_h::Library* lib, const std::string& prefix) { + return Omega_h::adios::read(filepath, lib, prefix); + }, + py::arg("filepath"), + py::arg("library"), + py::arg("prefix") = "", + "Read mesh from ADIOS2 file", + py::return_value_policy::move); + + m.def("write_mesh_adios2", + [](const std::string& filepath, Omega_h::Mesh& mesh, const std::string& prefix) { + Omega_h::adios::write(filepath, &mesh, prefix); + }, + py::arg("filepath"), + py::arg("mesh"), + py::arg("prefix") = "", + "Write mesh to ADIOS2 file"); +#endif + + // Other utilities + py::enum_(m, "Family") + .value("SIMPLEX", OMEGA_H_SIMPLEX) + .value("HYPERCUBE", OMEGA_H_HYPERCUBE) + .value("MIXED", OMEGA_H_MIXED) + .export_values(); // This allows using values without the class prefix + + // Bind Omega_h_Parting enum + py::enum_(m, "OmegaHParting") + .value("ELEM_BASED", OMEGA_H_ELEM_BASED) + .value("VERT_BASED", OMEGA_H_VERT_BASED) + .value("GHOSTED", OMEGA_H_GHOSTED) + .export_values(); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_omega_h_field2.cpp b/src/pcms/pythonapi/bind_omega_h_field2.cpp new file mode 100644 index 00000000..46d79fb1 --- /dev/null +++ b/src/pcms/pythonapi/bind_omega_h_field2.cpp @@ -0,0 +1,127 @@ +#include +#include +#include +#include +#include "pcms/arrays.h" +#include "pcms/adapter/omega_h/omega_h_field2.h" +#include "pcms/adapter/omega_h/omega_h_field_layout.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_omega_h_field2(py::module& m) { + + // Bind OmegaHField2 class + py::class_, std::shared_ptr>(m, "OmegaHField2") + .def(py::init(), + py::arg("layout"), + "Constructor for OmegaHField2") + + .def("get_localization_hint", [](const OmegaHField2& self, + py::array_t coordinates, + const CoordinateSystem& coord_system) { + // Create CoordinateView from numpy array + auto coords_view = numpy_to_view_2d(coordinates); + CoordinateView coord_view(coord_system, coords_view); + return self.GetLocalizationHint(coord_view); + }, + py::arg("coordinates"), + py::arg("coordinate_system"), + "Get localization hint for given coordinates") + + .def("evaluate", [](const OmegaHField2& self, + const LocalizationHint& location, + py::array_t values, + const CoordinateSystem& coord_system) { + // Create FieldDataView + auto values_view = numpy_to_view(values); + FieldDataView field_data_view(values_view, coord_system); + self.Evaluate(location, field_data_view); + }, + py::arg("location"), + py::arg("values"), + py::arg("coordinate_system"), + "Evaluate field at locations specified by localization hint") + + .def("evaluate_gradient", [](OmegaHField2& self, + py::array_t gradients, + const CoordinateSystem& coord_system) { + auto gradients_view = numpy_to_view(gradients); + FieldDataView field_data_view(gradients_view, coord_system); + self.EvaluateGradient(field_data_view); + }, + py::arg("gradients"), + py::arg("coordinate_system"), + "Evaluate gradient of the field") + + .def("get_layout", &OmegaHField2::GetLayout, + py::return_value_policy::reference, + "Get the field layout") + + .def("can_evaluate_gradient", &OmegaHField2::CanEvaluateGradient, + "Check if gradient evaluation is supported") + + .def("serialize", [](const OmegaHField2& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + return self.Serialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Serialize field data into buffer") + + .def("deserialize", [](OmegaHField2& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + self.Deserialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Deserialize field data from buffer") + + .def("get_dof_holder_data", [](const OmegaHField2& self) { + auto const_data = self.GetDOFHolderData(); + // Create a numpy array that owns its own data + return view_to_numpy(const_data); + }, + "Get the DOF holder data") + + .def("set_dof_holder_data", [](OmegaHField2& self, + py::array_t data) { + // Ensure array is contiguous + auto contiguous_data = py::array_t(data); + auto data_view = numpy_to_view(contiguous_data); + // Create const view wrapper + Rank1View const_view(data_view.data_handle(), data_view.size()); + self.SetDOFHolderData(const_view); + }, + py::arg("data"), + "Set the DOF holder data"); + + // Helper functions for creating views (if needed for testing) + m.def("create_coordinate_view", [](py::array_t coordinates, + const CoordinateSystem& coord_system) { + auto coords_view = numpy_to_view_2d(coordinates); + return CoordinateView(coord_system, coords_view); + }, + py::arg("coordinates"), + py::arg("coordinate_system"), + "Create a CoordinateView from numpy array"); + + m.def("create_field_data_view", [](py::array_t values, + const CoordinateSystem& coord_system) { + auto values_view = numpy_to_view(values); + return FieldDataView(values_view, coord_system); + }, + py::arg("values"), + py::arg("coordinate_system"), + "Create a FieldDataView from numpy array"); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_omega_h_field_layout.cpp b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp new file mode 100644 index 00000000..e228a43e --- /dev/null +++ b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp @@ -0,0 +1,119 @@ +#include +#include +#include +#include "pcms/adapter/omega_h/omega_h_field_layout.h" +#include "pcms/field.h" +#include "pcms/field_layout.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_omega_h_field_layout_module(py::module& m) { + // Bind the OmegaHFieldLayout class + py::class_>(m, "OmegaHFieldLayout") + .def(py::init, int, CoordinateSystem, std::string>(), + py::arg("mesh"), + py::arg("nodes_per_dim"), + py::arg("num_components"), + py::arg("coordinate_system"), + py::arg("global_id_name") = "global", + "Constructor for OmegaHFieldLayout") + + .def("create_field", [](OmegaHFieldLayout& self) { + return std::shared_ptr>(self.CreateField()); + }, + "Create a field with this layout") + + .def("get_num_components", &OmegaHFieldLayout::GetNumComponents, + "Get the number of components in the field") + + .def("get_num_owned_dof_holder", &OmegaHFieldLayout::GetNumOwnedDofHolder, + "Get the number of owned DOF holders") + + .def("get_num_global_dof_holder", &OmegaHFieldLayout::GetNumGlobalDofHolder, + "Get the number of global DOF holders") + + .def("get_owned", [](const OmegaHFieldLayout& self) { + auto owned = self.GetOwned(); + // Convert to numpy array + py::array_t result(owned.extent(0)); + auto buf = result.request(); + bool* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < owned.extent(0); ++i) { + ptr[i] = owned[i]; + } + return result; + }, + "Get the owned mask array") + + .def("get_gids", [](const OmegaHFieldLayout& self) { + auto gids = self.GetGids(); + // Convert to numpy array + py::array_t result(gids.extent(0)); + auto buf = result.request(); + GO* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < gids.extent(0); ++i) { + ptr[i] = gids[i]; + } + return result; + }, + "Get the global IDs array") + + .def("get_dof_holder_coordinates", [](const OmegaHFieldLayout& self) { + auto coords = self.GetDOFHolderCoordinates().GetCoordinates(); + // Convert to numpy array (2D) + py::array_t result({coords.extent(0), coords.extent(1)}); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < coords.extent(0); ++i) { + for (size_t j = 0; j < coords.extent(1); ++j) { + ptr[i * coords.extent(1) + j] = coords(i, j); + } + } + return result; + }, + "Get the DOF holder coordinates") + + .def("is_distributed", &OmegaHFieldLayout::IsDistributed, + "Check if the field layout is distributed") + + .def("get_ent_offsets", [](const OmegaHFieldLayout& self) { + auto offsets = self.GetEntOffsets(); + // Convert std::array to list/tuple + py::list result; + for (size_t i = 0; i < offsets.size(); ++i) { + result.append(offsets[i]); + } + return result; + }, + "Get the entity offsets array") + + .def("get_nodes_per_dim", [](const OmegaHFieldLayout& self) { + auto nodes = self.GetNodesPerDim(); + py::list result; + for (size_t i = 0; i < nodes.size(); ++i) { + result.append(nodes[i]); + } + return result; + }, + "Get the nodes per dimension array") + + .def("get_num_ents", &OmegaHFieldLayout::GetNumEnts, + "Get the total number of entities") + + .def("get_mesh", [](OmegaHFieldLayout& self) -> Omega_h::Mesh& { + return self.GetMesh(); + }, + py::return_value_policy::reference, + "Get the underlying Omega_h mesh") + + .def("owned_size", &OmegaHFieldLayout::OwnedSize, + "Get the owned size (num_components * num_owned_dof_holder)") + + .def("global_size", &OmegaHFieldLayout::GlobalSize, + "Get the global size (num_components * num_global_dof_holder)"); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_transfer_field2.cpp b/src/pcms/pythonapi/bind_transfer_field2.cpp new file mode 100644 index 00000000..3f9c28b7 --- /dev/null +++ b/src/pcms/pythonapi/bind_transfer_field2.cpp @@ -0,0 +1,62 @@ +#include +#include +#include "pcms/transfer_field2.h" +#include "pcms/field.h" + +namespace py = pybind11; + +namespace pcms { + +template +void bind_transfer_field2(py::module& m, const std::string& type_suffix) { + std::string copy_name = "copy_field2_" + type_suffix; + std::string interpolate_name = "interpolate_field2_" + type_suffix; + + m.def(copy_name.c_str(), + [](const FieldT& source, FieldT& target) { + copy_field2(source, target); + }, + py::arg("source"), + py::arg("target"), + "Copy field data from source to target. Source and target must be of the same type."); + + m.def(interpolate_name.c_str(), + [](const FieldT& source, FieldT& target) { + interpolate_field2(source, target); + }, + py::arg("source"), + py::arg("target"), + "Interpolate field from source to target. Coordinate systems must match."); +} + +void bind_transfer_field2_module(py::module& m) { + // Bind for Real type (typically double) + bind_transfer_field2(m, "Real"); + + // If Real is not double, also bind double explicitly + if constexpr (!std::is_same_v) { + bind_transfer_field2(m, "Double"); + } + + // Bind for float if needed + bind_transfer_field2(m, "Float"); + + // Add convenience aliases for the most common case (Real) + m.def("interpolate_field", + [](const FieldT& source, FieldT& target) { + interpolate_field2(source, target); + }, + py::arg("source"), + py::arg("target"), + "Interpolate field from source to target. Coordinate systems must match."); + + m.def("copy_field", + [](const FieldT& source, FieldT& target) { + copy_field2(source, target); + }, + py::arg("source"), + py::arg("target"), + "Copy field data from source to target. Source and target must be of the same type."); +} + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/bind_uniform_grid_field.cpp b/src/pcms/pythonapi/bind_uniform_grid_field.cpp new file mode 100644 index 00000000..3f9436f5 --- /dev/null +++ b/src/pcms/pythonapi/bind_uniform_grid_field.cpp @@ -0,0 +1,230 @@ +#include +#include +#include +#include +#include "pcms/arrays.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_uniform_grid_field_module(py::module& m) { + + // Bind UniformGridField class for 2D + py::class_, FieldT, std::shared_ptr>>( + m, "UniformGridField2D") + .def(py::init&>(), + py::arg("layout"), + "Constructor for UniformGridField2D") + + .def("get_localization_hint", [](const UniformGridField<2>& self, + py::array_t coordinates, + const CoordinateSystem& coord_system) { + // Create CoordinateView from numpy array + auto coords_view = numpy_to_view_2d(coordinates); + CoordinateView coord_view(coord_system, coords_view); + return self.GetLocalizationHint(coord_view); + }, + py::arg("coordinates"), + py::arg("coordinate_system"), + "Get localization hint for given coordinates") + + .def("evaluate", [](const UniformGridField<2>& self, + const LocalizationHint& location, + py::array_t values, + const CoordinateSystem& coord_system) { + // Create FieldDataView + auto values_view = numpy_to_view(values); + FieldDataView field_data_view(values_view, coord_system); + self.Evaluate(location, field_data_view); + }, + py::arg("location"), + py::arg("values"), + py::arg("coordinate_system"), + "Evaluate field at locations specified by localization hint") + + .def("evaluate_gradient", [](UniformGridField<2>& self, + py::array_t gradients, + const CoordinateSystem& coord_system) { + auto gradients_view = numpy_to_view(gradients); + FieldDataView field_data_view(gradients_view, coord_system); + self.EvaluateGradient(field_data_view); + }, + py::arg("gradients"), + py::arg("coordinate_system"), + "Evaluate gradient of the field") + + .def("get_layout", &UniformGridField<2>::GetLayout, + py::return_value_policy::reference, + "Get the field layout") + + .def("can_evaluate_gradient", &UniformGridField<2>::CanEvaluateGradient, + "Check if gradient evaluation is supported") + + .def("serialize", [](const UniformGridField<2>& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + return self.Serialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Serialize field data into buffer") + + .def("deserialize", [](UniformGridField<2>& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + self.Deserialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Deserialize field data from buffer") + + .def("get_dof_holder_data", [](const UniformGridField<2>& self) { + auto const_data = self.GetDOFHolderData(); + // Create a numpy array that owns its own data + return view_to_numpy(const_data); + }, + "Get the DOF holder data") + + .def("set_dof_holder_data", [](UniformGridField<2>& self, + py::array_t data) { + // Ensure array is contiguous + auto contiguous_data = py::array_t(data); + auto data_view = numpy_to_view(contiguous_data); + // Create const view wrapper + Rank1View const_view(data_view.data_handle(), data_view.size()); + self.SetDOFHolderData(const_view); + }, + py::arg("data"), + "Set the DOF holder data") + + .def("to_mdspan", [](const UniformGridField<2>& self) { + return mdspan_view_to_numpy(self.to_mdspan()); + }, + "Get the DOF holder data as a 2D numpy array (copy)"); + + // Bind UniformGridField class for 3D + py::class_, FieldT, std::shared_ptr>>( + m, "UniformGridField3D") + .def(py::init&>(), + py::arg("layout"), + "Constructor for UniformGridField3D") + + .def("get_localization_hint", [](const UniformGridField<3>& self, + py::array_t coordinates, + const CoordinateSystem& coord_system) { + // Create CoordinateView from numpy array + auto coords_view = numpy_to_view_2d(coordinates); + CoordinateView coord_view(coord_system, coords_view); + return self.GetLocalizationHint(coord_view); + }, + py::arg("coordinates"), + py::arg("coordinate_system"), + "Get localization hint for given coordinates") + + .def("evaluate", [](const UniformGridField<3>& self, + const LocalizationHint& location, + py::array_t values, + const CoordinateSystem& coord_system) { + // Create FieldDataView + auto values_view = numpy_to_view(values); + FieldDataView field_data_view(values_view, coord_system); + self.Evaluate(location, field_data_view); + }, + py::arg("location"), + py::arg("values"), + py::arg("coordinate_system"), + "Evaluate field at locations specified by localization hint") + + .def("evaluate_gradient", [](UniformGridField<3>& self, + py::array_t gradients, + const CoordinateSystem& coord_system) { + auto gradients_view = numpy_to_view(gradients); + FieldDataView field_data_view(gradients_view, coord_system); + self.EvaluateGradient(field_data_view); + }, + py::arg("gradients"), + py::arg("coordinate_system"), + "Evaluate gradient of the field") + + .def("get_layout", &UniformGridField<3>::GetLayout, + py::return_value_policy::reference, + "Get the field layout") + + .def("can_evaluate_gradient", &UniformGridField<3>::CanEvaluateGradient, + "Check if gradient evaluation is supported") + + .def("serialize", [](const UniformGridField<3>& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + return self.Serialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Serialize field data into buffer") + + .def("deserialize", [](UniformGridField<3>& self, + py::array_t buffer, + py::array_t permutation) { + auto buffer_view = numpy_to_view(buffer); + auto perm_view = numpy_to_view(permutation); + self.Deserialize(buffer_view, perm_view); + }, + py::arg("buffer"), + py::arg("permutation"), + "Deserialize field data from buffer") + + .def("get_dof_holder_data", [](const UniformGridField<3>& self) { + auto const_data = self.GetDOFHolderData(); + // Create a numpy array that owns its own data + return view_to_numpy(const_data); + }, + "Get the DOF holder data") + + .def("set_dof_holder_data", [](UniformGridField<3>& self, + py::array_t data) { + // Ensure array is contiguous + auto contiguous_data = py::array_t(data); + auto data_view = numpy_to_view(contiguous_data); + // Create const view wrapper + Rank1View const_view(data_view.data_handle(), data_view.size()); + self.SetDOFHolderData(const_view); + }, + py::arg("data"), + "Set the DOF holder data") + + .def("to_mdspan", [](const UniformGridField<3>& self) { + return mdspan_view_to_numpy(self.to_mdspan()); + }, + "Get the DOF holder data as a 3D numpy array (copy)"); + + // Helper functions for creating views (if needed for testing) + m.def("create_coordinate_view", [](py::array_t coordinates, + const CoordinateSystem& coord_system) { + auto coords_view = numpy_to_view_2d(coordinates); + return CoordinateView(coord_system, coords_view); + }, + py::arg("coordinates"), + py::arg("coordinate_system"), + "Create a CoordinateView from numpy array"); + + m.def("create_field_data_view", [](py::array_t values, + const CoordinateSystem& coord_system) { + auto values_view = numpy_to_view(values); + return FieldDataView(values_view, coord_system); + }, + py::arg("values"), + py::arg("coordinate_system"), + "Create a FieldDataView from numpy array"); +} + +} // namespace pcms diff --git a/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp new file mode 100644 index 00000000..2d8329ad --- /dev/null +++ b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp @@ -0,0 +1,222 @@ +#include +#include +#include +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/field.h" +#include "pcms/field_layout.h" +#include "pcms/uniform_grid.h" +#include "numpy_array_transform.h" + +namespace py = pybind11; + +namespace pcms { + +void bind_uniform_grid_field_layout_module(py::module& m) { + // Bind UniformGrid structure for 2D + py::class_>(m, "UniformGrid2D") + .def(py::init<>()) + .def_readwrite("edge_length", &UniformGrid<2>::edge_length, + "Edge length in each dimension") + .def_readwrite("bot_left", &UniformGrid<2>::bot_left, + "Bottom-left corner coordinates") + .def_readwrite("divisions", &UniformGrid<2>::divisions, + "Number of divisions in each dimension") + .def("get_num_cells", &UniformGrid<2>::GetNumCells, + "Get the total number of cells in the grid") + .def("closest_cell_id", [](const UniformGrid<2>& self, py::array_t point) { + if (point.size() != 2) { + throw std::runtime_error("Point must have 2 coordinates for 2D grid"); + } + auto buf = point.request(); + Real* ptr = static_cast(buf.ptr); + Omega_h::Vector<2> omega_point({ptr[0], ptr[1]}); + return self.ClosestCellID(omega_point); + }, + py::arg("point"), + "Get the cell ID that contains or is closest to the given point") + .def("get_cell_bbox", &UniformGrid<2>::GetCellBBOX, + py::arg("cell_index"), + "Get the bounding box of a cell"); + + // Bind UniformGrid structure for 3D + py::class_>(m, "UniformGrid3D") + .def(py::init<>()) + .def_readwrite("edge_length", &UniformGrid<3>::edge_length, + "Edge length in each dimension") + .def_readwrite("bot_left", &UniformGrid<3>::bot_left, + "Bottom-left corner coordinates") + .def_readwrite("divisions", &UniformGrid<3>::divisions, + "Number of divisions in each dimension") + .def("get_num_cells", &UniformGrid<3>::GetNumCells, + "Get the total number of cells in the grid") + .def("closest_cell_id", [](const UniformGrid<3>& self, py::array_t point) { + if (point.size() != 3) { + throw std::runtime_error("Point must have 3 coordinates for 3D grid"); + } + auto buf = point.request(); + Real* ptr = static_cast(buf.ptr); + Omega_h::Vector<3> omega_point({ptr[0], ptr[1], ptr[2]}); + return self.ClosestCellID(omega_point); + }, + py::arg("point"), + "Get the cell ID that contains or is closest to the given point") + .def("get_cell_bbox", &UniformGrid<3>::GetCellBBOX, + py::arg("cell_index"), + "Get the bounding box of a cell"); + + // Bind the UniformGridFieldLayout class for 2D + py::class_, FieldLayout, std::shared_ptr>>( + m, "UniformGridFieldLayout2D") + .def(py::init&, int, CoordinateSystem>(), + py::arg("grid"), + py::arg("num_components"), + py::arg("coordinate_system"), + "Constructor for UniformGridFieldLayout2D") + + .def("create_field", [](UniformGridFieldLayout<2>& self) { + return std::shared_ptr>(self.CreateField()); + }, + "Create a field with this layout") + + .def("get_num_components", &UniformGridFieldLayout<2>::GetNumComponents, + "Get the number of components in the field") + + .def("get_num_owned_dof_holder", &UniformGridFieldLayout<2>::GetNumOwnedDofHolder, + "Get the number of owned DOF holders") + + .def("get_num_global_dof_holder", &UniformGridFieldLayout<2>::GetNumGlobalDofHolder, + "Get the number of global DOF holders") + + .def("get_owned", [](const UniformGridFieldLayout<2>& self) { + auto owned = self.GetOwned(); + // Convert to numpy array + py::array_t result(owned.extent(0)); + auto buf = result.request(); + bool* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < owned.extent(0); ++i) { + ptr[i] = owned[i]; + } + return result; + }, + "Get the owned mask array") + + .def("get_gids", [](const UniformGridFieldLayout<2>& self) { + auto gids = self.GetGids(); + // Convert to numpy array + py::array_t result(gids.extent(0)); + auto buf = result.request(); + GO* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < gids.extent(0); ++i) { + ptr[i] = gids[i]; + } + return result; + }, + "Get the global IDs array") + + .def("get_dof_holder_coordinates", [](const UniformGridFieldLayout<2>& self) { + auto coords = self.GetDOFHolderCoordinates().GetCoordinates(); + // Convert to numpy array (2D) + py::array_t result({coords.extent(0), coords.extent(1)}); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < coords.extent(0); ++i) { + for (size_t j = 0; j < coords.extent(1); ++j) { + ptr[i * coords.extent(1) + j] = coords(i, j); + } + } + return result; + }, + "Get the DOF holder coordinates") + + .def("is_distributed", &UniformGridFieldLayout<2>::IsDistributed, + "Check if the field layout is distributed") + + .def("get_grid", &UniformGridFieldLayout<2>::GetGrid, + py::return_value_policy::reference, + "Get the underlying uniform grid") + + .def("get_num_cells", &UniformGridFieldLayout<2>::GetNumCells, + "Get the number of cells in the grid") + + .def("get_num_vertices", &UniformGridFieldLayout<2>::GetNumVertices, + "Get the number of vertices in the grid"); + + // Bind the UniformGridFieldLayout class for 3D + py::class_, FieldLayout, std::shared_ptr>>( + m, "UniformGridFieldLayout3D") + .def(py::init&, int, CoordinateSystem>(), + py::arg("grid"), + py::arg("num_components"), + py::arg("coordinate_system"), + "Constructor for UniformGridFieldLayout3D") + + .def("create_field", [](UniformGridFieldLayout<3>& self) { + return std::shared_ptr>(self.CreateField()); + }, + "Create a field with this layout") + + .def("get_num_components", &UniformGridFieldLayout<3>::GetNumComponents, + "Get the number of components in the field") + + .def("get_num_owned_dof_holder", &UniformGridFieldLayout<3>::GetNumOwnedDofHolder, + "Get the number of owned DOF holders") + + .def("get_num_global_dof_holder", &UniformGridFieldLayout<3>::GetNumGlobalDofHolder, + "Get the number of global DOF holders") + + .def("get_owned", [](const UniformGridFieldLayout<3>& self) { + auto owned = self.GetOwned(); + // Convert to numpy array + py::array_t result(owned.extent(0)); + auto buf = result.request(); + bool* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < owned.extent(0); ++i) { + ptr[i] = owned[i]; + } + return result; + }, + "Get the owned mask array") + + .def("get_gids", [](const UniformGridFieldLayout<3>& self) { + auto gids = self.GetGids(); + // Convert to numpy array + py::array_t result(gids.extent(0)); + auto buf = result.request(); + GO* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < gids.extent(0); ++i) { + ptr[i] = gids[i]; + } + return result; + }, + "Get the global IDs array") + + .def("get_dof_holder_coordinates", [](const UniformGridFieldLayout<3>& self) { + auto coords = self.GetDOFHolderCoordinates().GetCoordinates(); + // Convert to numpy array (2D) + py::array_t result({coords.extent(0), coords.extent(1)}); + auto buf = result.request(); + Real* ptr = static_cast(buf.ptr); + for (size_t i = 0; i < coords.extent(0); ++i) { + for (size_t j = 0; j < coords.extent(1); ++j) { + ptr[i * coords.extent(1) + j] = coords(i, j); + } + } + return result; + }, + "Get the DOF holder coordinates") + + .def("is_distributed", &UniformGridFieldLayout<3>::IsDistributed, + "Check if the field layout is distributed") + + .def("get_grid", &UniformGridFieldLayout<3>::GetGrid, + py::return_value_policy::reference, + "Get the underlying uniform grid") + + .def("get_num_cells", &UniformGridFieldLayout<3>::GetNumCells, + "Get the number of cells in the grid") + + .def("get_num_vertices", &UniformGridFieldLayout<3>::GetNumVertices, + "Get the number of vertices in the grid"); +} + +} // namespace pcms diff --git a/src/pcms/pythonapi/numpy_array_transform.h b/src/pcms/pythonapi/numpy_array_transform.h new file mode 100644 index 00000000..b1d6203d --- /dev/null +++ b/src/pcms/pythonapi/numpy_array_transform.h @@ -0,0 +1,232 @@ +#include +#include +#include +#include +#include "pcms/arrays.h" +#include "pcms/types.h" + + +namespace py = pybind11; + + +namespace pcms { + + +// Helper function to convert numpy array to Rank1 View +template +Rank1View numpy_to_view(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 1) { + throw std::runtime_error("Number of dimensions must be 1"); + } + T* data_ptr = reinterpret_cast(buf.ptr); + Rank1View view(data_ptr, buf.shape[0]); + return view; +} + + +// Helper function to convert numpy array to Kokkos View +template +Kokkos::View numpy_to_kokkos_view(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 1) { + throw std::runtime_error("Number of dimensions must be 1"); + } + Kokkos::View view( + reinterpret_cast(buf.ptr), buf.shape[0]); + return view; +} + + +// Helper function to convert Rank1 View to array (creates a copy) +template +py::array_t::type> view_to_numpy(Rank1View view) { + // Create a copy to avoid lifetime issues and mdspan type registration + using NonConstT = typename std::remove_const::type; + py::array_t result(view.size()); + py::buffer_info buf = result.request(); + NonConstT* ptr = static_cast(buf.ptr); + // Use element-wise copy to ensure proper mdspan access + for (size_t i = 0; i < view.size(); ++i) { + ptr[i] = view[i]; + } + return result; +} + + +// Helper function to convert Kokkos View to array (creates a reference) +template +py::array_t kokkos_view_to_numpy(Kokkos::View view) { + return py::array_t( + {static_cast(view.extent(0))}, // shape + {sizeof(T)}, // strides + view.data(), // data pointer + py::cast(view) // base object to manage lifetime + ); +} + + +// Helper to convert 2D numpy array to Kokkos View +template +Kokkos::View numpy_to_kokkos_view_2d(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 2) { + throw std::runtime_error("Number of dimensions must be 2"); + } + Kokkos::View view( + reinterpret_cast(buf.ptr), buf.shape[0], buf.shape[1]); + return view; +} + + +// Helper to convert 2D Kokkos View to numpy array (creates a reference) +template +py::array_t kokkos_view_2d_to_numpy(Kokkos::View view) { + return py::array_t( + {static_cast(view.extent(0)), static_cast(view.extent(1))}, // shape + {sizeof(T) * view.extent(1), sizeof(T)}, // strides (row-major) + view.data(), // data pointer + py::cast(view) // base object to manage lifetime + ); +} + + +// Helper to convert 1D numpy array to Omega_h::Read +template +Omega_h::Read numpy_to_omega_h_read(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 1) { + throw std::runtime_error("Number of dimensions must be 1"); + } + Kokkos::View> view( + reinterpret_cast(buf.ptr), buf.shape[0]); + Omega_h::Write write_view(view); + Omega_h::Read read_view(write_view); + return read_view; +} + + +// Helper to convert Omega_h::Read to numpy array (creates a copy) +template +py::array_t omega_h_read_to_numpy(Omega_h::Read read_view) { + auto read_view_host = Omega_h::HostRead(read_view); + py::array_t result(read_view_host.size()); + py::buffer_info buf = result.request(); + T* ptr = static_cast(buf.ptr); + for (Omega_h::LO i = 0; i < read_view_host.size(); ++i) { + ptr[i] = read_view_host[i]; + } + return result; +} + + +// Helper to convert 1D numpy array to Omega_h::Write +template +Omega_h::Write numpy_to_omega_h_write(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 1) { + throw std::runtime_error("Number of dimensions must be 1"); + } + // Get host mirror and copy data + auto write_view_host = Omega_h::HostWrite(buf.shape[0]); + T* ptr = static_cast(buf.ptr); + for (Omega_h::LO i = 0; i < buf.shape[0]; ++i) { + write_view_host[i] = ptr[i]; + } + auto write_view = Omega_h::Write(write_view_host); + return write_view; +} + + +// Helper to convert Omega_h::Write to numpy array (creates a reference) +template +py::array_t omega_h_write_to_numpy(Omega_h::Write write_view) { + return py::array_t( + {static_cast(write_view.size())}, // shape + {sizeof(T)}, // strides + write_view.data(), // data pointer + py::cast(write_view) // base object to manage lifetime + ); +} + + +// Helper function to convert numpy array to Rank2 View +template +Rank2View numpy_to_view_2d(py::array_t arr) { + py::buffer_info buf = arr.request(); + if (buf.ndim != 2) { + throw std::runtime_error("Number of dimensions must be 2"); + } + Rank2View view( + reinterpret_cast(buf.ptr), buf.shape[0], buf.shape[1]); + return view; +} + + +// Helper function to convert Rank2 View to array (creates a reference) +template +py::array_t view_2d_to_numpy(Rank2View view) { + return py::array_t( + {static_cast(view.extent(0)), static_cast(view.extent(1))}, // shape + {sizeof(T) * view.extent(1), sizeof(T)}, // strides (row-major) + view.data_handle(), // data pointer + py::cast(view) // base object to manage lifetime + ); +} + + +// Helper function to convert Rank3 View to array (creates a reference) +template +py::array_t view_3d_to_numpy(Rank3View view) { + return py::array_t( + {static_cast(view.extent(0)), + static_cast(view.extent(1)), + static_cast(view.extent(2))}, // shape + {sizeof(T) * view.extent(1) * view.extent(2), + sizeof(T) * view.extent(2), + sizeof(T)}, // strides (row-major) + view.data_handle(), // data pointer + py::cast(view) // base object to manage lifetime + ); +} + + +// Helper function to convert mdspan View to array (creates a reference) +template +py::array_t::type> mdspan_view_to_numpy( + View<2, T, HostMemorySpace> view) { + using NonConstT = typename std::remove_const::type; + py::array_t result( + {static_cast(view.extent(0)), + static_cast(view.extent(1))}); + auto out = result.template mutable_unchecked<2>(); + for (py::ssize_t i = 0; i < out.shape(0); ++i) { + for (py::ssize_t j = 0; j < out.shape(1); ++j) { + out(i, j) = view(i, j); + } + } + return result; +} + +template +py::array_t::type> mdspan_view_to_numpy( + View<3, T, HostMemorySpace> view) { + using NonConstT = typename std::remove_const::type; + py::array_t result( + {static_cast(view.extent(0)), + static_cast(view.extent(1)), + static_cast(view.extent(2))}); + auto out = result.template mutable_unchecked<3>(); + for (py::ssize_t i = 0; i < out.shape(0); ++i) { + for (py::ssize_t j = 0; j < out.shape(1); ++j) { + for (py::ssize_t k = 0; k < out.shape(2); ++k) { + out(i, j, k) = view(i, j, k); + } + } + } + return result; +} + + +} // namespace pcms \ No newline at end of file diff --git a/src/pcms/pythonapi/pythonapi.cpp b/src/pcms/pythonapi/pythonapi.cpp new file mode 100644 index 00000000..7fb274ed --- /dev/null +++ b/src/pcms/pythonapi/pythonapi.cpp @@ -0,0 +1,53 @@ + +#include + +namespace py = pybind11; + +namespace pcms { + +void bind_omega_h_field2(py::module& m); + +void bind_transfer_field2_module(py::module& m); + +void bind_omega_h_mesh_module(py::module& m); + +// void bind_interpolation_base_module(py::module& m); + +void bind_coordinate_system_module(py::module& m); + +void bind_coordinate_module(py::module& m); + +void bind_create_field_module(py::module& m); + +void bind_field_layout_module(py::module& m); + +void bind_field_module(py::module& m); + +void bind_omega_h_field_layout_module(py::module& m); + +void bind_uniform_grid_field_layout_module(py::module& m); + +void bind_uniform_grid_field_module(py::module& m); + +} +PYBIND11_MODULE(py_pcms, m) { + // Bind fundamental types first (coordinate systems, etc.) + pcms::bind_coordinate_system_module(m); + pcms::bind_coordinate_module(m); + + // Bind mesh and field infrastructure + pcms::bind_omega_h_mesh_module(m); + pcms::bind_field_layout_module(m); + pcms::bind_field_module(m); + pcms::bind_omega_h_field_layout_module(m); + pcms::bind_uniform_grid_field_layout_module(m); + pcms::bind_create_field_module(m); + + // Bind field types + pcms::bind_omega_h_field2(m); + pcms::bind_uniform_grid_field_module(m); + + // Bind field operations + pcms::bind_transfer_field2_module(m); + // pcms::bind_interpolation_base_module(m); +} \ No newline at end of file diff --git a/src/pcms/pythonapi/test_field_copy.py b/src/pcms/pythonapi/test_field_copy.py new file mode 100644 index 00000000..36d2ee23 --- /dev/null +++ b/src/pcms/pythonapi/test_field_copy.py @@ -0,0 +1,84 @@ +import py_pcms +import numpy as np + +def test_copy(world, dim, order, num_components): + """Test copying omega_h_field2 data""" + nx = 100 + ny = 100 if dim > 1 else 0 + nz = 100 if dim > 2 else 0 + print(f"\nStarting test: dim={dim}, order={order}, num_components={num_components}") + + # Build mesh + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, + nx, ny, nz, + False + ) + print(f" Built {dim}D mesh with {nx}x{ny}x{nz} elements") + print(f" Mesh type: {type(mesh)}") + print(f" Mesh object: {mesh}") + + # Create layout + print(" About to create layout...") + layout = py_pcms.create_lagrange_layout( + mesh, + order, + num_components, + py_pcms.CoordinateSystem.Cartesian + ) + print(f" Layout created successfully") + print(f"Testing dim={dim}, order={order}, num_components={num_components}...") + + # Get number of data points + ndata = layout.get_num_owned_dof_holder() * num_components + print(f" Number of data points: {ndata}") + + # Create sequential array of IDs + ids = np.arange(ndata, dtype=np.float64) + print(f" Created array of IDs from 0 to {ndata-1}") + + # Create original field and set data + original = layout.create_field() + print(" Created original field") + original.set_dof_holder_data(ids) + print(" Set data in original field") + + # Create copied field and copy data + copied = layout.create_field() + py_pcms.copy_field2_Real(original, copied) + print(" Copied data to new field") + + # Get copied data + copied_array = copied.get_dof_holder_data() + + # Verify the copy + assert len(copied_array) == ndata, f"Expected {ndata} elements, got {len(copied_array)}" + + # Check that all values match + differences = np.abs(ids - copied_array) + num_matches = np.sum(differences < 1e-12) + + assert num_matches == ndata, f"Only {num_matches}/{ndata} elements matched" + + print(f"✓ Test passed: dim={dim}, order={order}, num_components={num_components}") + +def main(): + """Run all test cases""" + print("Testing copy omega_h_field2 data...") + + # Initialize Omega_h library + lib = py_pcms.OmegaHLibrary() + world = lib.world() + print("Initialized Omega_h library and world") + + # Run test cases + test_copy(world, 2, 1, 1) + print("first passed") + test_copy(world, 2, 2, 1) + + print("\nAll tests passed!") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/pcms/pythonapi/test_file_io.py b/src/pcms/pythonapi/test_file_io.py new file mode 100644 index 00000000..a688e7a4 --- /dev/null +++ b/src/pcms/pythonapi/test_file_io.py @@ -0,0 +1,452 @@ +#!/usr/bin/env python3 +""" +Test Omega_h file I/O Python bindings +""" + +import py_pcms +import numpy as np +import os +import shutil + +def test_binary_io(lib, world): + """Test binary file format I/O""" + print("\n=== Testing Binary Format I/O ===") + + # Build a simple 3D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, # x, y, z dimensions + 4, 4, 4, # nx, ny, nz divisions + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test binary write/read + binary_file = os.path.join(test_dir, "test_mesh.osh") + print(f"Writing to binary file: {binary_file}") + py_pcms.write_mesh_binary(binary_file, mesh) + + print(f"Reading from binary file: {binary_file}") + mesh_read = py_pcms.read_mesh_binary(binary_file, lib) + + # Verify mesh properties + assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" + assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" + assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" + + print("✓ Binary I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_gmsh_io(lib, world): + """Test Gmsh file format I/O""" + print("\n=== Testing Gmsh Format I/O ===") + + # Build a simple 2D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 2.0, 1.0, 0.0, # x, y, z dimensions (z=0 for 2D) + 5, 3, 0, # nx, ny, nz divisions (nz=0 for 2D) + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test Gmsh write/read + gmsh_file = os.path.join(test_dir, "test_mesh.msh") + print(f"Writing to Gmsh file: {gmsh_file}") + py_pcms.write_mesh_gmsh(gmsh_file, mesh) + + print(f"Reading from Gmsh file: {gmsh_file}") + mesh_read = py_pcms.read_mesh_gmsh(gmsh_file, world) + + # Verify mesh properties + assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" + assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" + assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" + + print("✓ Gmsh I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_vtk_io(lib, world): + """Test VTK file format I/O""" + print("\n=== Testing VTK Format I/O ===") + + # Build a simple 3D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.5, 1.0, 0.5, # x, y, z dimensions + 3, 3, 2, # nx, ny, nz divisions + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + coords_orig = mesh.coords() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + print(f"Coordinates shape: {coords_orig.shape}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test VTU write (VTU is write-only in the API, typically for visualization) + vtu_file = os.path.join(test_dir, "test_mesh.vtu") + print(f"Writing to VTU file: {vtu_file}") + py_pcms.write_mesh_vtu(vtu_file, mesh, compress=False) + + # Verify file was created + assert os.path.exists(vtu_file), f"VTU file was not created: {vtu_file}" + file_size = os.path.getsize(vtu_file) + print(f"VTU file created successfully (size: {file_size} bytes)") + + # Test compressed VTU write + vtu_compressed = os.path.join(test_dir, "test_mesh_compressed.vtu") + print(f"Writing compressed VTU file: {vtu_compressed}") + py_pcms.write_mesh_vtu(vtu_compressed, mesh, compress=True) + + assert os.path.exists(vtu_compressed), f"Compressed VTU file was not created" + compressed_size = os.path.getsize(vtu_compressed) + print(f"Compressed VTU file created (size: {compressed_size} bytes)") + + print("✓ VTK I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_meshb_io(lib, world): + """Test MESHB file format I/O""" + print("\n=== Testing MESHB Format I/O ===") + + # Check if MESHB support is available + if not hasattr(py_pcms, 'write_mesh_meshb'): + print("⊘ MESHB support not available (OMEGA_H_USE_LIBMESHB not enabled)") + return + + # Build a simple 3D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, # x, y, z dimensions + 3, 3, 3, # nx, ny, nz divisions + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test MESHB write/read + meshb_file = os.path.join(test_dir, "test_mesh.mesh") + print(f"Writing to MESHB file: {meshb_file}") + py_pcms.write_mesh_meshb(mesh, meshb_file, version=2) + + print(f"Reading from MESHB file: {meshb_file}") + # MESHB read requires pre-created mesh object + mesh_read = py_pcms.OmegaHMesh(lib) + mesh_read.set_comm(world) + py_pcms.read_mesh_meshb(mesh_read, meshb_file) + + # Verify mesh properties + assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" + assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" + assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" + + print("✓ MESHB I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_exodus_io(lib, world): + """Test Exodus file format I/O""" + print("\n=== Testing Exodus Format I/O ===") + + # Check if Exodus support is available + if not hasattr(py_pcms, 'write_mesh_exodus'): + print("⊘ Exodus support not available (OMEGA_H_USE_SEACASEXODUS not enabled)") + return + + # Build a simple 3D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, # x, y, z dimensions + 3, 3, 3, # nx, ny, nz divisions + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test Exodus write + exodus_file = os.path.join(test_dir, "test_mesh.exo") + print(f"Writing to Exodus file: {exodus_file}") + py_pcms.write_mesh_exodus(exodus_file, mesh, verbose=False) + + # Test Exodus file handle API + if hasattr(py_pcms, 'exodus_open'): + print(f"Testing Exodus file handle API with: {exodus_file}") + + # Open the file + exo_handle = py_pcms.exodus_open(exodus_file, verbose=False) + print(f"Opened Exodus file with handle: {exo_handle}") + + # Get number of time steps + num_steps = py_pcms.exodus_get_num_time_steps(exo_handle) + print(f"Number of time steps: {num_steps}") + + # Read mesh using file handle + mesh_from_handle = py_pcms.OmegaHMesh(lib) + mesh_from_handle.set_comm(world) + py_pcms.read_mesh_exodus(exo_handle, mesh_from_handle, verbose=False) + + # Verify mesh properties + assert mesh_from_handle.dim() == dim_orig, f"Dimension mismatch: {mesh_from_handle.dim()} != {dim_orig}" + assert mesh_from_handle.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_from_handle.nverts()} != {nverts_orig}" + assert mesh_from_handle.nelems() == nelems_orig, f"Element count mismatch: {mesh_from_handle.nelems()} != {nelems_orig}" + + # Close the file + py_pcms.exodus_close(exo_handle) + + print("✓ Exodus I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_adios2_io(lib, world): + """Test ADIOS2 file format I/O""" + print("\n=== Testing ADIOS2 Format I/O ===") + + # Check if ADIOS2 support is available + if not hasattr(py_pcms, 'write_mesh_adios2'): + print("⊘ ADIOS2 support not available (OMEGA_H_USE_ADIOS2 not enabled)") + return + + # Build a simple 3D box mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, # x, y, z dimensions + 3, 3, 3, # nx, ny, nz divisions + False # symmetric + ) + + # Get initial mesh properties + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + dim_orig = mesh.dim() + print(f"Original mesh: dim={dim_orig}, nverts={nverts_orig}, nelems={nelems_orig}") + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Test ADIOS2 write/read + adios2_file = os.path.join(test_dir, "test_mesh.bp") + print(f"Writing to ADIOS2 file: {adios2_file}") + py_pcms.write_mesh_adios2(adios2_file, mesh, prefix="") + + print(f"Reading from ADIOS2 file: {adios2_file}") + mesh_read = py_pcms.read_mesh_adios2(adios2_file, lib, prefix="") + + # Verify mesh properties + assert mesh_read.dim() == dim_orig, f"Dimension mismatch: {mesh_read.dim()} != {dim_orig}" + assert mesh_read.nverts() == nverts_orig, f"Vertex count mismatch: {mesh_read.nverts()} != {nverts_orig}" + assert mesh_read.nelems() == nelems_orig, f"Element count mismatch: {mesh_read.nelems()} != {nelems_orig}" + + print("✓ ADIOS2 I/O test passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + + +def test_read_mesh_file_auto_detect(lib, world): + """Test automatic format detection with read_mesh_file""" + print("\n=== Testing Automatic Format Detection ===") + + # Build a test mesh + print("Creating test mesh...") + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, + 3, 3, 3, + False + ) + + nverts_orig = mesh.nverts() + nelems_orig = mesh.nelems() + + # Create temporary directory for test files + test_dir = "test_data" + os.makedirs(test_dir, exist_ok=True) + try: + # Write in different formats + binary_file = os.path.join(test_dir, "mesh.osh") + gmsh_file = os.path.join(test_dir, "mesh.msh") + + py_pcms.write_mesh_binary(binary_file, mesh) + py_pcms.write_mesh_gmsh(gmsh_file, mesh) + + # Test auto-detection for binary format + print(f"Auto-detecting binary file: {binary_file}") + mesh_binary = py_pcms.read_mesh_file(binary_file, world) + assert mesh_binary.nverts() == nverts_orig + assert mesh_binary.nelems() == nelems_orig + print("✓ Binary auto-detection passed") + + # Test auto-detection for Gmsh format + print(f"Auto-detecting Gmsh file: {gmsh_file}") + mesh_gmsh = py_pcms.read_mesh_file(gmsh_file, world) + assert mesh_gmsh.nverts() == nverts_orig + assert mesh_gmsh.nelems() == nelems_orig + print("✓ Gmsh auto-detection passed") + + finally: + # Clean up temporary files + if os.path.exists(test_dir): + for item in os.listdir(test_dir): + item_path = os.path.join(test_dir, item) + if os.path.isfile(item_path): + os.remove(item_path) + elif os.path.isdir(item_path): + shutil.rmtree(item_path) + os.rmdir(test_dir) + print(f"Cleaned up test directory: {test_dir}") + +if __name__ == "__main__": + print("=" * 60) + print("Omega_h File I/O Python Binding Tests") + print("=" * 60) + + # Create library and world communicator once for all tests + lib = py_pcms.OmegaHLibrary() + world = lib.world() + + try: + test_binary_io(lib, world) + test_gmsh_io(lib, world) + test_vtk_io(lib, world) + test_meshb_io(lib, world) + test_exodus_io(lib, world) + test_adios2_io(lib, world) + test_read_mesh_file_auto_detect(lib, world) + + print("\n" + "=" * 60) + print("✓ All tests passed!") + print("=" * 60) + + except Exception as e: + print("\n" + "=" * 60) + print(f"✗ Test failed with error: {e}") + print("=" * 60) + import traceback + traceback.print_exc() + exit(1) + finally: + # Explicitly delete objects to avoid an MPI finalizing issue + del world + del lib diff --git a/src/pcms/pythonapi/test_omega_h_field.py b/src/pcms/pythonapi/test_omega_h_field.py new file mode 100644 index 00000000..c0c2053a --- /dev/null +++ b/src/pcms/pythonapi/test_omega_h_field.py @@ -0,0 +1,276 @@ +#!/usr/bin/env python3 +""" +Test OmegaHFieldLayout Python bindings +""" + +import py_pcms +import numpy as np + +def test_layout_methods(world, dim, order, num_components): + """Test OmegaHFieldLayout methods""" + nx = 10 + ny = 10 if dim > 1 else 0 + nz = 10 if dim > 2 else 0 + + print(f"\nTesting layout: dim={dim}, order={order}, num_components={num_components}") + + # Build mesh + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, + nx, ny, nz, + False + ) + + # Create layout + layout = py_pcms.create_lagrange_layout( + mesh, + order, + num_components, + py_pcms.CoordinateSystem.Cartesian + ) + + # Test layout methods + num_comp = layout.get_num_components() + assert num_comp == num_components, f"Expected {num_components} components, got {num_comp}" + + num_owned = layout.get_num_owned_dof_holder() + print(f" Number of owned DOF holders: {num_owned}") + assert num_owned > 0, "Expected at least one owned DOF holder" + + num_global = layout.get_num_global_dof_holder() + print(f" Number of global DOF holders: {num_global}") + assert num_global >= num_owned, "Global DOF holders should be >= owned" + + # Get ownership mask + owned = layout.get_owned() + print(f" Ownership mask length: {len(owned)}") + assert len(owned) > 0, "Ownership mask should not be empty" + + # Get global IDs + gids = layout.get_gids() + print(f" Global IDs length: {len(gids)}") + assert len(gids) > 0, "Global IDs should not be empty" + + # Get coordinates of DOF holders (2D parametric coordinates) + coords = layout.get_dof_holder_coordinates() + print(f" Coordinates shape: {coords.shape}") + assert coords.shape[0] > 0, "Should have coordinates" + assert coords.shape[1] == 2, "Coordinates should be 2D (parametric)" + + # Check if distributed + is_distributed = layout.is_distributed() + print(f" Is distributed: {is_distributed}") + + # Get entity offsets + ent_offsets = layout.get_ent_offsets() + print(f" Entity offsets: {ent_offsets}") + assert len(ent_offsets) == 5, f"Expected 5 entity offsets" + + # Get nodes per dimension + nodes = layout.get_nodes_per_dim() + print(f" Nodes per dim: {nodes}") + assert len(nodes) == 4, "Should have 4 entries (verts, edges, faces, regions)" + + # Get number of entities + num_ents = layout.get_num_ents() + print(f" Number of entities: {num_ents}") + assert num_ents > 0, "Should have at least one entity" + + # Get sizes + owned_size = layout.owned_size() + global_size = layout.global_size() + print(f" Owned size: {owned_size}") + print(f" Global size: {global_size}") + assert owned_size == num_owned * num_components + assert global_size == num_global * num_components + + # Create a field from the layout + field = layout.create_field() + print(f" Created field: {type(field)}") + + # Test setting and getting data + ndata = num_owned * num_components + test_data = np.arange(ndata, dtype=np.float64) + field.set_dof_holder_data(test_data) + + retrieved_data = field.get_dof_holder_data() + assert len(retrieved_data) == ndata, f"Expected {ndata} elements, got {len(retrieved_data)}" + + # Verify data matches + # debug print + for i in range(min(10, ndata)): + print(f" Data[{i}]: set={test_data[i]}, got={retrieved_data[i]}") + differences = np.abs(test_data - retrieved_data) + assert np.all(differences < 1e-12), "Data mismatch after set/get" + + print(f"✓ Test passed: dim={dim}, order={order}, num_components={num_components}") + +def test_field_evaluation(world, dim, order, num_components): + """Test OmegaHField evaluation at points""" + nx = 10 + ny = 10 if dim > 1 else 0 + nz = 10 if dim > 2 else 0 + + print(f"\nTesting field evaluation: dim={dim}, order={order}, num_components={num_components}") + + # Build mesh + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 1.0, + nx, ny, nz, + False + ) + + # Create layout + layout = py_pcms.create_lagrange_layout( + mesh, + order, + num_components, + py_pcms.CoordinateSystem.Cartesian + ) + + # Create field + field = layout.create_field() + + # Set up test data - use a simple function: f(x,y,z) = sin(x*y) for component 0, etc. + print(f" Setting up field data...") + # Get mesh vertex coordinates to set field values + mesh_coords = mesh.coords() + num_verts = mesh.nverts() + print(f" Mesh has {num_verts} vertices") + + # For order 1, we only need vertex values + # For order 2, we also need edge midpoint values + num_owned = layout.get_num_owned_dof_holder() + ndata = num_owned * num_components + print(f" Setting up field data for {ndata} DOFs") + + # Define test function + def test_func(x, y, z, component): + return np.sin(5.0 * x * y) + float(component) + + # Create test data array + test_data = np.zeros(ndata, dtype=np.float64) + + # Set vertex values + for i in range(min(num_verts, num_owned)): + x = mesh_coords[i * dim + 0] if dim >= 1 else 0.0 + y = mesh_coords[i * dim + 1] if dim >= 2 else 0.0 + z = mesh_coords[i * dim + 2] if dim >= 3 else 0.0 + for c in range(num_components): + idx = i * num_components + c + test_data[idx] = test_func(x, y, z, c) + + print(f" Set vertex values for field data") + + # For order 2, set edge midpoint values (simplified - would need proper edge coordinates) + if order == 2 and num_owned > num_verts: + for i in range(num_verts, num_owned): + for c in range(num_components): + idx = i * num_components + c + # Use simplified values for edges + test_data[idx] = float(c) + 0.5 + + field.set_dof_holder_data(test_data) + print(f" Set field data with test function") + + # Define evaluation points (physical coordinates) + if dim == 2: + eval_coords = np.array([ + [0.5, 0.5], + [0.25, 0.25], + [0.75, 0.75], + [0.1, 0.9], + [0.9, 0.1] + ], dtype=np.float64) + else: # dim == 3 + eval_coords = np.array([ + [0.5, 0.5, 0.5], + [0.25, 0.25, 0.25], + [0.75, 0.75, 0.75], + [0.1, 0.9, 0.1], + [0.9, 0.1, 0.9] + ], dtype=np.float64) + + num_eval_points = eval_coords.shape[0] + print(f" Evaluating at {num_eval_points} points") + + # Get localization hint for the coordinates + coord_system = py_pcms.CoordinateSystem.Cartesian + location_hint = field.get_localization_hint(eval_coords, coord_system) + print(f" Got localization hint") + + # Create output buffer for evaluation results + eval_values = np.zeros(num_eval_points * num_components, dtype=np.float64) + + # Evaluate the field + field.evaluate(location_hint, eval_values, coord_system) + print(f" Field evaluated successfully") + + # Print and verify results + for i in range(num_eval_points): + coords_str = ", ".join([f"{eval_coords[i, j]:.3f}" for j in range(dim)]) + values_str = ", ".join([f"{eval_values[i * num_components + c]:.4f}" + for c in range(num_components)]) + print(f" Point {i} ({coords_str}): values = [{values_str}]") + + # Verify we got valid values (not NaN or infinity) + for c in range(num_components): + val = eval_values[i * num_components + c] + assert not np.isnan(val), f"Got NaN value at point {i}, component {c}" + assert not np.isinf(val), f"Got inf value at point {i}, component {c}" + + # Test gradient evaluation if available + if field.can_evaluate_gradient(): + print(f" Testing gradient evaluation...") + gradient_values = np.zeros(num_eval_points * num_components * dim, dtype=np.float64) + + try: + field.evaluate_gradient(gradient_values, coord_system) + print(f" Gradient evaluated successfully") + + # Print gradient results + for i in range(min(3, num_eval_points)): + for c in range(num_components): + grad_components = [] + for d in range(dim): + idx = i * num_components * dim + c * dim + d + grad_components.append(f"{gradient_values[idx]:.4f}") + grad_str = ", ".join(grad_components) + print(f" Point {i}, component {c}: gradient = [{grad_str}]") + except Exception as e: + print(f" Gradient evaluation failed: {e}") + else: + print(f" Gradient evaluation not supported") + + print(f"✓ Field evaluation test passed: dim={dim}, order={order}, num_components={num_components}") + +def main(): + """Run all test cases""" + print("Testing OmegaHFieldLayout Python bindings...") + + # Initialize Omega_h library + lib = py_pcms.OmegaHLibrary() + world = lib.world() + print("Initialized Omega_h library and world") + + # Test different configurations + test_layout_methods(world, 2, 1, 1) + test_layout_methods(world, 2, 2, 1) + # test_layout_methods(world, 2, 1, 3) + # test_layout_methods(world, 3, 1, 1) + + # Test field evaluation + print("\n" + "="*60) + print("Testing field evaluation...") + print("="*60) + test_field_evaluation(world, 2, 1, 1) + test_field_evaluation(world, 2, 2, 1) + + print("\n✓ All tests passed!") + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/pcms/pythonapi/test_uniform_grid_field.py b/src/pcms/pythonapi/test_uniform_grid_field.py new file mode 100644 index 00000000..06109c78 --- /dev/null +++ b/src/pcms/pythonapi/test_uniform_grid_field.py @@ -0,0 +1,476 @@ +""" +Example test file for uniform grid field Python bindings. + +This demonstrates how to use the UniformGridField and UniformGridFieldLayout +classes in Python, similar to how OmegaHField is used. +""" +import numpy as np +import py_pcms + + +def test_uniform_grid_field_creation(): + """Test creating a 2D uniform grid field.""" + # Create a 2D uniform grid + grid = py_pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [4, 4] # 4x4 cells + + print(f"Grid cells: {grid.get_num_cells()}") # Should be 16 + + # Create field layout with 1 component (scalar field) + layout = py_pcms.UniformGridFieldLayout2D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + + print(f"Num components: {layout.get_num_components()}") + print(f"Num owned DOF holders: {layout.get_num_owned_dof_holder()}") + print(f"Num vertices: {layout.get_num_vertices()}") # 5x5 = 25 vertices + + # Create field + field = layout.create_field() + print(f"Field created: {field is not None}") + + return grid, layout, field + + +def test_uniform_grid_field_data_operations(): + """Test setting and getting field data.""" + # Create a simple grid + grid = py_pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [2, 2] # 2x2 cells, 3x3 vertices + + layout = py_pcms.UniformGridFieldLayout2D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + ug_field = layout.create_field() + + # Get the number of DOF holders (vertices) + num_vertices = layout.get_num_vertices() + print(f"Number of vertices: {num_vertices}") # Should be 9 + + # Create data for vertices (initialize with values) + data = np.arange(num_vertices, dtype=np.float64) + print(f"Setting data: {data}") + + # Set field data + ug_field.set_dof_holder_data(data) + + # Get field data back + retrieved_data = ug_field.get_dof_holder_data() + print(f"Retrieved data: {retrieved_data}") + + # Verify they match + assert np.allclose(data, retrieved_data) + print("Data successfully set and retrieved!") + + +def test_uniform_grid_field_mdspan_2d(): + """Test 2D mdspan output as numpy array.""" + grid = py_pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [2, 3] # 3x4 vertices + + layout = py_pcms.UniformGridFieldLayout2D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + ug_field = layout.create_field() + + num_vertices = layout.get_num_vertices() + data = np.arange(num_vertices, dtype=np.float64) + ug_field.set_dof_holder_data(data) + + mdspan = ug_field.to_mdspan() + expected = data.reshape((grid.divisions[0] + 1, grid.divisions[1] + 1)) + assert type(mdspan) == np.ndarray + + assert mdspan.shape == expected.shape + assert np.allclose(mdspan, expected) + print("2D mdspan output verified!") + + +def test_uniform_grid_field_evaluation(): + """Test evaluating field at specific coordinates.""" + # Create grid + grid = py_pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [2, 2] + + layout = py_pcms.UniformGridFieldLayout2D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + ug_field = layout.create_field() + + # Set simple linear field data: f(x,y) = x + y + num_vertices = layout.get_num_vertices() + coords = layout.get_dof_holder_coordinates() + data = np.zeros(num_vertices) + for i in range(num_vertices): + data[i] = coords[i, 0] + coords[i, 1] + + ug_field.set_dof_holder_data(data) + + # Evaluate at some points + eval_coords = np.array([[2.5, 2.5], [7.5, 7.5]]) + + # Get localization hint + hint = ug_field.get_localization_hint( + eval_coords, py_pcms.CoordinateSystem.Cartesian + ) + + # Evaluate field + results = np.zeros(len(eval_coords)) + ug_field.evaluate( + hint, results, py_pcms.CoordinateSystem.Cartesian + ) + + print(f"Evaluation results: {results}") + print(f"Expected (approximately): [5.0, 15.0]") + + +def test_uniform_grid_closest_cell(): + """Test finding closest cell to a point.""" + grid = py_pcms.UniformGrid2D() + grid.bot_left = [0.0, 0.0] + grid.edge_length = [10.0, 10.0] + grid.divisions = [4, 4] + + # Test point in middle of grid + point = np.array([5.0, 5.0]) + cell_id = grid.closest_cell_id(point) + print(f"Point {point} is in cell {cell_id}") + + # Test point outside grid (should clamp to closest cell) + point_outside = np.array([-1.0, -1.0]) + cell_id_outside = grid.closest_cell_id(point_outside) + print(f"Point {point_outside} (outside) maps to cell {cell_id_outside}") + + +def test_3d_uniform_grid(): + """Test 3D uniform grid field.""" + grid = py_pcms.UniformGrid3D() + grid.bot_left = [0.0, 0.0, 0.0] + grid.edge_length = [10.0, 10.0, 10.0] + grid.divisions = [2, 2, 2] # 2x2x2 cells + + print(f"3D Grid cells: {grid.get_num_cells()}") # Should be 8 + + layout = py_pcms.UniformGridFieldLayout3D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + + print(f"3D Grid vertices: {layout.get_num_vertices()}") # 3x3x3 = 27 + + ug_field = layout.create_field() + + # Set and get data + num_vertices = layout.get_num_vertices() + data = np.ones(num_vertices) * 42.0 + ug_field.set_dof_holder_data(data) + + retrieved = ug_field.get_dof_holder_data() + assert np.allclose(data, retrieved) + print("3D field data successfully set and retrieved!") + + +def test_uniform_grid_field_mdspan_3d(): + """Test 3D mdspan output as numpy array.""" + grid = py_pcms.UniformGrid3D() + grid.bot_left = [0.0, 0.0, 0.0] + grid.edge_length = [10.0, 10.0, 10.0] + grid.divisions = [2, 1, 3] # 3x2x4 vertices + + layout = py_pcms.UniformGridFieldLayout3D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + ug_field = layout.create_field() + + num_vertices = layout.get_num_vertices() + data = np.arange(num_vertices, dtype=np.float64) + ug_field.set_dof_holder_data(data) + + mdspan = ug_field.to_mdspan() + expected = data.reshape( + (grid.divisions[0] + 1, grid.divisions[1] + 1, grid.divisions[2] + 1) + ) + + assert mdspan.shape == expected.shape + assert np.allclose(mdspan, expected) + print("3D mdspan output verified!") + + +def test_uniform_grid_workflow(world): + """ + Test complete UniformGrid workflow with Omega_h mesh and field interpolation. + + This test: + 1. Creates an Omega_h mesh + 2. Creates a uniform grid from the mesh + 3. Creates a binary mask field + 4. Creates and initializes an Omega_h field with f(x,y) = x + 2*y + 5. Transfers the field to uniform grid via interpolation + 6. Verifies the transferred values + """ + # Create a simple 2D box mesh: 1.0 x 1.0 domain with 4x4 elements + mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 0.0, + 4, 4, 0, + False + ) + + # Create uniform grid from mesh with 4x4 divisions + grid = py_pcms.create_uniform_grid_from_mesh(mesh, [4, 4]) + print(f"Created uniform grid with {grid.get_num_cells()} cells") + + # Create binary mask field (returns tuple of (layout, field)) + mask_layout, mask_field = py_pcms.create_uniform_grid_binary_field(mesh, [4, 4]) + mask_data = mask_field.get_dof_holder_data() + print(f"Created mask field with {len(mask_data)} vertices") + + # Create Omega_h field layout with linear elements + omega_h_layout = py_pcms.create_lagrange_layout( + mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian + ) + omega_h_field = omega_h_layout.create_field() + + # Initialize omega_h field with f(x,y) = x + 2*y + coords = omega_h_layout.get_dof_holder_coordinates() + num_nodes = omega_h_layout.get_num_owned_dof_holder() + omega_h_data = np.zeros(num_nodes) + + for i in range(num_nodes): + x = coords[i, 0] + y = coords[i, 1] + omega_h_data[i] = x + 2.0 * y + + omega_h_field.set_dof_holder_data(omega_h_data) + print(f"Initialized Omega_h field with {num_nodes} nodes") + + # Create uniform grid field layout + ug_layout = py_pcms.UniformGridFieldLayout2D( + grid, 1, py_pcms.CoordinateSystem.Cartesian + ) + ug_field = ug_layout.create_field() + + # Transfer from omega_h field to uniform grid field using interpolation + py_pcms.interpolate_field(omega_h_field, ug_field) + print("Field interpolation completed") + + # Get uniform grid field data and coordinates + ug_field_data = ug_field.get_dof_holder_data() + ug_coords = ug_layout.get_dof_holder_coordinates() + + # Verify ug_field values + print("\nVerifying uniform grid field values:") + num_errors = 0 + for j in range(grid.divisions[1] + 1): + for i in range(grid.divisions[0] + 1): + vertex_id = j * (grid.divisions[0] + 1) + i + x = ug_coords[vertex_id, 0] + y = ug_coords[vertex_id, 1] + expected = x + 2.0 * y + actual = ug_field_data[vertex_id] + + error = abs(expected - actual) + if error > 1e-10: + num_errors += 1 + print(f" Vertex ({i}, {j}) at ({x:.2f}, {y:.2f}): " + f"expected {expected:.4f}, got {actual:.4f}, error {error:.2e}") + + if num_errors == 0: + print("All uniform grid field values verified successfully!") + else: + print(f"Found {num_errors} verification errors") + + # Verify mask field values (all should be 1 for vertices inside the mesh) + print("\nVerifying mask field values:") + mask_errors = 0 + nx, ny = grid.divisions[0] + 1, grid.divisions[1] + 1 + print(f"\nBinary mask field ({nx}x{ny} vertices) visualization:") + print("(Each position shows the mask value at that vertex)") + for j in range(ny - 1, -1, -1): # Print from top to bottom + row = [] + for i in range(nx): + vertex_id = j * nx + i + mask_value = mask_data[vertex_id] + row.append(str(int(mask_value))) + if mask_value != 1.0: + mask_errors += 1 + print(f" Row {j}: [{' '.join(row)}]") + + if mask_errors == 0: + print("All mask field values verified successfully!") + else: + print(f"Found {mask_errors} mask errors") + + # Serialize the uniform grid field + print("\nSerializing uniform grid field:") + num_vertices = ug_layout.get_num_vertices() + buffer = np.zeros(num_vertices) + permutation = np.arange(num_vertices, dtype=np.int32) # Identity permutation + + bytes_written = ug_field.serialize(buffer, permutation) + print(f"Serialized {bytes_written} values to buffer") + + print(f"\nSerialized uniform grid field data ({grid.divisions[0]+1}x{grid.divisions[1]+1} vertices) visualization:") + print("(Each vertex shows its field value)") + for j in range(grid.divisions[1], -1, -1): # Print from top to bottom + row_values = [] + for i in range(grid.divisions[0] + 1): + vertex_id = j * (grid.divisions[0] + 1) + i + value = buffer[vertex_id] + row_values.append(f"{value:6.3f}") + print(f" Row {j}: [{' '.join(row_values)}]") + + # Also print with coordinates for reference + print(f"\nDetailed vertex information:") + for j in range(grid.divisions[1] + 1): + for i in range(grid.divisions[0] + 1): + vertex_id = j * (grid.divisions[0] + 1) + i + x = ug_coords[vertex_id, 0] + y = ug_coords[vertex_id, 1] + value = buffer[vertex_id] + print(f" V[{i},{j}] (id={vertex_id:2d}) at ({x:.3f}, {y:.3f}): value = {value:.6f}") + + # Test deserialization + print("\nTesting deserialization:") + ug_field_copy = ug_layout.create_field() + ug_field_copy.deserialize(buffer, permutation) + + ug_field_copy_data = ug_field_copy.get_dof_holder_data() + if np.allclose(ug_field_data, ug_field_copy_data): + print("✓ Deserialization successful - data matches!") + else: + print("✗ Deserialization failed - data mismatch!") + print(f" Max difference: {np.max(np.abs(ug_field_data - ug_field_copy_data))}") + + assert num_errors == 0, f"Field verification failed with {num_errors} errors" + assert mask_errors == 0, f"Mask verification failed with {mask_errors} errors" + + +def test_omega_h_to_omega_h_transfer_workflow(world): + """ + Test field transfer from one Omega_h mesh to another Omega_h mesh. + + This test: + 1. Creates a source Omega_h mesh + 2. Creates a target Omega_h mesh + 3. Creates Lagrange layouts and fields on both + 4. Initializes the source with f(x,y) = x + 2*y + 5. Transfers the field to the target mesh + 6. Verifies the transferred values at target DOF holders + """ + # Source mesh (coarser) + src_mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 0.0, + 4, 4, 0, + False + ) + + # Target mesh (finer) + tgt_mesh = py_pcms.build_box( + world, + py_pcms.Family.SIMPLEX, + 1.0, 1.0, 0.0, + 8, 8, 0, + False + ) + + # Create Lagrange layouts and fields + src_layout = py_pcms.create_lagrange_layout( + src_mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian + ) + tgt_layout = py_pcms.create_lagrange_layout( + tgt_mesh, 1, 1, py_pcms.CoordinateSystem.Cartesian + ) + src_field = src_layout.create_field() + tgt_field = tgt_layout.create_field() + + # Initialize source field with f(x,y) = x + 2*y + src_coords = src_layout.get_dof_holder_coordinates() + src_num_nodes = src_layout.get_num_owned_dof_holder() + src_data = np.zeros(src_num_nodes) + for i in range(src_num_nodes): + x = src_coords[i, 0] + y = src_coords[i, 1] + src_data[i] = x + 2.0 * y + src_field.set_dof_holder_data(src_data) + + # Transfer field to target mesh + py_pcms.interpolate_field(src_field, tgt_field) + + # Verify target field values at target DOF holders + tgt_coords = tgt_layout.get_dof_holder_coordinates() + tgt_data = tgt_field.get_dof_holder_data() + errors = 0 + for i in range(tgt_layout.get_num_owned_dof_holder()): + x = tgt_coords[i, 0] + y = tgt_coords[i, 1] + expected = x + 2.0 * y + actual = tgt_data[i] + if abs(expected - actual) > 1e-10: + errors += 1 + assert errors == 0, f"Omega_h transfer verification failed with {errors} errors" + + +if __name__ == "__main__": + lib = py_pcms.OmegaHLibrary() + world = lib.world() + print("=" * 60) + print("Testing UniformGrid Field Creation") + print("=" * 60) + test_uniform_grid_field_creation() + + print("\n" + "=" * 60) + print("Testing UniformGrid Field Data Operations") + print("=" * 60) + test_uniform_grid_field_data_operations() + + print("\n" + "=" * 60) + print("Testing UniformGrid Field mdspan (2D)") + print("=" * 60) + test_uniform_grid_field_mdspan_2d() + + print("\n" + "=" * 60) + print("Testing UniformGrid Field Evaluation") + print("=" * 60) + test_uniform_grid_field_evaluation() + + print("\n" + "=" * 60) + print("Testing Closest Cell ID") + print("=" * 60) + test_uniform_grid_closest_cell() + + print("\n" + "=" * 60) + print("Testing 3D UniformGrid") + print("=" * 60) + test_3d_uniform_grid() + + print("\n" + "=" * 60) + print("Testing UniformGrid Field mdspan (3D)") + print("=" * 60) + test_uniform_grid_field_mdspan_3d() + + print("\n" + "=" * 60) + print("Testing UniformGrid Workflow") + print("=" * 60) + test_uniform_grid_workflow(world) + + print("\n" + "=" * 60) + print("Testing Omega_h to Omega_h Transfer Workflow") + print("=" * 60) + test_omega_h_to_omega_h_transfer_workflow(world) + + print("\n" + "=" * 60) + print("All tests passed!") + print("=" * 60) + + del world diff --git a/src/pcms/pythonapi/uniform_grid_workflow_guide.md b/src/pcms/pythonapi/uniform_grid_workflow_guide.md new file mode 100644 index 00000000..8614e4da --- /dev/null +++ b/src/pcms/pythonapi/uniform_grid_workflow_guide.md @@ -0,0 +1,264 @@ +# Unstructured mesh field transfer to structured mesh field workflow guide + +This guide demonstrates how to transfer unstructured mesh fields to uniform grid fields in PCMS Python API, including mesh creation, field interpolation, and data serialization. + +## Overview + +The workflow allows you to: +1. Create a structured uniform grid as the bounding box as an unstructured Omega_h mesh +2. Create binary mask fields indicating which cells are inside the mesh +3. Interpolate field data from Omega_h mesh fields to uniform grid fields +4. Serialize and deserialize field data for I/O operations + +## Complete Workflow Example + +### Step 1: Initialize Omega_h Library + +Initialize the Omega_h library which manages parallel communication. + +```python +import py_pcms +import numpy as np + +# Create library and world communicator +lib = py_pcms.OmegaHLibrary() +world = lib.world() +``` + + +--- + +### Step 2: Create an Omega_h Mesh + +You can either create a mesh or read an existing mesh from a file. + +#### Option A: Create Mesh Programmatically + +```python +# Create a 2D box mesh: 1.0 x 1.0 domain with 4x4 elements +mesh = py_pcms.build_box( + world, # Communicator + py_pcms.Family.SIMPLEX, # Element type + 1.0, 1.0, 0.0, # Domain size: x, y, z + 4, 4, 0, # Number of elements: nx, ny, nz + False # Symmetric flag +) +``` + +#### Option B: Read Mesh from File + +```python +# Auto-detect file format and read +mesh = py_pcms.read_mesh_file("my_mesh.osh", world) + +# Or use format-specific readers for explicit control: + +# Binary format (.osh) +mesh = py_pcms.read_mesh_binary("mesh.osh", lib) + +# Gmsh format (.msh) +mesh = py_pcms.read_mesh_gmsh("mesh.msh", world) + +# VTK format (.pvtu for parallel files) +mesh = py_pcms.read_mesh_parallel_vtk("mesh.pvtu", world) + +# Exodus format (.exo) - if SEACAS support is enabled +# Using file handle API +exo_handle = py_pcms.exodus_open("mesh.exo", verbose=False) +num_steps = py_pcms.exodus_get_num_time_steps(exo_handle) +mesh = py_pcms.OmegaHMesh(lib) +mesh.set_comm(world) +py_pcms.read_mesh_exodus(exo_handle, mesh, verbose=False) +py_pcms.exodus_close(exo_handle) + +# ADIOS2 format (.bp) - if ADIOS2 support is enabled +mesh = py_pcms.read_mesh_adios2("mesh.bp", lib, prefix="") + +# MESHB format (.mesh) - if LIBMESHB support is enabled +mesh = py_pcms.OmegaHMesh(lib) +mesh.set_comm(world) +py_pcms.read_mesh_meshb(mesh, "mesh.mesh") +# Read solution/resolution data if available +py_pcms.read_meshb_sol(mesh, "mesh.sol", sol_name="resolution") +``` + +**PS**: File formats such as Exodus or libMeshB require that PCMS be built with the respective libraries enabled. + +Checking Mesh Properties: + +```python +print(f"Mesh dimension: {mesh.dim()}") +print(f"Number of vertices: {mesh.nverts()}") +print(f"Number of elements: {mesh.nelems()}") +print(f"Mesh family: {mesh.family()}") +``` + +--- + +### Step 3: Create Uniform Grid from Mesh + +Generate a structured uniform grid that covers the bounding box of the mesh. + +```python +# Create uniform grid with 4x4 cell divisions +grid = py_pcms.create_uniform_grid_from_mesh(mesh, [4, 4]) +print(f"Created uniform grid with {grid.get_num_cells()} cells") +``` + +--- + +### Step 4: Create Binary Mask Field + +Create a mask where each uniform grid vertex is marked as 1 (inside mesh) or 0 (outside mesh). + +```python +# Create binary mask indicating which vertices are inside the mesh +# Returns tuple of (layout, field) - layout must be kept alive while using field +mask_layout, mask_field = py_pcms.create_uniform_grid_binary_field(mesh, [4, 4]) +print(f"Created mask field with {mask_layout.get_num_vertices()} vertices") +``` + +**Accessing Values**: +```python +# Access mask value for vertex (i, j) +vertex_id = j * (grid.divisions[0] + 1) + i +mask_data = mask_field.get_dof_holder_data() +mask_value = mask_data[vertex_id] # Returns 0.0 or 1.0 +``` + +--- + +### Step 5: Create and Initialize Omega_h Field + +Create a field on the unstructured mesh and initialize it with data. + +```python +# Create field layout with linear (order=1) elements and 1 component (scalar) +omega_h_layout = py_pcms.create_lagrange_layout( + mesh, + 1, + 1, + py_pcms.CoordinateSystem.Cartesian +) +omega_h_field = omega_h_layout.create_field() + +# Initialize field with f(x,y) = x + 2*y +coords = omega_h_layout.get_dof_holder_coordinates() +num_nodes = omega_h_layout.get_num_owned_dof_holder() +omega_h_data = np.zeros(num_nodes) + +for i in range(num_nodes): + x = coords[i, 0] + y = coords[i, 1] + omega_h_data[i] = x + 2.0 * y + +omega_h_field.set_dof_holder_data(omega_h_data) +``` + +--- + +### Step 6: Create Uniform Grid Field Layout + +Set up the data structure for storing field values on the uniform grid vertices. + +```python +# Create uniform grid field layout +ug_layout = py_pcms.UniformGridFieldLayout2D( + grid, + 1, # Number of components + py_pcms.CoordinateSystem.Cartesian +) +ug_field = ug_layout.create_field() +``` + +--- + +### Step 7: Interpolate Field Data + +Interpolate field values from the unstructured mesh to the structured uniform grid vertices. + +```python +# Transfer field from Omega_h mesh to uniform grid +py_pcms.interpolate_field(omega_h_field, ug_field) +``` + +--- + +### Step 8: Access Field Data + +Retrieve interpolated values and verify correctness. + +```python +# Get field data and coordinates +ug_field_data = ug_field.get_dof_holder_data() +ug_coords = ug_layout.get_dof_holder_coordinates() + +# Access data at vertex (i, j) +vertex_id = j * (grid.divisions[0] + 1) + i +value = ug_field_data[vertex_id] +x, y = ug_coords[vertex_id, 0], ug_coords[vertex_id, 1] +``` + +--- + +### Step 9: Export Field Data with `to_mdspan` + +Convert field data to a structured $x \times y$ (or $x \times y \times z$) array +for downstream analyses. + +```python +# Convert field data to a structured array +grid_values = ug_field.to_mdspan() # 2D: (nx+1, ny+1), 3D: (nx+1, ny+1, nz+1) +mask_values = mask_field.to_mdspan() # 2D: (nx+1, ny+1), 3D: (nx+1, ny+1, nz+1) + +# Save to file (example) +np.save('field_data.npy', grid_values) +``` + +--- + +## API Reference Summary + +### Grid Creation +- `create_uniform_grid_from_mesh(mesh, divisions)` - Create grid from mesh +- `create_uniform_grid_binary_field(mesh, divisions)` - Create inside/outside mask + +### Field Layout +- `create_lagrange_layout(mesh, order, num_components, coord_system)` - Omega_h field layout +- `UniformGridFieldLayout2D(grid, num_components, coord_system)` - 2D grid layout +- `UniformGridFieldLayout3D(grid, num_components, coord_system)` - 3D grid layout + +### Field Operations +- `layout.create_field()` - Create field from layout +- `field.set_dof_holder_data(data)` - Set field values +- `field.get_dof_holder_data()` - Get field values +- `field.to_mdspan()` - Get field values as a structured array +- `interpolate_field(source_field, target_field)` - Interpolate between fields + +### Exodus I/O (if OMEGA_H_USE_SEACASEXODUS enabled) +- `exodus_open(filepath, verbose=False)` - Open Exodus file, returns file handle +- `exodus_close(exodus_file)` - Close Exodus file handle +- `exodus_get_num_time_steps(exodus_file)` - Get number of time steps in file +- `read_mesh_exodus(exodus_file, mesh, verbose=False, classify_with=...)` - Read mesh from file handle +- `read_exodus_nodal_fields(exodus_file, mesh, time_step, prefix="", postfix="", verbose=False)` - Read nodal fields +- `read_exodus_element_fields(exodus_file, mesh, time_step, prefix="", postfix="", verbose=False)` - Read element fields +- `read_mesh_exodus_sliced(filepath, comm, verbose=False, classify_with=..., time_step=-1)` - Read mesh in parallel +- `write_mesh_exodus(filepath, mesh, verbose=False, classify_with=...)` - Write mesh to Exodus file +- `ExodusClassifyWith.NODE_SETS` - Classify using node sets +- `ExodusClassifyWith.SIDE_SETS` - Classify using side sets + +### MESHB I/O (if OMEGA_H_USE_LIBMESHB enabled) +- `read_mesh_meshb(mesh, filepath)` - Read mesh from MESHB file +- `write_mesh_meshb(mesh, filepath, version)` - Write mesh to MESHB file +- `read_meshb_sol(mesh, filepath, sol_name)` - Read solution/resolution from .sol file +- `write_meshb_sol(mesh, filepath, sol_name, version)` - Write solution/resolution to .sol file + +### Grid Properties +- `grid.get_num_cells()` - Total number of cells +- `grid.divisions` - Cell divisions in each dimension +- `layout.get_num_vertices()` - Total number of vertices +- `layout.get_dof_holder_coordinates()` - Vertex coordinates + +--- + +For more examples and advanced usage, see the test files in the `pcms/pythonapi/` directory. diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 0ce22eba..f04ec029 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -54,15 +54,15 @@ void VerifyUniformGridFieldValues( // Helper function to verify mask field values void VerifyMaskFieldValues( const pcms::UniformGrid<2>& grid, - const std::vector& mask_field) { + const pcms::UniformGridField<2>& mask_field) { fprintf(stderr, "\nVerifying mask_field values:\n"); - auto mask_data = mask_field.data(); - for (int j = 0; j < grid.divisions[1]; ++j) { - for (int i = 0; i < grid.divisions[0]; ++i) { - int cell_id = j * grid.divisions[0] + i; - int mask_value = mask_data[cell_id]; - fprintf(stderr, "Cell (%d, %d): mask = %d\n", i, j, mask_value); - REQUIRE(mask_value == 1); + auto mask_data = mask_field.GetDOFHolderData(); + for (int j = 0; j <= grid.divisions[1]; ++j) { + for (int i = 0; i <= grid.divisions[0]; ++i) { + int vertex_id = j * (grid.divisions[0] + 1) + i; + pcms::Real mask_value = mask_data(vertex_id); + fprintf(stderr, "Vertex (%d, %d): mask = %.0f\n", i, j, mask_value); + REQUIRE(mask_value == 1.0); } } } @@ -371,20 +371,24 @@ TEST_CASE("Create binary field from uniform grid") auto lib = Omega_h::Library{}; auto world = lib.world(); - SECTION("Simple 2D box mesh - all cells inside") + SECTION("Simple 2D box mesh - all vertices inside") { // Create a mesh that fills the domain [0,1] x [0,1] auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, 10, 0, false); // Create a 5x5 grid (coarser than mesh) - auto field = CreateUniformGridBinaryField<2>(mesh, {5, 5}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {5, 5}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 25); + REQUIRE(field_data.extent(0) == 36); // (5+1) * (5+1) = 36 vertices - // All grid cells should be inside the mesh - int sum = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(sum == 25); + // All grid vertices should be inside the mesh + pcms::Real sum = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + sum += field_data(i); + } + REQUIRE(sum == 36.0); } SECTION("Binary field with custom divisions") @@ -392,14 +396,18 @@ TEST_CASE("Create binary field from uniform grid") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 8, 8, 0, false); - auto field = CreateUniformGridBinaryField<2>(mesh, {10, 8}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 8}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 80); + REQUIRE(field_data.extent(0) == 99); // (10+1) * (8+1) = 99 vertices - // Most cells should be inside - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); - REQUIRE(inside_count <= 80); + // Most vertices should be inside + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); + REQUIRE(inside_count <= 99.0); } SECTION("Binary field with equal divisions convenience function") @@ -407,12 +415,16 @@ TEST_CASE("Create binary field from uniform grid") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); - auto field = CreateUniformGridBinaryField<2>(mesh, 8); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 8); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 64); + REQUIRE(field_data.extent(0) == 81); // (8+1) * (8+1) = 81 vertices - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); } SECTION("Fine grid over coarse mesh") @@ -423,18 +435,19 @@ TEST_CASE("Create binary field from uniform grid") // Create a fine grid (20x20) auto grid = CreateUniformGridFromMesh<2>(mesh, {20, 20}); - auto field = CreateUniformGridBinaryField<2>(mesh, {20, 20}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {20, 20}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 400); + REQUIRE(field_data.extent(0) == 441); // (20+1) * (20+1) = 441 vertices - // Verify consistency: check some specific cells - // Center cells should be inside - auto center_idx = grid.GetCellIndex({10, 10}); - REQUIRE(field[center_idx] == 1); + // Verify consistency: check some specific vertices + // Center vertex should be inside (vertex at i=10, j=10) + int center_idx = 10 * 21 + 10; // 21 vertices per row + REQUIRE(field_data(center_idx) == 1.0); - // Corner cells should be inside - auto corner_idx = grid.GetCellIndex({0, 0}); - REQUIRE(field[corner_idx] == 1); + // Corner vertex should be inside + int corner_idx = 0; // vertex (0, 0) + REQUIRE(field_data(corner_idx) == 1.0); } SECTION("Verify field values are binary") @@ -442,15 +455,17 @@ TEST_CASE("Create binary field from uniform grid") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); - auto field = CreateUniformGridBinaryField<2>(mesh, 10); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10); + auto field_data = field->GetDOFHolderData(); // All values should be 0 or 1 - for (const auto& val : field) { - REQUIRE((val == 0 || val == 1)); + for (size_t i = 0; i < field_data.extent(0); ++i) { + pcms::Real val = field_data(i); + REQUIRE((val == 0.0 || val == 1.0)); } } - SECTION("Grid extends beyond mesh - cells outside should be marked 0") + SECTION("Grid extends beyond mesh - vertices outside should be marked 0") { // Create a small mesh in the center of a domain auto mesh = @@ -459,13 +474,17 @@ TEST_CASE("Create binary field from uniform grid") // The mesh occupies [0, 0.5] x [0, 0.5] // Create a grid that would cover this auto grid = CreateUniformGridFromMesh<2>(mesh, {10, 10}); - auto field = CreateUniformGridBinaryField<2>(mesh, {10, 10}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 10}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 100); + REQUIRE(field_data.extent(0) == 121); // (10+1) * (10+1) = 121 vertices - // All cells should be inside since grid is exactly on mesh bbox - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + // All vertices should be inside since grid is exactly on mesh bbox + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); } SECTION("Test with different aspect ratio") @@ -474,16 +493,20 @@ TEST_CASE("Create binary field from uniform grid") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 3.0, 1.0, 0.0, 12, 4, 0, false); - auto field = CreateUniformGridBinaryField<2>(mesh, {30, 10}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {30, 10}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 300); + REQUIRE(field_data.extent(0) == 341); // (30+1) * (10+1) = 341 vertices - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); // Calculate percentage inside - double inside_percent = 100.0 * inside_count / field.size(); - // Most cells should be inside + double inside_percent = 100.0 * inside_count / field_data.extent(0); + // Most vertices should be inside REQUIRE(inside_percent > 50.0); } } @@ -496,56 +519,73 @@ TEST_CASE("Binary field integration with grid methods") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 10, 10, 0, false); - SECTION("Query field value at specific grid cell") + SECTION("Query field value at specific grid vertex") { auto grid = CreateUniformGridFromMesh<2>(mesh, {8, 8}); - auto field = CreateUniformGridBinaryField<2>(mesh, {8, 8}); - - // Get field value for a specific cell - pcms::LO cell_id = grid.GetCellIndex({4, 4}); // Middle cell - REQUIRE(field[cell_id] == 1); // Should be inside - - // Get bbox of that cell - auto bbox = grid.GetCellBBOX(cell_id); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {8, 8}); + auto field_data = field->GetDOFHolderData(); + + // Get field value for a specific vertex (middle vertex at i=4, j=4) + pcms::LO vertex_id = 4 * 9 + 4; // 9 = (8+1) vertices per row + REQUIRE(field_data(vertex_id) == 1.0); // Should be inside + + // Compute vertex position + pcms::Real dx = grid.edge_length[0] / grid.divisions[0]; + pcms::Real dy = grid.edge_length[1] / grid.divisions[1]; + pcms::Real x = grid.bot_left[0] + 4 * dx; + pcms::Real y = grid.bot_left[1] + 4 * dy; + // Verify it's roughly in the middle - REQUIRE(bbox.center[0] == Catch::Approx(0.5).margin(0.1)); - REQUIRE(bbox.center[1] == Catch::Approx(0.5).margin(0.1)); + REQUIRE(x == Catch::Approx(0.5).margin(0.1)); + REQUIRE(y == Catch::Approx(0.5).margin(0.1)); } - SECTION("Count cells by region") + SECTION("Count vertices by region") { auto grid = CreateUniformGridFromMesh<2>(mesh, 10); - auto field = CreateUniformGridBinaryField<2>(mesh, 10); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10); + auto field_data = field->GetDOFHolderData(); - // Count cells in different quadrants + // Count vertices in different quadrants int q1 = 0, q2 = 0, q3 = 0, q4 = 0; // quadrants - for (pcms::LO i = 0; i < grid.GetNumCells(); ++i) { - if (field[i] == 1) { - auto bbox = grid.GetCellBBOX(i); - if (bbox.center[0] < 0.5 && bbox.center[1] < 0.5) - q1++; - else if (bbox.center[0] >= 0.5 && bbox.center[1] < 0.5) - q2++; - else if (bbox.center[0] < 0.5 && bbox.center[1] >= 0.5) - q3++; - else - q4++; + pcms::Real dx = grid.edge_length[0] / grid.divisions[0]; + pcms::Real dy = grid.edge_length[1] / grid.divisions[1]; + + for (int j = 0; j <= grid.divisions[1]; ++j) { + for (int i = 0; i <= grid.divisions[0]; ++i) { + pcms::LO vertex_id = j * (grid.divisions[0] + 1) + i; + if (field_data(vertex_id) == 1.0) { + pcms::Real x = grid.bot_left[0] + i * dx; + pcms::Real y = grid.bot_left[1] + j * dy; + + if (x < 0.5 && y < 0.5) + q1++; + else if (x >= 0.5 && y < 0.5) + q2++; + else if (x < 0.5 && y >= 0.5) + q3++; + else + q4++; + } } } - // All quadrants should have some inside cells + // All quadrants should have some inside vertices REQUIRE(q1 > 0); REQUIRE(q2 > 0); REQUIRE(q3 > 0); REQUIRE(q4 > 0); - // Distribution should be roughly equal + // For 11x11 vertices, boundary at x=0.5, y=0.5 splits asymmetrically: + // q1 (x<0.5, y<0.5): 5x5=25, q2 (x>=0.5, y<0.5): 6x5=30 + // q3 (x<0.5, y>=0.5): 5x6=30, q4 (x>=0.5, y>=0.5): 6x6=36 int total = q1 + q2 + q3 + q4; - REQUIRE(q1 == Catch::Approx(total / 4.0).margin(5)); - REQUIRE(q2 == Catch::Approx(total / 4.0).margin(5)); - REQUIRE(q3 == Catch::Approx(total / 4.0).margin(5)); - REQUIRE(q4 == Catch::Approx(total / 4.0).margin(5)); + REQUIRE(total == 121); // All vertices inside + REQUIRE(q1 == 25); + REQUIRE(q2 == 30); + REQUIRE(q3 == 30); + REQUIRE(q4 == 36); } } @@ -560,12 +600,16 @@ TEST_CASE("Performance and edge cases") Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 5, 5, 0, false); // Create a very fine grid - auto field = CreateUniformGridBinaryField<2>(mesh, 50); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 50); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 2500); + REQUIRE(field_data.extent(0) == 2601); // (50+1) * (50+1) = 2601 vertices - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); } SECTION("Coarse grid") @@ -574,13 +618,17 @@ TEST_CASE("Performance and edge cases") 10, 0, false); // Very coarse grid - auto field = CreateUniformGridBinaryField<2>(mesh, {2, 2}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {2, 2}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 4); + REQUIRE(field_data.extent(0) == 9); // (2+1) * (2+1) = 9 vertices - // All 4 cells should be inside for this configuration - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + // All vertices should be inside for this configuration + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); } SECTION("Non-square domain") @@ -588,15 +636,19 @@ TEST_CASE("Performance and edge cases") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 5.0, 2.0, 0.0, 20, 8, 0, false); - auto field = CreateUniformGridBinaryField<2>(mesh, {25, 10}); + auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {25, 10}); + auto field_data = field->GetDOFHolderData(); - REQUIRE(field.size() == 250); + REQUIRE(field_data.extent(0) == 286); // (25+1) * (10+1) = 286 vertices - int inside_count = std::accumulate(field.begin(), field.end(), 0); - REQUIRE(inside_count > 0); + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + REQUIRE(inside_count > 0.0); } - SECTION("Grid larger than mesh - cells outside marked as 0") + SECTION("Grid larger than mesh - vertices outside marked as 0") { // Create a mesh covering [0, 0.5] x [0, 0.5] auto mesh = @@ -609,40 +661,38 @@ TEST_CASE("Performance and edge cases") grid.divisions = {10, 10}; // Create binary field on the larger grid - auto field = pcms::CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid); - - REQUIRE(field.size() == 100); - - // Count cells inside and outside - int inside_count = std::accumulate(field.begin(), field.end(), 0); - int outside_count = field.size() - inside_count; - - // Should have both inside (1) and outside (0) cells - REQUIRE(inside_count > 0); - REQUIRE(outside_count > 0); - - // Check specific cells - // Bottom-left corner [0, 0.1] x [0, 0.1] should be inside (mesh covers [0, - // 0.5]) - auto corner_id = grid.GetCellIndex({0, 0}); - auto bbox = grid.GetCellBBOX(corner_id); - REQUIRE(field[corner_id] == 1); - - // Top-right corner [0.9, 1.0] x [0.9, 1.0] should be outside (mesh ends at - // 0.5) - corner_id = grid.GetCellIndex({9, 9}); - bbox = grid.GetCellBBOX(corner_id); - REQUIRE(field[corner_id] == 0); - - // Cell at [0.5, 0.6] x [0.5, 0.6] should be outside - corner_id = grid.GetCellIndex({6, 6}); - bbox = grid.GetCellBBOX(corner_id); - REQUIRE(field[corner_id] == 0); - - // Cell at [0.2, 0.3] x [0.2, 0.3] should be inside - auto center_id = grid.GetCellIndex({2, 2}); - bbox = grid.GetCellBBOX(center_id); - REQUIRE(field[center_id] == 1); + auto [layout, field] = pcms::CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid); + auto field_data = field->GetDOFHolderData(); + + REQUIRE(field_data.extent(0) == 121); // (10+1) * (10+1) = 121 vertices + + // Count vertices inside and outside + pcms::Real inside_count = 0.0; + for (size_t i = 0; i < field_data.extent(0); ++i) { + inside_count += field_data(i); + } + pcms::Real outside_count = field_data.extent(0) - inside_count; + + // Should have both inside (1) and outside (0) vertices + REQUIRE(inside_count > 0.0); + REQUIRE(outside_count > 0.0); + + // Check specific vertices + // Bottom-left corner should be inside (mesh covers [0, 0.5]) + int corner_id = 0 * 11 + 0; // vertex (0, 0) + REQUIRE(field_data(corner_id) == 1.0); + + // Top-right corner should be outside (mesh ends at 0.5) + corner_id = 10 * 11 + 10; // vertex (10, 10) + REQUIRE(field_data(corner_id) == 0.0); + + // Vertex at i=6, j=6 (coords ~0.6, ~0.6) should be outside + corner_id = 6 * 11 + 6; + REQUIRE(field_data(corner_id) == 0.0); + + // Vertex at i=2, j=2 (coords ~0.2, ~0.2) should be inside + auto center_id = 2 * 11 + 2; + REQUIRE(field_data(center_id) == 1.0); } } @@ -655,7 +705,7 @@ TEST_CASE("UniformGrid workflow") auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 4, 4, 0, false); auto grid = pcms::CreateUniformGridFromMesh<2>(mesh, {4, 4}); - auto mask_field = pcms::CreateUniformGridBinaryField<2>(mesh, {4, 4}); + auto [mask_layout, mask_field] = pcms::CreateUniformGridBinaryField<2>(mesh, {4, 4}); // Create OmegaH field layout with linear elements auto omega_h_layout = @@ -686,5 +736,5 @@ TEST_CASE("UniformGrid workflow") VerifyUniformGridFieldValues(grid, ug_coords, ug_field_data); // Verify mask field values - VerifyMaskFieldValues(grid, mask_field); + VerifyMaskFieldValues(grid, *mask_field); }