diff --git a/.github/workflows/clang-tidy.yml b/.github/workflows/clang-tidy.yml index 3204f6b1..d05507ea 100644 --- a/.github/workflows/clang-tidy.yml +++ b/.github/workflows/clang-tidy.yml @@ -63,7 +63,7 @@ jobs: repo-path: 'kokkos/kokkos' repo-ref: '4.6.01' cache: true - options: '-DCMAKE_CXX_STANDARD=17 + options: '-DCMAKE_CXX_STANDARD=20 -DBUILD_SHARED_LIBS=OFF -DKokkos_ENABLE_SERIAL=ON -DKokkos_ENABLE_OPENMP=OFF @@ -78,7 +78,7 @@ jobs: repo-path: 'kokkos/kokkos-kernels' repo-ref: '4.6.01' cache: true - options: '-DCMAKE_CXX_STANDARD=17 + options: '-DCMAKE_CXX_STANDARD=20 -DBUILD_SHARED_LIBS=OFF -DKokkos_DIR=${{ runner.temp }}/build-kokkos-openmpi/install/lib/cmake/Kokkos' diff --git a/.github/workflows/cmake-test.yml b/.github/workflows/cmake-test.yml index 591284a3..de0a6e97 100644 --- a/.github/workflows/cmake-test.yml +++ b/.github/workflows/cmake-test.yml @@ -25,6 +25,23 @@ jobs: compiler: [g++] language: ['cpp'] python_api: [OFF, ON] + meshfields: [ON] + petsc: [ON] + include: + - build_type: Release + memory_test: OFF + compiler: g++ + language: 'cpp' + python_api: OFF + meshfields: OFF + petsc: ON + - build_type: Release + memory_test: OFF + compiler: g++ + language: 'cpp' + python_api: OFF + meshfields: ON + petsc: OFF exclude: - build_type: Release memory_test: ON @@ -68,7 +85,7 @@ jobs: repo-ref: '4.6.01' cache: true cache-suffix: ${{ matrix.python_api == 'ON' && '-shared' || '' }} - options: '-DCMAKE_CXX_STANDARD=17 + options: '-DCMAKE_CXX_STANDARD=20 -DBUILD_SHARED_LIBS=${{ matrix.python_api }} -DKokkos_ENABLE_SERIAL=ON -DKokkos_ENABLE_OPENMP=OFF @@ -84,7 +101,7 @@ jobs: repo-ref: '4.6.01' cache: true cache-suffix: ${{ matrix.python_api == 'ON' && '-shared' || '' }} - options: '-DCMAKE_CXX_STANDARD=17 + options: '-DCMAKE_CXX_STANDARD=20 -DBUILD_SHARED_LIBS=${{ matrix.python_api }} -DKokkos_DIR=${{ runner.temp }}/build-kokkos/install/lib/cmake/Kokkos' @@ -106,6 +123,7 @@ jobs: -DKokkos_DIR=${{ runner.temp }}/build-kokkos/install/lib/cmake/Kokkos' - name: build meshFields + if: matrix.meshfields == 'ON' uses: ./.github/actions/install-repo with: repo-name: 'meshFields' @@ -154,6 +172,7 @@ jobs: -Dperfstubs_DIR=${{ runner.temp }}/build-perfstubs/install/lib/cmake' - name: clone petsc + if: matrix.petsc == 'ON' id: clone-petsc run: | cd ${{ runner.temp }} @@ -162,6 +181,7 @@ jobs: echo "petsc-commit-hash=$(git rev-parse HEAD)" >> $GITHUB_OUTPUT - name: Cache PETSc + if: matrix.petsc == 'ON' id: cache-petsc uses: actions/cache@v3 with: @@ -169,7 +189,7 @@ jobs: key: build-petsc-${{ steps.clone-petsc.outputs.petsc-commit-hash }} - name: build petsc - if: steps.cache-petsc.outputs.cache-hit != 'true' + if: matrix.petsc == 'ON' && steps.cache-petsc.outputs.cache-hit != 'true' run: | cd ${{ runner.temp }}/petsc ./configure \ @@ -181,6 +201,15 @@ jobs: --download-fblaslapack make all check + - name: Set PCMS PETSc options + run: | + echo "PCMS_ENABLE_PETSC=-DPCMS_ENABLE_PETSC=${{ matrix.petsc }}" >> $GITHUB_ENV + if [ "${{ matrix.petsc }}" = "ON" ]; then + echo "PCMS_PETSC_OPTIONS=-DPETSC_LINK_STATIC=ON -DPETSC_DIR=${{ runner.temp }}/petsc -DPETSC_ARCH=ubuntu-kokkos" >> $GITHUB_ENV + else + echo "PCMS_PETSC_OPTIONS=" >> $GITHUB_ENV + fi + - name: checkout pcms_testcases uses: actions/checkout@v3 with: @@ -193,7 +222,20 @@ jobs: - name: Set LD_LIBRARY_PATH for shared libraries if: matrix.python_api == 'ON' run: | - echo "LD_LIBRARY_PATH=${{ runner.temp }}/build-kokkos/install/lib:${{ runner.temp }}/build-omega_h/install/lib:${{ runner.temp }}/build-meshFields/install/lib:${{ runner.temp }}/build-redev/install/lib:${{ runner.temp }}/build-ADIOS2/install/lib:${{ runner.temp }}/build-perfstubs/install/lib:$LD_LIBRARY_PATH" >> $GITHUB_ENV + meshfields_lib="" + if [ "${{ matrix.meshfields }}" = "ON" ]; then + meshfields_lib="${{ runner.temp }}/build-meshFields/install/lib:" + fi + echo "LD_LIBRARY_PATH=${{ runner.temp }}/build-kokkos/install/lib:${{ runner.temp }}/build-omega_h/install/lib:${meshfields_lib}${{ runner.temp }}/build-redev/install/lib:${{ runner.temp }}/build-ADIOS2/install/lib:${{ runner.temp }}/build-perfstubs/install/lib:$LD_LIBRARY_PATH" >> $GITHUB_ENV + + - name: Set PCMS MeshFields options + run: | + echo "PCMS_ENABLE_MESHFIELDS=-DPCMS_ENABLE_MESHFIELDS=${{ matrix.meshfields }}" >> $GITHUB_ENV + if [ "${{ matrix.meshfields }}" = "ON" ]; then + echo "PCMS_MESHFIELDS_DIR=-Dmeshfields_DIR=${{ runner.temp }}/build-meshFields/install/lib/cmake/meshfields" >> $GITHUB_ENV + else + echo "PCMS_MESHFIELDS_DIR=" >> $GITHUB_ENV + fi - name: build pcms uses: ./.github/actions/install-repo @@ -207,16 +249,15 @@ jobs: -DCMAKE_CXX_COMPILER=`which mpicxx` -DCMAKE_Fortran_COMPILER=`which mpifort` -DMPIEXEC_PREFLAGS="--oversubscribe" - -DPCMS_ENABLE_PETSC=ON - -DPETSC_LINK_STATIC=ON - -DPETSC_DIR=${{ runner.temp }}/petsc - -DPETSC_ARCH=ubuntu-kokkos + ${{ env.PCMS_ENABLE_MESHFIELDS }} + ${{ env.PCMS_ENABLE_PETSC }} + ${{ env.PCMS_PETSC_OPTIONS }} -DPCMS_TIMEOUT=10 -DPCMS_ENABLE_SPDLOG=OFF -DPCMS_ENABLE_Python=${{ matrix.python_api }} -DCatch2_DIR=${{ runner.temp }}/build-Catch2/install/lib/cmake/Catch2 -DOmega_h_DIR=${{ runner.temp }}/build-omega_h/install/lib/cmake/Omega_h - -Dmeshfields_DIR=${{ runner.temp }}/build-meshFields/install/lib/cmake/meshfields + ${{ env.PCMS_MESHFIELDS_DIR }} -Dredev_DIR=${{ runner.temp }}/build-redev/install/lib/cmake/redev -DMPIEXEC_EXECUTABLE=`which mpirun` -DADIOS2_DIR=${{ runner.temp }}/build-ADIOS2/install/lib/cmake/adios2 @@ -264,8 +305,6 @@ jobs: -B ${{github.workspace}}/examples/external-usage-example/build \ -S ${{github.workspace}}/examples/external-usage-example/ \ -Dpcms_DIR=${{ runner.temp }}/build-pcms/install/lib/cmake/pcms \ - -DPETSC_LINK_STATIC=ON \ - -DPETSC_DIR=${{ runner.temp }}/petsc \ - -DPETSC_ARCH=ubuntu-kokkos \ + ${{ matrix.petsc == 'ON' && format('-DPETSC_LINK_STATIC=ON -DPETSC_DIR={0}/petsc -DPETSC_ARCH=ubuntu-kokkos', runner.temp) || '' }} \ --debug-output cmake --build ${{github.workspace}}/examples/external-usage-example/build diff --git a/.github/workflows/perlmutter/install.sh b/.github/workflows/perlmutter/install.sh index 3840592c..0b9ceb62 100644 --- a/.github/workflows/perlmutter/install.sh +++ b/.github/workflows/perlmutter/install.sh @@ -45,7 +45,7 @@ build-pcms/install # cmake -S kokkos-kernels -B build-kokkos-kernels \ # -DCMAKE_INSTALL_PREFIX=$PWD/build-kokkos-kernels/install \ # -DCMAKE_CXX_COMPILER=CC \ -# -DCMAKE_CXX_STANDARD=17 +# -DCMAKE_CXX_STANDARD=20 # cmake --build build-kokkos-kernels -j24 --target install @@ -131,4 +131,4 @@ cmake -S pcms -B build-pcms \ -DPCMS_TIMEOUT=100 \ -DCMAKE_CXX_STANDARD=20 \ -DPCMS_TEST_DATA_DIR=$PWD/pcms_testcases -cmake --build build-pcms -j8 \ No newline at end of file +cmake --build build-pcms -j8 diff --git a/README.md b/README.md index 656e48a3..cfd571df 100644 --- a/README.md +++ b/README.md @@ -32,7 +32,7 @@ module load cuda/12.1.1-zxa4msk git clone --branch 4.6.01 --depth 1 git@github.com:kokkos/kokkos.git cmake -S kokkos -B build-kokkos \ -DCMAKE_INSTALL_PREFIX=build-kokkos/install \ - -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ENABLE_SERIAL=ON \ -DKokkos_ENABLE_OPENMP=OFF \ -DKokkos_ENABLE_CUDA=OFF \ @@ -44,7 +44,7 @@ cmake --build build-kokkos --target install git clone --branch 4.6.01 --depth 1 git@github.com:kokkos/kokkos-kernels.git cmake -S kokkos-kernels -B build-kokkos-kernels \ -DCMAKE_INSTALL_PREFIX=build-kokkos-kernels/install \ - -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ROOT=$PWD/build-kokkos/install/lib64/cmake \ -DBUILD_SHARED_LIBS=OFF cmake --build build-kokkos-kernels --target install @@ -97,7 +97,7 @@ cmake --build build-redev --target install git clone --branch 4.6.01 --depth 1 git@github.com:kokkos/kokkos.git cmake -S kokkos -B build-kokkos \ -DCMAKE_INSTALL_PREFIX=build-kokkos/install \ - -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_STANDARD=20 \ -DCMAKE_BUILD_TYPE="Release" \ -DCMAKE_CXX_COMPILER=$PWD/kokkos/bin/nvcc_wrapper \ -DKokkos_ARCH_AMPERE80=ON \ @@ -112,7 +112,7 @@ cmake --build build-kokkos --target install git clone --branch 4.6.01 --depth 1 git@github.com:kokkos/kokkos-kernels.git cmake -S kokkos-kernels -B build-kokkos-kernels \ -DCMAKE_INSTALL_PREFIX=build-kokkos-kernels/install \ - -DCMAKE_CXX_STANDARD=17 \ + -DCMAKE_CXX_STANDARD=20 \ -DKokkos_ROOT=$PWD/build-kokkos/install/lib64/cmake \ -DBUILD_SHARED_LIBS=off cmake --build build-kokkos-kernels --target install diff --git a/examples/external-usage-example/main.cpp b/examples/external-usage-example/main.cpp index 8cc23e43..3af974b3 100644 --- a/examples/external-usage-example/main.cpp +++ b/examples/external-usage-example/main.cpp @@ -1,6 +1,5 @@ #include #include -#include #include int main() { diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2e18096b..c4edad3f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,7 +6,6 @@ set( pcms/coordinate_systems.h pcms/coordinate_transform.h pcms/field.h - pcms/create_field.h pcms/field_communicator.h pcms/field_communicator2.h pcms/field_evaluation_methods.h @@ -17,25 +16,23 @@ set( pcms/coordinate_system.h pcms/field_layout.h pcms/field_layout_communicator.h + pcms/nodal_field_factory.h ) set( PCMS_SOURCES pcms.cpp - pcms/create_field.cpp pcms/coupler2.cpp pcms/field_layout_communicator.cpp + pcms/nodal_field_factory.cpp pcms/adapter/point_cloud/point_cloud_layout.cpp pcms/adapter/point_cloud/point_cloud.cpp - pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp ) set( ADAPTER_HEADERS pcms/adapter/dummy_field_adapter.h pcms/adapter/point_cloud/point_cloud_layout.h pcms/adapter/point_cloud/point_cloud.h - pcms/adapter/meshfields/mesh_fields_adapter_layout.h - pcms/adapter/meshfields/mesh_fields_adapter2.h pcms/adapter/xgc/xgc_field_adapter.h ) @@ -52,21 +49,36 @@ if(PCMS_ENABLE_XGC) list(APPEND PCMS_SOURCES pcms/adapter/xgc/xgc_reverse_classification.cpp) list(APPEND ADAPTER_HEADERS pcms/adapter/xgc/xgc_reverse_classification.h) endif() +if(PCMS_ENABLE_MESHFIELDS) + list(APPEND PCMS_SOURCES + pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp + ) + list(APPEND ADAPTER_HEADERS + pcms/adapter/meshfields/mesh_fields_adapter_layout.h + pcms/adapter/meshfields/mesh_fields_adapter2.h + pcms/adapter/meshfields/mesh_fields_adapter.h + ) +endif() if(PCMS_ENABLE_OMEGA_H) list(APPEND PCMS_SOURCES + pcms/lagrange_field_factory.cpp + pcms/adapter/omega_h/omega_h_lagrange_layout.cpp pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp pcms/adapter/uniform_grid/uniform_grid_field.cpp ) list( APPEND PCMS_HEADERS + pcms/lagrange_field_factory.h pcms/transfer_field.h pcms/transfer_field2.h ) list(APPEND ADAPTER_HEADERS - pcms/adapter/meshfields/mesh_fields_adapter.h + pcms/adapter/omega_h/omega_h_lagrange_layout.h + pcms/adapter/omega_h/omega_h_lagrange_field.h pcms/adapter/uniform_grid/uniform_grid_field_layout.h pcms/adapter/uniform_grid/uniform_grid_field.h + pcms/adapter/uniform_grid/uniform_grid_binary_field.h ) endif() @@ -79,11 +91,14 @@ set_target_properties( core ) add_library(pcms::core ALIAS pcms_core) -target_compile_features(pcms_core PUBLIC cxx_std_17) +target_compile_features(pcms_core PUBLIC cxx_std_20) target_link_libraries( - pcms_core PUBLIC meshfields::meshfields redev::redev + pcms_core PUBLIC redev::redev MPI::MPI_CXX Kokkos::kokkos perfstubs pcms::utility pcms::localization ) +if(PCMS_ENABLE_MESHFIELDS) + target_link_libraries(pcms_core PUBLIC meshfields::meshfields) +endif() if(PCMS_ENABLE_OMEGA_H) target_link_libraries(pcms_core PUBLIC Omega_h::omega_h) endif() diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h index 291b47de..8a3e9543 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter2.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter2.h @@ -288,7 +288,7 @@ template class MeshFieldsAdapter2 : public FieldT { public: - MeshFieldsAdapter2(const MeshFieldsAdapterLayout& layout); + MeshFieldsAdapter2(std::shared_ptr layout); LocalizationHint GetLocalizationHint( CoordinateView coordinate_view) const override; @@ -316,7 +316,7 @@ class MeshFieldsAdapter2 : public FieldT ~MeshFieldsAdapter2() noexcept = default; private: - const MeshFieldsAdapterLayout& layout_; + std::shared_ptr layout_; Omega_h::Mesh& mesh_; std::unique_ptr> mesh_field_; GridPointSearch2D search_; @@ -328,16 +328,16 @@ class MeshFieldsAdapter2 : public FieldT */ template inline MeshFieldsAdapter2::MeshFieldsAdapter2( - const MeshFieldsAdapterLayout& layout) - : layout_(layout), - mesh_(layout.GetMesh()), + std::shared_ptr layout) + : layout_(std::move(layout)), + mesh_(layout_->GetMesh()), search_(mesh_, 10, 10), - dof_holder_data_("dof_holder_data", static_cast(layout.OwnedSize())) + dof_holder_data_("dof_holder_data", static_cast(layout_->OwnedSize())) { if (mesh_.dim() == 3) { throw pcms_error("MeshFieldsAdapter2 does not support 3D meshes"); } - auto nodes_per_dim = layout.GetNodesPerDim(); + auto nodes_per_dim = layout_->GetNodesPerDim(); if (nodes_per_dim[2] == 0 && nodes_per_dim[3] == 0) { if (nodes_per_dim[0] == 1 && nodes_per_dim[1] == 0) { switch (mesh_.dim()) { @@ -368,8 +368,8 @@ inline Rank1View MeshFieldsAdapter2::GetDOFHolderData() const { PCMS_FUNCTION_TIMER; - auto nodes_per_dim = layout_.GetNodesPerDim(); - auto num_components = layout_.GetNumComponents(); + auto nodes_per_dim = layout_->GetNodesPerDim(); + auto num_components = layout_->GetNumComponents(); size_t offset = 0; for (int i = 0; i <= mesh_.dim(); ++i) { if (nodes_per_dim[i]) { @@ -392,10 +392,10 @@ inline void MeshFieldsAdapter2::SetDOFHolderData( { PCMS_FUNCTION_TIMER; - auto nodes_per_dim = layout_.GetNodesPerDim(); - auto num_components = layout_.GetNumComponents(); + auto nodes_per_dim = layout_->GetNodesPerDim(); + auto num_components = layout_->GetNumComponents(); PCMS_ALWAYS_ASSERT(static_cast(data.size()) == - layout_.GetNumOwnedDofHolder() * num_components); + layout_->GetNumOwnedDofHolder() * num_components); size_t offset = 0; for (int i = 0; i <= mesh_.dim(); ++i) { if (nodes_per_dim[i]) { @@ -418,7 +418,7 @@ inline LocalizationHint MeshFieldsAdapter2::GetLocalizationHint( // TODO decide if we want to implicitly perform the coordinate transformations // when possible if (coordinate_view.GetCoordinateSystem() != - layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { throw pcms_error("Coordinate system mismatch"); } @@ -445,7 +445,7 @@ inline void MeshFieldsAdapter2::Evaluate( // TODO decide if we want to implicitly perform the coordinate transformations // when possible if (results.GetCoordinateSystem() != - layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { throw pcms_error("Coordinate system mismatch"); } @@ -491,7 +491,7 @@ inline void MeshFieldsAdapter2::EvaluateGradient( template inline const FieldLayout& MeshFieldsAdapter2::GetLayout() const { - return layout_; + return *layout_; } template @@ -510,7 +510,7 @@ inline int MeshFieldsAdapter2::Serialize( // host copy of filtered field data array const auto array_h = GetDOFHolderData(); if (buffer.size() > 0) { - auto owned = layout_.GetOwned(); + auto owned = layout_->GetOwned(); for (size_t i = 0; i < array_h.size(); i++) { if (owned[i]) buffer[permutation[i]] = array_h[i]; @@ -526,7 +526,7 @@ inline void MeshFieldsAdapter2::Deserialize( { PCMS_FUNCTION_TIMER; Omega_h::HostWrite sorted_buffer(permutation.size()); - auto owned = layout_.GetOwned(); + auto owned = layout_->GetOwned(); for (LO i = 0; i < sorted_buffer.size(); ++i) { if (owned[i]) sorted_buffer[i] = buffer[permutation[i]]; diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp b/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp index 39c0cb92..b8a6ee3e 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.cpp @@ -1,5 +1,4 @@ #include "mesh_fields_adapter.h" -#include "mesh_fields_adapter2.h" #include "pcms/adapter/meshfields/mesh_fields_adapter_layout.h" #include "mesh_fields_adapter_layout.h" #include "pcms/utility/inclusive_scan.h" @@ -198,11 +197,6 @@ MeshFieldsAdapterLayout::MeshFieldsAdapterLayout( gids_host_ = Omega_h::HostWrite(gids_); } -std::unique_ptr> MeshFieldsAdapterLayout::CreateFieldReal() const -{ - return std::make_unique>(*this); -} - int MeshFieldsAdapterLayout::GetNumComponents() const { return num_components_; diff --git a/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.h b/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.h index 79ad3cee..cc0c9e6c 100644 --- a/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.h +++ b/src/pcms/adapter/meshfields/mesh_fields_adapter_layout.h @@ -20,8 +20,6 @@ class MeshFieldsAdapterLayout : public FieldLayout CoordinateSystem coordinate_system, std::string global_id_name = "global"); - std::unique_ptr> CreateFieldReal() const override; - int GetNumComponents() const override; // nodes for standard lagrange FEM LO GetNumOwnedDofHolder() const override; diff --git a/src/pcms/adapter/omega_h/omega_h_lagrange_field.h b/src/pcms/adapter/omega_h/omega_h_lagrange_field.h new file mode 100644 index 00000000..decca851 --- /dev/null +++ b/src/pcms/adapter/omega_h/omega_h_lagrange_field.h @@ -0,0 +1,234 @@ +#ifndef PCMS_ADAPTER_OMEGA_H_LAGRANGE_FIELD_H +#define PCMS_ADAPTER_OMEGA_H_LAGRANGE_FIELD_H + +#include "pcms/adapter/omega_h/omega_h_lagrange_layout.h" +#include "pcms/field.h" +#include "pcms/localization/point_search.h" +#include "pcms/utility/assert.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/profile.h" + +#include +#include + +namespace pcms { +struct OmegaHLagrangeLocHint { + // Per valid query point: + Kokkos::View elem_ids; // containing element + Kokkos::View bary; // [n_valid x (dim+1)] + Kokkos::View orig_indices; // original query index + // Points that fell outside the mesh: + Kokkos::View missing_indices; + OutOfBoundsMode mode; +}; + +namespace detail { +template OmegaHLagrangeLocHint BuildLagrangeLocHint( + int mesh_dim, + Kokkos::View::Result *, + HostMemorySpace> results_h, + OutOfBoundsMode mode) { + LO n = static_cast(results_h.size()); + std::vector valid, missing; + for (LO i = 0; i < n; ++i) { + bool out = (static_cast(results_h(i).dimensionality) != mesh_dim) || + (results_h(i).element_id < 0); + if (out) + missing.push_back(i); + else + valid.push_back(i); + } + + if (mode == OutOfBoundsMode::ERROR) { + PCMS_ALWAYS_ASSERT(missing.empty() && + "Points found outside mesh domain"); + } else if (mode == OutOfBoundsMode::NEAREST_BOUNDARY) { + PCMS_ALWAYS_ASSERT(false && + "NEAREST_BOUNDARY mode not yet implemented"); + } + + LO nv = static_cast(valid.size()); + LO nm = static_cast(missing.size()); + + Kokkos::View elem_ids("elem_ids", nv); + Kokkos::View bary("bary", nv, mesh_dim + 1); + Kokkos::View orig_indices("orig_indices", nv); + Kokkos::View missing_indices("missing_indices", nm); + + for (LO k = 0; k < nv; ++k) { + LO i = valid[k]; + elem_ids(k) = results_h(i).element_id; + orig_indices(k) = i; + for (int d = 0; d <= mesh_dim; ++d) + bary(k, d) = results_h(i).parametric_coords[d]; + } + for (LO k = 0; k < nm; ++k) + missing_indices(k) = missing[k]; + + return OmegaHLagrangeLocHint{ + elem_ids, bary, orig_indices, missing_indices, + mode + }; +} +} // namespace detail + +namespace detail { +inline std::variant MakeSearch( + Omega_h::Mesh &mesh) { + if (mesh.dim() == 2) + return GridPointSearch2D(mesh, 10, 10); + else if (mesh.dim() == 3) + return GridPointSearch3D(mesh, 10, 10, 10); + throw std::invalid_argument( + "OmegaHLagrangeField: only 2D and 3D meshes are supported"); +} +} // namespace detail + +template class OmegaHLagrangeField : public FieldT { +public: + explicit OmegaHLagrangeField( + std::shared_ptr layout) + : layout_(std::move(layout)), + mesh_(layout_->GetMesh()), + search_(detail::MakeSearch(mesh_)), + dof_holder_data_("dof_holder_data", + static_cast(layout_->OwnedSize())) { + } + + LocalizationHint GetLocalizationHint( + CoordinateView coordinate_view) const override { + PCMS_FUNCTION_TIMER; + if (coordinate_view.GetCoordinateSystem() != + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { + throw pcms_error("Coordinate system mismatch"); + } + + auto coords = coordinate_view.GetCoordinates(); + LO n_pts = static_cast(coords.extent(0)); + int mesh_dim = mesh_.dim(); + + auto hint = + std::make_shared(std::visit( + [&](auto &search) { + using SearchT = std::decay_t; + constexpr int Dim = SearchT::DIM; + + Kokkos::View coords_d("coords_d", n_pts); + auto coords_h = Kokkos::View( + coords.data_handle(), + n_pts, + Dim); + deep_copy_mismatch_layouts(coords_d, coords_h); + + auto results_d = search(coords_d); + Kokkos::View::Result *, + HostMemorySpace> + results_h("results_h", results_d.size()); + Kokkos::deep_copy(results_h, results_d); + + return detail::BuildLagrangeLocHint(mesh_dim, + results_h, + this->out_of_bounds_mode_); + }, + search_)); + + return LocalizationHint{hint}; + } + + void Evaluate(LocalizationHint location, + FieldDataView results) const override { + PCMS_FUNCTION_TIMER; + if (results.GetCoordinateSystem() != + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { + throw pcms_error("Coordinate system mismatch"); + } + + const auto &hint = + *reinterpret_cast(location.data.get()); + auto values = results.GetValues(); + LO n_valid = static_cast(hint.elem_ids.size()); + + if (layout_->GetOrder() == 0) { + // One DOF per element — return the element's value directly + for (LO k = 0; k < n_valid; ++k) + values[hint.orig_indices(k)] = dof_holder_data_(hint.elem_ids(k)); + } else { + // Order-1: barycentric interpolation over element vertices + int mesh_dim = mesh_.dim(); + int nvpe = mesh_dim + 1; // vertices per element (simplex) + auto elem_verts = Omega_h::HostRead(mesh_.ask_elem_verts()); + for (LO k = 0; k < n_valid; ++k) { + LO elem = hint.elem_ids(k); + T val = T{}; + for (int v = 0; v < nvpe; ++v) + val += static_cast(hint.bary(k, v)) * + dof_holder_data_(elem_verts[elem * nvpe + v]); + values[hint.orig_indices(k)] = val; + } + } + + if (hint.mode == OutOfBoundsMode::FILL) { + T fill_val = static_cast(this->fill_value_); + for (LO k = 0; k < static_cast(hint.missing_indices.size()); ++k) + values[hint.missing_indices(k)] = fill_val; + } + } + + void EvaluateGradient(FieldDataView /*unused*/) override { + throw pcms_error( + "EvaluateGradient not implemented for OmegaHLagrangeField"); + } + + bool CanEvaluateGradient() override { return false; } + + Rank1View GetDOFHolderData() const override { + return make_const_array_view(dof_holder_data_); + } + + void SetDOFHolderData(Rank1View data) override { + PCMS_ALWAYS_ASSERT(static_cast(data.size()) == + layout_->GetNumOwnedDofHolder() * + layout_->GetNumComponents()); + for (size_t i = 0; i < data.size(); ++i) + dof_holder_data_(i) = data[i]; + } + + const FieldLayout &GetLayout() const override { return *layout_; } + + int Serialize( + Rank1View buffer, + Rank1View permutation) const override { + PCMS_FUNCTION_TIMER; + auto owned = layout_->GetOwned(); + LO n = static_cast(dof_holder_data_.size()); + if (buffer.size() > 0) { + for (LO i = 0; i < n; ++i) { + if (owned[i]) + buffer[permutation[i]] = dof_holder_data_(i); + } + } + return n; + } + + void Deserialize( + Rank1View buffer, + Rank1View permutation) override { + PCMS_FUNCTION_TIMER; + auto owned = layout_->GetOwned(); + LO n = static_cast(dof_holder_data_.size()); + for (LO i = 0; i < n; ++i) { + if (owned[i]) + dof_holder_data_(i) = buffer[permutation[i]]; + } + } + + ~OmegaHLagrangeField() noexcept = default; + +private: + std::shared_ptr layout_; + Omega_h::Mesh &mesh_; + std::variant search_; + Kokkos::View dof_holder_data_; +}; +} // namespace pcms +#endif // PCMS_ADAPTER_OMEGA_H_LAGRANGE_FIELD_H diff --git a/src/pcms/adapter/omega_h/omega_h_lagrange_layout.cpp b/src/pcms/adapter/omega_h/omega_h_lagrange_layout.cpp new file mode 100644 index 00000000..b6fbec51 --- /dev/null +++ b/src/pcms/adapter/omega_h/omega_h_lagrange_layout.cpp @@ -0,0 +1,182 @@ +#include "pcms/adapter/omega_h/omega_h_lagrange_layout.h" +#include "pcms/partition.h" +#include "pcms/utility/assert.h" +#include "pcms/utility/mesh_geometry.h" +#include "pcms/utility/profile.h" +#include + +namespace pcms +{ + +namespace +{ +int EntityDimForOrder(int order, int mesh_dim) +{ + switch (order) { + case 0: return mesh_dim; // one DOF per element + case 1: return 0; // one DOF per vertex + default: + throw std::invalid_argument( + "OmegaHLagrangeLayout: only order 0 and 1 are supported"); + } +} + +template +Omega_h::Write GetGidsHelper(Omega_h::Mesh& mesh, int entity_dim, + const std::string& global_id_name) +{ + auto dim_gids = mesh.get_array(entity_dim, global_id_name); + Omega_h::Write gids(dim_gids.size()); + Omega_h::parallel_for( + dim_gids.size(), + OMEGA_H_LAMBDA(int i) { gids[i] = dim_gids[i]; }); + return gids; +} + +Omega_h::HostWrite BuildGids(Omega_h::Mesh& mesh, int entity_dim, + const std::string& global_id_name) +{ + auto tag = mesh.get_tagbase(entity_dim, global_id_name); + Omega_h::Write gids; + if (Omega_h::is(tag)) { + gids = GetGidsHelper(mesh, entity_dim, global_id_name); + } else if (Omega_h::is(tag)) { + gids = GetGidsHelper(mesh, entity_dim, global_id_name); + } else { + std::cerr << "Weird tag type for global arrays.\n"; + std::abort(); + } + return Omega_h::HostWrite(gids); +} + +Kokkos::View BuildOwned(Omega_h::Mesh& mesh, + int entity_dim) +{ + int n = mesh.nents(entity_dim); + Kokkos::View owned("owned", n); + auto src = Omega_h::HostRead(mesh.owned(entity_dim)); + for (int i = 0; i < n; ++i) + owned(i) = static_cast(src[i]); + return owned; +} +} // namespace + +OmegaHLagrangeLayout::OmegaHLagrangeLayout(Omega_h::Mesh& mesh, int order, + int num_components, + CoordinateSystem coordinate_system, + std::string global_id_name) + : mesh_(mesh), + order_(order), + num_components_(num_components), + coordinate_system_(coordinate_system), + global_id_name_(std::move(global_id_name)) +{ + PCMS_FUNCTION_TIMER; + int entity_dim = EntityDimForOrder(order_, mesh_.dim()); + + gids_host_ = BuildGids(mesh_, entity_dim, global_id_name_); + coords_host_ = + Omega_h::HostRead(get_entity_centroids(mesh_, entity_dim)); + owned_host_ = BuildOwned(mesh_, entity_dim); + + class_ids_host_ = Omega_h::HostRead( + mesh_.get_array(entity_dim, "class_id")); + class_dims_host_ = Omega_h::HostRead( + mesh_.get_array(entity_dim, "class_dim")); +} + +int OmegaHLagrangeLayout::GetNumComponents() const +{ + return num_components_; +} + +LO OmegaHLagrangeLayout::GetNumOwnedDofHolder() const +{ + return mesh_.nents(EntityDimForOrder(order_, mesh_.dim())); +} + +GO OmegaHLagrangeLayout::GetNumGlobalDofHolder() const +{ + return mesh_.nglobal_ents(EntityDimForOrder(order_, mesh_.dim())); +} + +Rank1View OmegaHLagrangeLayout::GetOwned() const +{ + return make_const_array_view(owned_host_); +} + +GlobalIDView OmegaHLagrangeLayout::GetGids() const +{ + return GlobalIDView(gids_host_.data(), gids_host_.size()); +} + +CoordinateView +OmegaHLagrangeLayout::GetDOFHolderCoordinates() const +{ + int n = mesh_.nents(EntityDimForOrder(order_, mesh_.dim())); + int dim = mesh_.dim(); + Rank2View coords_view(coords_host_.data(), n, + dim); + return CoordinateView{coordinate_system_, coords_view}; +} + +bool OmegaHLagrangeLayout::IsDistributed() +{ + return true; +} + +EntOffsetsArray OmegaHLagrangeLayout::GetEntOffsets() const +{ + // Slot i holds the starting DOF index for entity dimension i. + // Slot 4 holds the total DOF count. + // For order-1 (entity_dim=0): offsets = {0, n, n, n, n} + // For order-0 (entity_dim=mesh.dim()): offsets = {0,..,0, n, n} + EntOffsetsArray offsets{}; + offsets.fill(0); + int entity_dim = EntityDimForOrder(order_, mesh_.dim()); + LO n = mesh_.nents(entity_dim); + for (int i = entity_dim + 1; i < ent_offsets_len; ++i) + offsets[i] = n; + return offsets; +} + +ReversePartitionMap2 OmegaHLagrangeLayout::GetReversePartitionMap( + const redev::Partition& partition) const +{ + PCMS_FUNCTION_TIMER; + int entity_dim = EntityDimForOrder(order_, mesh_.dim()); + int n = mesh_.nents(entity_dim); + int dim = mesh_.dim(); + + ReversePartitionMap2 reverse_partition; + for (int i = 0; i < n; ++i) { + if (!owned_host_(i)) + continue; + + std::array coord{}; + for (int d = 0; d < dim; ++d) + coord[d] = coords_host_[i * dim + d]; + + auto dr = std::visit( + GetRank{class_ids_host_[i], class_dims_host_[i], coord}, partition); + reverse_partition[dr].indices.emplace_back(i); + + const auto n_slots = reverse_partition[dr].ent_offsets.size(); + for (size_t e = static_cast(entity_dim) + 1; e < n_slots; ++e) + reverse_partition[dr].ent_offsets[e] += 1; + } + + return reverse_partition; +} + +int OmegaHLagrangeLayout::GetOrder() const +{ + return order_; +} + +Omega_h::Mesh& OmegaHLagrangeLayout::GetMesh() const +{ + return mesh_; +} + +} // namespace pcms diff --git a/src/pcms/adapter/omega_h/omega_h_lagrange_layout.h b/src/pcms/adapter/omega_h/omega_h_lagrange_layout.h new file mode 100644 index 00000000..d6ad313c --- /dev/null +++ b/src/pcms/adapter/omega_h/omega_h_lagrange_layout.h @@ -0,0 +1,57 @@ +#ifndef PCMS_ADAPTER_OMEGA_H_LAGRANGE_LAYOUT_H +#define PCMS_ADAPTER_OMEGA_H_LAGRANGE_LAYOUT_H + +#include +#include "pcms/utility/arrays.h" +#include "pcms/field_layout.h" +#include "pcms/coordinate_system.h" +#include "pcms/field.h" + +namespace pcms +{ + +// Layout for native Omega_h Lagrange fields. +// order 0: one DOF holder per element (centroid coordinates) +// order 1: one DOF holder per vertex (barycentric interpolation) +// Throws std::invalid_argument for any other order. +class OmegaHLagrangeLayout : public FieldLayout +{ +public: + OmegaHLagrangeLayout(Omega_h::Mesh& mesh, int order, int num_components, + CoordinateSystem coordinate_system, + std::string global_id_name = "global"); + + 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; + + int GetOrder() const; + Omega_h::Mesh& GetMesh() const; + +private: + Omega_h::Mesh& mesh_; + int order_; + int num_components_; + CoordinateSystem coordinate_system_; + std::string global_id_name_; + + Omega_h::HostWrite gids_host_; + Omega_h::HostRead coords_host_; // flat row-major: [n_dofs * mesh_dim] + Kokkos::View owned_host_; + Omega_h::HostRead class_ids_host_; + Omega_h::HostRead class_dims_host_; +}; + +} // namespace pcms +#endif // PCMS_ADAPTER_OMEGA_H_LAGRANGE_LAYOUT_H diff --git a/src/pcms/adapter/point_cloud/point_cloud.cpp b/src/pcms/adapter/point_cloud/point_cloud.cpp index c0de9102..2c149dbb 100644 --- a/src/pcms/adapter/point_cloud/point_cloud.cpp +++ b/src/pcms/adapter/point_cloud/point_cloud.cpp @@ -39,10 +39,10 @@ struct PointCloudLocalizationHint Kokkos::View coordinates_; }; -PointCloud::PointCloud(const PointCloudLayout& layout) - : layout_(layout), - data_("", layout_.GetDOFHolderCoordinates().GetCoordinates().extent(0)), - data_host_("", layout_.GetDOFHolderCoordinates().GetCoordinates().extent(0)) +PointCloud::PointCloud(std::shared_ptr layout) + : layout_(std::move(layout)), + data_("", layout_->GetDOFHolderCoordinates().GetCoordinates().extent(0)), + data_host_("", layout_->GetDOFHolderCoordinates().GetCoordinates().extent(0)) { } @@ -84,7 +84,7 @@ void PointCloud::EvaluateGradient( const FieldLayout& PointCloud::GetLayout() const { - return layout_; + return *layout_; } bool PointCloud::CanEvaluateGradient() diff --git a/src/pcms/adapter/point_cloud/point_cloud.h b/src/pcms/adapter/point_cloud/point_cloud.h index f19298c4..638bb7b8 100644 --- a/src/pcms/adapter/point_cloud/point_cloud.h +++ b/src/pcms/adapter/point_cloud/point_cloud.h @@ -4,13 +4,14 @@ #include "pcms/field.h" #include "pcms/utility/arrays.h" #include "point_cloud_layout.h" +#include namespace pcms { class PointCloud : public FieldT { public: - PointCloud(const PointCloudLayout& layout); + PointCloud(std::shared_ptr layout); LocalizationHint GetLocalizationHint( CoordinateView coordinate_view) const override; @@ -36,7 +37,7 @@ class PointCloud : public FieldT void SetDOFHolderData(Rank1View data) override; private: - const PointCloudLayout& layout_; + std::shared_ptr layout_; Kokkos::View data_; Kokkos::View data_host_; }; diff --git a/src/pcms/adapter/point_cloud/point_cloud_layout.cpp b/src/pcms/adapter/point_cloud/point_cloud_layout.cpp index 93c7754c..04bcc5d8 100644 --- a/src/pcms/adapter/point_cloud/point_cloud_layout.cpp +++ b/src/pcms/adapter/point_cloud/point_cloud_layout.cpp @@ -1,5 +1,4 @@ #include "point_cloud_layout.h" -#include "point_cloud.h" #include #include @@ -23,11 +22,6 @@ PointCloudLayout::PointCloudLayout(int dim, Kokkos::View coords, iota_view(gids_); } -std::unique_ptr> PointCloudLayout::CreateFieldReal() const -{ - return std::make_unique(*this); -} - int PointCloudLayout::GetNumComponents() const { return components_; diff --git a/src/pcms/adapter/point_cloud/point_cloud_layout.h b/src/pcms/adapter/point_cloud/point_cloud_layout.h index 13897e20..c902682b 100644 --- a/src/pcms/adapter/point_cloud/point_cloud_layout.h +++ b/src/pcms/adapter/point_cloud/point_cloud_layout.h @@ -12,8 +12,6 @@ class PointCloudLayout : public FieldLayout PointCloudLayout(int dim, Kokkos::View coords, CoordinateSystem coordinate_system); - std::unique_ptr> CreateFieldReal() const override; - int GetNumComponents() const override; // nodes for standard lagrange FEM LO GetNumOwnedDofHolder() const override; diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_binary_field.h b/src/pcms/adapter/uniform_grid/uniform_grid_binary_field.h new file mode 100644 index 00000000..82397eb2 --- /dev/null +++ b/src/pcms/adapter/uniform_grid/uniform_grid_binary_field.h @@ -0,0 +1,92 @@ +#ifndef PCMS_UNIFORM_GRID_BINARY_FIELD_H +#define PCMS_UNIFORM_GRID_BINARY_FIELD_H + +#include "pcms/adapter/uniform_grid/uniform_grid_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/localization/point_search.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" +#include "pcms/utility/uniform_grid.h" + +#include +#include +#include +#include +#include + +namespace pcms +{ + +/** + * \brief Create a binary (inside/outside) mask field on a uniform grid. + * + * Each DOF in the returned order-1 (vertex-centered) UniformGridField is set + * to 1.0 if the vertex lies inside the mesh and 0.0 otherwise. A single + * point search pass over the grid vertices is performed; no intermediate + * source field is required. + * + * \tparam Dim Spatial dimension (2 or 3). + * \param mesh The Omega_h mesh defining the domain. + * \param grid Uniform grid whose vertices are tested. + * \return Pair of (layout, field) with binary mask values. + */ +template +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid& grid) +{ + auto layout = std::make_shared>( + grid, 1, CoordinateSystem::Cartesian); + auto field = std::make_unique>(layout); + + auto coord_view = layout->GetDOFHolderCoordinates(); + auto coords = coord_view.GetCoordinates(); + LO n = layout->GetNumOwnedDofHolder(); + + // Copy DOF coordinates from host mdspan to a Kokkos device view for search + Kokkos::View coords_d("coords_d", n); + auto coords_h = Kokkos::create_mirror_view(coords_d); + for (LO i = 0; i < n; ++i) + for (unsigned d = 0; d < Dim; ++d) + coords_h(i, d) = coords(i, d); + Kokkos::deep_copy(coords_d, coords_h); + + // Run point-in-mesh search + Kokkos::View::Result*> results_d; + if constexpr (Dim == 2) { + GridPointSearch2D search(mesh, grid.divisions[0], grid.divisions[1]); + results_d = search(coords_d); + } else { + GridPointSearch3D search(mesh, grid.divisions[0], grid.divisions[1], + grid.divisions[2]); + results_d = search(coords_d); + } + auto results_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, results_d); + + Kokkos::View data("binary_mask", n); + for (LO i = 0; i < n; ++i) + data(i) = (results_h(i).element_id >= 0) ? 1.0 : 0.0; + + field->SetDOFHolderData( + Rank1View(data.data(), n)); + + return {std::move(layout), std::move(field)}; +} + +/** + * \brief Convenience overload: derives the grid from the mesh bounding box. + */ +template +std::pair>, + std::unique_ptr>> +CreateUniformGridBinaryField(Omega_h::Mesh& mesh, + const std::array& divisions) +{ + return CreateUniformGridBinaryField( + mesh, CreateUniformGridFromMesh(mesh, divisions)); +} + +} // namespace pcms + +#endif // PCMS_UNIFORM_GRID_BINARY_FIELD_H diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp index ea5cb4be..37d0027c 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.cpp @@ -37,10 +37,10 @@ struct UniformGridFieldLocalizationHint template UniformGridField::UniformGridField( - const UniformGridFieldLayout& layout) - : layout_(layout), - grid_(layout.GetGrid()), - dof_holder_data_("dof_holder_data", static_cast(layout.OwnedSize())) + std::shared_ptr> layout) + : layout_(std::move(layout)), + grid_(layout_->GetGrid()), + dof_holder_data_("dof_holder_data", static_cast(layout_->OwnedSize())) { PCMS_FUNCTION_TIMER; // Default to NEAREST_BOUNDARY for uniform grid fields @@ -73,9 +73,18 @@ View UniformGridField::to_mdspan() PCMS_FUNCTION_TIMER; if constexpr (Dim == 2) { + if (layout_->GetOrder() == 0) { + return View( + dof_holder_data_.data(), grid_.divisions[0], grid_.divisions[1]); + } return View( dof_holder_data_.data(), grid_.divisions[0] + 1, grid_.divisions[1] + 1); } else if constexpr (Dim == 3) { + if (layout_->GetOrder() == 0) { + return View( + dof_holder_data_.data(), grid_.divisions[0], grid_.divisions[1], + grid_.divisions[2]); + } return View( dof_holder_data_.data(), grid_.divisions[0] + 1, grid_.divisions[1] + 1, grid_.divisions[2] + 1); @@ -91,9 +100,18 @@ View UniformGridField::to_mdspan() const PCMS_FUNCTION_TIMER; if constexpr (Dim == 2) { + if (layout_->GetOrder() == 0) { + return View( + dof_holder_data_.data(), grid_.divisions[0], grid_.divisions[1]); + } return View( dof_holder_data_.data(), grid_.divisions[0] + 1, grid_.divisions[1] + 1); } else if constexpr (Dim == 3) { + if (layout_->GetOrder() == 0) { + return View( + dof_holder_data_.data(), grid_.divisions[0], grid_.divisions[1], + grid_.divisions[2]); + } return View( dof_holder_data_.data(), grid_.divisions[0] + 1, grid_.divisions[1] + 1, grid_.divisions[2] + 1); @@ -110,7 +128,7 @@ LocalizationHint UniformGridField::GetLocalizationHint( PCMS_FUNCTION_TIMER; if (coordinate_view.GetCoordinateSystem() != - layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { throw std::runtime_error( "Coordinate system mismatch in GetLocalizationHint"); } @@ -160,7 +178,7 @@ void UniformGridField::Evaluate( PCMS_FUNCTION_TIMER; if (results.GetCoordinateSystem() != - layout_.GetDOFHolderCoordinates().GetCoordinateSystem()) { + layout_->GetDOFHolderCoordinates().GetCoordinateSystem()) { throw std::runtime_error("Coordinate system mismatch in Evaluate"); } @@ -173,6 +191,18 @@ void UniformGridField::Evaluate( auto cell_indices = hint.cell_indices_; LO num_points = coordinates.extent(0); + if (layout_->GetOrder() == 0) { + auto evaluated_values = results.GetValues(); + for (LO i = 0; i < evaluated_values.size(); ++i) { + if (hint.is_out_of_bounds_[i] && hint.mode_ == OutOfBoundsMode::FILL) { + evaluated_values[i] = fill_value_; + } else { + evaluated_values[i] = values[CellIdToDofIndex(cell_indices[i])]; + } + } + return; + } + Kokkos::View cell_dimensioned_indices( "cell_dimensioned_indices", num_points, Dim); // Convert cell indices to multi-dimensional indices @@ -184,7 +214,7 @@ void UniformGridField::Evaluate( } // Dimensions for vertex grid (m+1 vertices per dimension for m cells) - auto cell_divisions = layout_.GetGrid().divisions; + 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) { @@ -252,7 +282,7 @@ void UniformGridField::EvaluateGradient( template const FieldLayout& UniformGridField::GetLayout() const { - return layout_; + return *layout_; } template @@ -261,6 +291,12 @@ bool UniformGridField::CanEvaluateGradient() return false; } +template +LO UniformGridField::CellIdToDofIndex(LO cell_id) const +{ + return cell_id; +} + template int UniformGridField::Serialize( Rank1View buffer, @@ -287,7 +323,7 @@ void UniformGridField::Deserialize( Kokkos::View sorted_buffer("sorted_buffer", permutation.size()); - auto owned = layout_.GetOwned(); + auto owned = layout_->GetOwned(); for (LO i = 0; i < sorted_buffer.size(); ++i) { PCMS_ALWAYS_ASSERT(owned[i]); diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field.h b/src/pcms/adapter/uniform_grid/uniform_grid_field.h index c1fed28b..62bb6e96 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field.h +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field.h @@ -14,7 +14,7 @@ template class UniformGridField : public FieldT { public: - UniformGridField(const UniformGridFieldLayout& layout); + UniformGridField(std::shared_ptr> layout); LocalizationHint GetLocalizationHint( CoordinateView coordinate_view) const override; @@ -45,8 +45,10 @@ class UniformGridField : public FieldT ~UniformGridField() noexcept = default; private: - const UniformGridFieldLayout& layout_; - UniformGrid& grid_; + LO CellIdToDofIndex(LO cell_id) const; + + std::shared_ptr> layout_; + const UniformGrid& grid_; Kokkos::View dof_holder_data_; }; diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp index 93b26e18..6df8d953 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.cpp @@ -1,6 +1,6 @@ #include "uniform_grid_field_layout.h" -#include "uniform_grid_field.h" #include "pcms/utility/profile.h" +#include "pcms/utility/assert.h" #include namespace pcms @@ -8,67 +8,80 @@ namespace pcms template UniformGridFieldLayout::UniformGridFieldLayout( - UniformGrid& grid, int num_components, - CoordinateSystem coordinate_system) - : grid_(grid), + UniformGrid grid, int num_components, + CoordinateSystem coordinate_system, int order) + : grid_(std::move(grid)), num_components_(num_components), coordinate_system_(coordinate_system), - gids_("gids", GetNumVertices()), - dof_holder_coords_("dof_holder_coords", GetNumVertices(), Dim), - owned_("owned", GetNumVertices()) + order_(order), + gids_("gids", GetNumDofHolders()), + dof_holder_coords_("dof_holder_coords", GetNumDofHolders(), Dim), + owned_("owned", GetNumDofHolders()) { PCMS_FUNCTION_TIMER; + PCMS_ALWAYS_ASSERT(order_ == 0 || order_ == 1); - LO num_vertices = GetNumVertices(); + LO num_dofs = GetNumDofHolders(); // Initialize global IDs and ownership - for (LO i = 0; i < num_vertices; ++i) { + for (LO i = 0; i < num_dofs; ++i) { gids_[i] = static_cast(i); owned_[i] = true; } - // Initialize DOF holder coordinates at grid vertices + // Initialize DOF holder coordinates at grid vertices or cell centers. 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; + LO dof_idx = 0; + if (order_ == 1) { + for (LO j = 0; j <= grid_.divisions[1]; ++j) { + for (LO i = 0; i <= grid_.divisions[0]; ++i) { + dof_holder_coords_(dof_idx, 0) = grid_.bot_left[0] + i * vertex_spacing[0]; + dof_holder_coords_(dof_idx, 1) = grid_.bot_left[1] + j * vertex_spacing[1]; + ++dof_idx; + } + } + } else { + for (LO j = 0; j < grid_.divisions[1]; ++j) { + for (LO i = 0; i < grid_.divisions[0]; ++i) { + dof_holder_coords_(dof_idx, 0) = grid_.bot_left[0] + (i + 0.5) * vertex_spacing[0]; + dof_holder_coords_(dof_idx, 1) = grid_.bot_left[1] + (j + 0.5) * vertex_spacing[1]; + ++dof_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; + LO dof_idx = 0; + if (order_ == 1) { + 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_(dof_idx, 0) = grid_.bot_left[0] + i * vertex_spacing[0]; + dof_holder_coords_(dof_idx, 1) = grid_.bot_left[1] + j * vertex_spacing[1]; + dof_holder_coords_(dof_idx, 2) = grid_.bot_left[2] + k * vertex_spacing[2]; + ++dof_idx; + } + } + } + } else { + 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_(dof_idx, 0) = grid_.bot_left[0] + (i + 0.5) * vertex_spacing[0]; + dof_holder_coords_(dof_idx, 1) = grid_.bot_left[1] + (j + 0.5) * vertex_spacing[1]; + dof_holder_coords_(dof_idx, 2) = grid_.bot_left[2] + (k + 0.5) * vertex_spacing[2]; + ++dof_idx; + } } } } } } -template -std::unique_ptr> UniformGridFieldLayout::CreateFieldReal() - const -{ - return std::make_unique>(*this); -} - template int UniformGridFieldLayout::GetNumComponents() const { @@ -78,13 +91,13 @@ int UniformGridFieldLayout::GetNumComponents() const template LO UniformGridFieldLayout::GetNumOwnedDofHolder() const { - return GetNumVertices(); + return GetNumDofHolders(); } template GO UniformGridFieldLayout::GetNumGlobalDofHolder() const { - return GetNumVertices(); + return GetNumDofHolders(); } template @@ -116,7 +129,7 @@ bool UniformGridFieldLayout::IsDistributed() } template -UniformGrid& UniformGridFieldLayout::GetGrid() const +const UniformGrid& UniformGridFieldLayout::GetGrid() const { return grid_; } @@ -137,15 +150,22 @@ LO UniformGridFieldLayout::GetNumVertices() const return num_vertices; } +template +LO UniformGridFieldLayout::GetNumDofHolders() const +{ + return order_ == 0 ? GetNumCells() : GetNumVertices(); +} + 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(); + offsets.fill(0); + LO n = GetNumDofHolders(); + int entity_dim = (order_ == 0) ? static_cast(Dim) : 0; + for (int i = entity_dim + 1; i < ent_offsets_len; ++i) { + offsets[i] = n; + } return offsets; } @@ -156,6 +176,12 @@ ReversePartitionMap2 UniformGridFieldLayout::GetReversePartitionMap( throw std::runtime_error("Unimplemented"); } +template +int UniformGridFieldLayout::GetOrder() const +{ + return order_; +} + // Explicit template instantiations template class UniformGridFieldLayout<2>; template class UniformGridFieldLayout<3>; diff --git a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h index 340d5daf..d01c922e 100644 --- a/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h +++ b/src/pcms/adapter/uniform_grid/uniform_grid_field_layout.h @@ -15,10 +15,8 @@ template class UniformGridFieldLayout : public FieldLayout { public: - UniformGridFieldLayout(UniformGrid& grid, int num_components, - CoordinateSystem coordinate_system); - - std::unique_ptr> CreateFieldReal() const override; + UniformGridFieldLayout(UniformGrid grid, int num_components, + CoordinateSystem coordinate_system, int order = 1); int GetNumComponents() const override; LO GetNumOwnedDofHolder() const override; @@ -35,14 +33,18 @@ class UniformGridFieldLayout : public FieldLayout ReversePartitionMap2 GetReversePartitionMap( const redev::Partition& partition) const override; - UniformGrid& GetGrid() const; + const UniformGrid& GetGrid() const; LO GetNumCells() const; LO GetNumVertices() const; + int GetOrder() const; private: - UniformGrid& grid_; + LO GetNumDofHolders() const; + + UniformGrid grid_; int num_components_; CoordinateSystem coordinate_system_; + int order_; Kokkos::View gids_; Kokkos::View dof_holder_coords_; Kokkos::View owned_; diff --git a/src/pcms/coupler.h b/src/pcms/coupler.h index 494bc0b6..8cd1e6bb 100644 --- a/src/pcms/coupler.h +++ b/src/pcms/coupler.h @@ -2,7 +2,10 @@ #define PCMS_COUPLER_H #include "pcms/utility/common.h" #include "pcms/field_communicator.h" +#include "pcms/configuration.h" +#ifdef PCMS_ENABLE_MESHFIELDS #include "pcms/adapter/meshfields/mesh_fields_adapter.h" +#endif #include "pcms/utility/profile.h" namespace pcms diff --git a/src/pcms/coupler2.cpp b/src/pcms/coupler2.cpp index cc25ff97..d35dfd74 100644 --- a/src/pcms/coupler2.cpp +++ b/src/pcms/coupler2.cpp @@ -17,7 +17,7 @@ FieldLayoutCommunicator& Application2::GetLayoutCommunicator( } const FieldLayout& Application2::AddLayout(std::string name, - std::unique_ptr layout) + std::shared_ptr layout) { layouts_.push_back(std::move(layout)); const FieldLayout& layout_ref = *layouts_.back(); @@ -38,7 +38,7 @@ void Application2::AddField(std::string name, OwnedFieldPtr field, FieldCommunicator2Ptr field_communicator = std::visit( [this, name](auto* field_ptr) -> FieldCommunicator2Ptr { - using T = std::remove_pointer_t::value_type; + using T = typename std::remove_pointer_t::value_type; FieldLayoutCommunicator& layout_communicator = GetLayoutCommunicator(field_ptr->GetLayout()); return std::make_unique>(layout_communicator, diff --git a/src/pcms/coupler2.h b/src/pcms/coupler2.h index 1ee0e763..c176b663 100644 --- a/src/pcms/coupler2.h +++ b/src/pcms/coupler2.h @@ -34,7 +34,7 @@ class Application2 } const FieldLayout& AddLayout(std::string name, - std::unique_ptr layout); + std::shared_ptr layout); // FIXME should take a file path for the parameters, not take adios2 params. // These fields are supposed to be agnostic to adios2... @@ -111,7 +111,7 @@ class Application2 MPI_Comm mpi_comm_; redev::Redev& redev_; redev::Channel channel_; - std::vector> layouts_; + std::vector> layouts_; std::vector fields_; // map is used rather than unordered_map because we give pointers to the // internal data and rehash of unordered_map can cause pointer invalidation. diff --git a/src/pcms/create_field.cpp b/src/pcms/create_field.cpp deleted file mode 100644 index f42a1a0f..00000000 --- a/src/pcms/create_field.cpp +++ /dev/null @@ -1,121 +0,0 @@ -#include "create_field.h" -#include "adapter/meshfields/mesh_fields_adapter2.h" -#include "adapter/meshfields/mesh_fields_adapter_layout.h" -#include "adapter/uniform_grid/uniform_grid_field.h" -#include "adapter/uniform_grid/uniform_grid_field_layout.h" -#include "utility/uniform_grid.h" -#include "localization/point_search.h" - -#include -#include - -namespace pcms -{ - -std::unique_ptr CreateLagrangeLayout( - Omega_h::Mesh& mesh, int order, int num_components, - CoordinateSystem coordinate_system, std::string global_id_name) -{ - - std::array nodes_per_dim; - - switch (order) { - case 1: nodes_per_dim = {1, 0, 0, 0}; break; - case 2: nodes_per_dim = {1, 1, 0, 0}; break; - default: throw std::runtime_error("Unimplemented order"); - } - - return std::make_unique( - mesh, nodes_per_dim, num_components, coordinate_system, global_id_name); -} - -template <> -std::pair>, - std::unique_ptr>> -CreateUniformGridBinaryFieldFromGrid<2>(Omega_h::Mesh& mesh, - UniformGrid<2>& grid) -{ - constexpr unsigned dim = 2; - - // Get total number of vertices - const LO num_vertices = (grid.divisions[0] + 1) * (grid.divisions[1] + 1); - - // Create GridPointSearch for point-in-mesh queries - GridPointSearch2D point_search(mesh, grid.divisions[0], grid.divisions[1]); - - // 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(vertices, vertices_h); - - // Perform point search - 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"); - - // 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).element_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 {std::move(layout), std::move(field)}; -} - -template <> -std::pair>, - std::unique_ptr>> -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::pair>, - std::unique_ptr>> -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 deleted file mode 100644 index 246e03f6..00000000 --- a/src/pcms/create_field.h +++ /dev/null @@ -1,90 +0,0 @@ -#ifndef CREATE_FIELD_H_ -#define CREATE_FIELD_H_ - -#include "adapter/meshfields/mesh_fields_adapter_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" -#include "utility/types.h" -#include "utility/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, std::string global_id_name = "global"); - -/** - * \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 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. - * - * \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::pair containing the layout and field (layout must outlive field) - * - * \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::pair>, - std::unique_ptr>> -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::pair containing the layout and field (layout must outlive field) - */ -template -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 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. - * - * \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::pair containing the layout and field (layout must outlive field) - */ -template -std::pair>, - std::unique_ptr>> -CreateUniformGridBinaryFieldFromGrid(Omega_h::Mesh& mesh, - UniformGrid& grid); - -} // namespace pcms - -#endif // CREATE_FIELD_H_ diff --git a/src/pcms/field_layout.h b/src/pcms/field_layout.h index 7614d1eb..d4da451c 100644 --- a/src/pcms/field_layout.h +++ b/src/pcms/field_layout.h @@ -29,8 +29,6 @@ class FieldT; class FieldLayout { public: - virtual std::unique_ptr> CreateFieldReal() const = 0; - // number of components int virtual GetNumComponents() const = 0; diff --git a/src/pcms/lagrange_field_factory.cpp b/src/pcms/lagrange_field_factory.cpp new file mode 100644 index 00000000..88052732 --- /dev/null +++ b/src/pcms/lagrange_field_factory.cpp @@ -0,0 +1,145 @@ +#include "pcms/lagrange_field_factory.h" +#include "pcms/configuration.h" +#include "pcms/utility/assert.h" +#ifdef PCMS_ENABLE_MESHFIELDS +#include "pcms/adapter/meshfields/mesh_fields_adapter_layout.h" +#include "pcms/adapter/meshfields/mesh_fields_adapter2.h" +#endif +#include "pcms/adapter/omega_h/omega_h_lagrange_layout.h" +#include "pcms/adapter/omega_h/omega_h_lagrange_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field.h" +#include "pcms/utility/uniform_grid.h" + +#include + +namespace pcms +{ + +LagrangeFieldFactory::LagrangeFieldFactory( + std::shared_ptr layout, + std::function>()> create_fn) noexcept + : layout_(std::move(layout)), create_fn_(std::move(create_fn)) +{ +} + +LagrangeFieldFactory LagrangeFieldFactory::FromMesh( + Omega_h::Mesh& mesh, int order, int num_components, + CoordinateSystem coordinate_system, std::string global_id_name, + Backend backend) +{ + if (backend == Backend::MeshFields) { +#ifdef PCMS_ENABLE_MESHFIELDS + std::array nodes_per_dim{}; + switch (order) { + case 1: nodes_per_dim = {1, 0, 0, 0}; break; + case 2: nodes_per_dim = {1, 1, 0, 0}; break; + default: throw pcms_error("Unimplemented Lagrange order"); + } + auto mesh_layout = std::make_shared( + mesh, nodes_per_dim, num_components, coordinate_system, + std::move(global_id_name)); + return LagrangeFieldFactory( + mesh_layout, + [mesh_layout]() { + return std::make_unique>(mesh_layout); + }); +#else + throw pcms_error( + "LagrangeFieldFactory::FromMesh: MeshFields backend requested but " + "PCMS_ENABLE_MESHFIELDS is not set"); +#endif + } + // OmegaH backend + if (order > 1) { + throw pcms_error( + "LagrangeFieldFactory::FromMesh: OmegaH backend only supports " + "Lagrange orders 0 and 1"); + } + auto layout = std::make_shared( + mesh, order, num_components, coordinate_system, std::move(global_id_name)); + return LagrangeFieldFactory( + layout, + [layout]() { + return std::make_unique>(layout); + }); +} + +namespace +{ +template +UniformGrid MakeUniformGrid( + const Rank1View edge_length, + const Rank1View bot_left, + const Rank1View divisions) +{ + PCMS_ALWAYS_ASSERT(edge_length.size() == Dim); + PCMS_ALWAYS_ASSERT(bot_left.size() == Dim); + PCMS_ALWAYS_ASSERT(divisions.size() == Dim); + + UniformGrid grid; + for (unsigned i = 0; i < Dim; ++i) { + grid.edge_length[i] = edge_length[i]; + grid.bot_left[i] = bot_left[i]; + grid.divisions[i] = divisions[i]; + } + return grid; +} +} // namespace + +LagrangeFieldFactory LagrangeFieldFactory::FromUniformGrid( + Rank1View edge_length, + Rank1View bot_left, + Rank1View divisions, + int num_components, + CoordinateSystem coordinate_system, + int order) +{ + if (order != 0 && order != 1) { + throw std::invalid_argument( + "LagrangeFieldFactory::FromUniformGrid: only orders 0 and 1 are supported"); + } + if (edge_length.size() != bot_left.size() || + edge_length.size() != divisions.size()) { + throw std::invalid_argument( + "edge_length, bot_left, and divisions must have the same size"); + } + switch (edge_length.size()) { + case 2: { + auto ug_layout = std::make_shared>( + MakeUniformGrid<2>(edge_length, bot_left, divisions), num_components, + coordinate_system, order); + return LagrangeFieldFactory( + ug_layout, + [ug_layout]() { + return std::make_unique>(ug_layout); + }); + } + case 3: { + auto ug_layout = std::make_shared>( + MakeUniformGrid<3>(edge_length, bot_left, divisions), num_components, + coordinate_system, order); + return LagrangeFieldFactory( + ug_layout, + [ug_layout]() { + return std::make_unique>(ug_layout); + }); + } + default: + throw std::invalid_argument( + "LagrangeFieldFactory::FromUniformGrid: only dim 2 and 3 are supported"); + } +} + +std::shared_ptr +LagrangeFieldFactory::GetLayout() const noexcept +{ + return layout_; +} + +std::unique_ptr> LagrangeFieldFactory::CreateFieldReal() const +{ + return create_fn_(); +} + +} // namespace pcms diff --git a/src/pcms/lagrange_field_factory.h b/src/pcms/lagrange_field_factory.h new file mode 100644 index 00000000..61a810aa --- /dev/null +++ b/src/pcms/lagrange_field_factory.h @@ -0,0 +1,61 @@ +#ifndef PCMS_LAGRANGE_FIELD_FACTORY_H +#define PCMS_LAGRANGE_FIELD_FACTORY_H + +#include +#include "pcms/configuration.h" +#include "pcms/field.h" +#include "pcms/field_layout.h" +#include "pcms/coordinate_system.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" +#include "pcms/utility/memory_spaces.h" + +#include +#include +#include + +namespace pcms +{ + +class LagrangeFieldFactory +{ +public: + enum class Backend { MeshFields, OmegaH }; + +#ifdef PCMS_ENABLE_MESHFIELDS + static constexpr Backend DefaultBackend = Backend::MeshFields; +#else + static constexpr Backend DefaultBackend = Backend::OmegaH; +#endif + + // Unstructured mesh — dispatches to MeshFields or native Omega_h backend + [[nodiscard]] static LagrangeFieldFactory FromMesh( + Omega_h::Mesh& mesh, int order, int num_components, + CoordinateSystem coordinate_system, + std::string global_id_name = "global", + Backend backend = DefaultBackend); + + // Structured uniform grid — order-1 H1-conforming nodal field on a regular grid + [[nodiscard]] static LagrangeFieldFactory FromUniformGrid( + Rank1View edge_length, + Rank1View bot_left, + Rank1View divisions, + int num_components, + CoordinateSystem coordinate_system, + int order = 1); + + [[nodiscard]] std::shared_ptr GetLayout() const noexcept; + [[nodiscard]] std::unique_ptr> CreateFieldReal() const; + +private: + explicit LagrangeFieldFactory( + std::shared_ptr layout, + std::function>()> create_fn) noexcept; + + std::shared_ptr layout_; + std::function>()> create_fn_; +}; + +} // namespace pcms + +#endif // PCMS_LAGRANGE_FIELD_FACTORY_H diff --git a/src/pcms/localization/CMakeLists.txt b/src/pcms/localization/CMakeLists.txt index c46aaad4..ecfc89e5 100644 --- a/src/pcms/localization/CMakeLists.txt +++ b/src/pcms/localization/CMakeLists.txt @@ -18,7 +18,7 @@ target_include_directories( "$" "$") target_link_libraries(pcms_localization PUBLIC Kokkos::kokkos pcms::utility Omega_h::omega_h) -target_compile_features(pcms_localization PUBLIC cxx_std_17) +target_compile_features(pcms_localization PUBLIC cxx_std_20) set_target_properties( pcms_localization PROPERTIES OUTPUT_NAME pcmslocalization EXPORT_NAME @@ -47,4 +47,4 @@ install( EXPORT pcms_localization-targets NAMESPACE pcms:: DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms -) \ No newline at end of file +) diff --git a/src/pcms/nodal_field_factory.cpp b/src/pcms/nodal_field_factory.cpp new file mode 100644 index 00000000..1878cb5a --- /dev/null +++ b/src/pcms/nodal_field_factory.cpp @@ -0,0 +1,44 @@ +#include "pcms/nodal_field_factory.h" +#include "pcms/adapter/point_cloud/point_cloud_layout.h" +#include "pcms/adapter/point_cloud/point_cloud.h" + +#include + +namespace pcms +{ + +NodalFieldFactory::NodalFieldFactory( + std::shared_ptr layout, + std::function>()> create_fn) noexcept + : layout_(std::move(layout)), create_fn_(std::move(create_fn)) +{ +} + +NodalFieldFactory NodalFieldFactory::Create( + Rank2View coords, + CoordinateSystem coordinate_system) +{ + int dim = static_cast(coords.extent(1)); + Kokkos::View host_view( + coords.data_handle(), coords.extent(0), coords.extent(1)); + auto device_view = + Kokkos::create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace{}, host_view); + auto pc_layout = + std::make_shared(dim, device_view, coordinate_system); + return NodalFieldFactory( + pc_layout, + [pc_layout]() { return std::make_unique(pc_layout); }); +} + +std::shared_ptr +NodalFieldFactory::GetLayout() const noexcept +{ + return layout_; +} + +std::unique_ptr> NodalFieldFactory::CreateFieldReal() const +{ + return create_fn_(); +} + +} // namespace pcms diff --git a/src/pcms/nodal_field_factory.h b/src/pcms/nodal_field_factory.h new file mode 100644 index 00000000..23cbf8fb --- /dev/null +++ b/src/pcms/nodal_field_factory.h @@ -0,0 +1,37 @@ +#ifndef PCMS_NODAL_FIELD_FACTORY_H +#define PCMS_NODAL_FIELD_FACTORY_H + +#include "pcms/field.h" +#include "pcms/field_layout.h" +#include "pcms/coordinate_system.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/memory_spaces.h" + +#include +#include + +namespace pcms +{ + +class NodalFieldFactory +{ +public: + [[nodiscard]] static NodalFieldFactory Create( + Rank2View coords, + CoordinateSystem coordinate_system); + + [[nodiscard]] std::shared_ptr GetLayout() const noexcept; + [[nodiscard]] std::unique_ptr> CreateFieldReal() const; + +private: + explicit NodalFieldFactory( + std::shared_ptr layout, + std::function>()> create_fn) noexcept; + + std::shared_ptr layout_; + std::function>()> create_fn_; +}; + +} // namespace pcms + +#endif // PCMS_NODAL_FIELD_FACTORY_H diff --git a/src/pcms/pythonapi/CMakeLists.txt b/src/pcms/pythonapi/CMakeLists.txt index 999ee1b5..7094123b 100644 --- a/src/pcms/pythonapi/CMakeLists.txt +++ b/src/pcms/pythonapi/CMakeLists.txt @@ -15,7 +15,7 @@ pybind11_add_module(pcms bind_mls_interpolation.cpp bind_mesh_utilities.cpp ) -set_target_properties(pcms PROPERTIES CXX_STANDARD 17) +target_compile_features(pcms PUBLIC cxx_std_20) target_link_libraries(pcms PRIVATE Kokkos::kokkos pybind11::module MPI::MPI_C pcms::core pcms::transfer) # if we are building as a python package from the pyproject.toml diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp index 1643e243..b13f7629 100644 --- a/src/pcms/pythonapi/bind_field_base.cpp +++ b/src/pcms/pythonapi/bind_field_base.cpp @@ -3,7 +3,10 @@ #include #include "pcms/coordinate_system.h" #include "pcms/coordinate.h" -#include "pcms/create_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_field_layout.h" +#include "pcms/adapter/uniform_grid/uniform_grid_binary_field.h" +#include "pcms/lagrange_field_factory.h" #include "pcms/utility/uniform_grid.h" #include "numpy_array_transform.h" @@ -131,17 +134,64 @@ void bind_coordinate_module(py::module& m) 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") = 1, - py::arg("coordinate_system") = CoordinateSystem::Cartesian); + // Bind LagrangeFieldFactory + py::class_(m, "LagrangeFieldFactory") + .def_static( + "from_mesh", + [](Omega_h::Mesh& mesh, int order, int num_components, + CoordinateSystem coordinate_system) { + return LagrangeFieldFactory::FromMesh(mesh, order, num_components, + coordinate_system); + }, + py::arg("mesh"), py::arg("order"), py::arg("num_components") = 1, + py::arg("coordinate_system") = CoordinateSystem::Cartesian, + "Create a LagrangeFieldFactory from an Omega_h mesh") + + .def_static( + "from_uniform_grid", + [](const UniformGrid<2>& grid, int num_components, CoordinateSystem cs, + int order) { + auto el = grid.edge_length; + auto bl = grid.bot_left; + auto div = grid.divisions; + Rank1View el_view(el.data(), 2); + Rank1View bl_view(bl.data(), 2); + Rank1View div_view(div.data(), 2); + return LagrangeFieldFactory::FromUniformGrid(el_view, bl_view, div_view, + num_components, cs, order); + }, + py::arg("grid"), py::arg("num_components") = 1, + py::arg("coordinate_system") = CoordinateSystem::Cartesian, + py::arg("order") = 1, + "Create a LagrangeFieldFactory from a 2D uniform grid") + + .def_static( + "from_uniform_grid", + [](const UniformGrid<3>& grid, int num_components, CoordinateSystem cs, + int order) { + auto el = grid.edge_length; + auto bl = grid.bot_left; + auto div = grid.divisions; + Rank1View el_view(el.data(), 3); + Rank1View bl_view(bl.data(), 3); + Rank1View div_view(div.data(), 3); + return LagrangeFieldFactory::FromUniformGrid(el_view, bl_view, div_view, + num_components, cs, order); + }, + py::arg("grid"), py::arg("num_components") = 1, + py::arg("coordinate_system") = CoordinateSystem::Cartesian, + py::arg("order") = 1, + "Create a LagrangeFieldFactory from a 3D uniform grid") + + .def("get_layout", &LagrangeFieldFactory::GetLayout, + "Get the field layout (shared_ptr)") + + .def( + "create_field_real", + [](const LagrangeFieldFactory& self) { + return std::shared_ptr>(self.CreateFieldReal()); + }, + "Create a real-valued field"); // Bind CreateUniformGridFromMesh for 2D m.def( @@ -161,21 +211,17 @@ void bind_create_field_module(py::module& m) 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))); + auto [layout, field] = + CreateUniformGridBinaryField<2>(mesh, divisions); + return py::make_tuple(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."); + "Create a 2D vertex mask field indicating inside/outside mesh"); + } template diff --git a/src/pcms/pythonapi/bind_field_layout.cpp b/src/pcms/pythonapi/bind_field_layout.cpp index a2277c12..9a6cfa9b 100644 --- a/src/pcms/pythonapi/bind_field_layout.cpp +++ b/src/pcms/pythonapi/bind_field_layout.cpp @@ -13,9 +13,6 @@ 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::CreateFieldReal, - "Create a field with this layout") - .def("get_num_components", &FieldLayout::GetNumComponents, "Get the number of components in the field") diff --git a/src/pcms/pythonapi/bind_omega_h_field2.cpp b/src/pcms/pythonapi/bind_omega_h_field2.cpp index 1a6a5428..b966affb 100644 --- a/src/pcms/pythonapi/bind_omega_h_field2.cpp +++ b/src/pcms/pythonapi/bind_omega_h_field2.cpp @@ -28,7 +28,7 @@ void bind_omega_h_field2(py::module& m) // Bind MeshFieldsAdapter2 class template instantiation for Real type py::class_, FieldT, std::shared_ptr>>(m, "MeshFieldsAdapter2") - .def(py::init(), py::arg("layout"), + .def(py::init>(), py::arg("layout"), "Constructor for MeshFieldsAdapter2") .def( diff --git a/src/pcms/pythonapi/bind_omega_h_field_layout.cpp b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp index b50fb7e8..ebda5dd2 100644 --- a/src/pcms/pythonapi/bind_omega_h_field_layout.cpp +++ b/src/pcms/pythonapi/bind_omega_h_field_layout.cpp @@ -23,13 +23,6 @@ void bind_omega_h_field_layout_module(py::module& m) py::arg("coordinate_system"), py::arg("global_id_name") = "global", "Constructor for MeshFieldsAdapterLayout") - .def( - "create_field", - [](MeshFieldsAdapterLayout& self) { - return std::shared_ptr>(self.CreateFieldReal()); - }, - "Create a field with this layout") - .def("get_num_components", &MeshFieldsAdapterLayout::GetNumComponents, "Get the number of components in the field") diff --git a/src/pcms/pythonapi/bind_uniform_grid_field.cpp b/src/pcms/pythonapi/bind_uniform_grid_field.cpp index c6940985..7b97b3ae 100644 --- a/src/pcms/pythonapi/bind_uniform_grid_field.cpp +++ b/src/pcms/pythonapi/bind_uniform_grid_field.cpp @@ -18,7 +18,7 @@ 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"), + .def(py::init>>(), py::arg("layout"), "Constructor for UniformGridField2D") .def( @@ -128,7 +128,7 @@ void bind_uniform_grid_field_module(py::module& m) // Bind UniformGridField class for 3D py::class_, FieldT, std::shared_ptr>>(m, "UniformGridField3D") - .def(py::init&>(), py::arg("layout"), + .def(py::init>>(), py::arg("layout"), "Constructor for UniformGridField3D") .def( diff --git a/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp index d913f3d1..cf464eef 100644 --- a/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp +++ b/src/pcms/pythonapi/bind_uniform_grid_field_layout.cpp @@ -72,17 +72,11 @@ void bind_uniform_grid_field_layout_module(py::module& m) py::class_, FieldLayout, std::shared_ptr>>( m, "UniformGridFieldLayout2D") - .def(py::init&, int, CoordinateSystem>(), py::arg("grid"), + .def(py::init, int, CoordinateSystem, int>(), py::arg("grid"), py::arg("num_components"), py::arg("coordinate_system"), + py::arg("order") = 1, "Constructor for UniformGridFieldLayout2D") - .def( - "create_field", - [](UniformGridFieldLayout<2>& self) { - return std::shared_ptr>(self.CreateFieldReal()); - }, - "Create a field with this layout") - .def("get_num_components", &UniformGridFieldLayout<2>::GetNumComponents, "Get the number of components in the field") @@ -151,23 +145,20 @@ void bind_uniform_grid_field_layout_module(py::module& m) "Get the number of cells in the grid") .def("get_num_vertices", &UniformGridFieldLayout<2>::GetNumVertices, - "Get the number of vertices in the grid"); + "Get the number of vertices in the grid") + + .def("get_order", &UniformGridFieldLayout<2>::GetOrder, + "Get the Lagrange order for the layout"); // Bind the UniformGridFieldLayout class for 3D py::class_, FieldLayout, std::shared_ptr>>( m, "UniformGridFieldLayout3D") - .def(py::init&, int, CoordinateSystem>(), py::arg("grid"), + .def(py::init, int, CoordinateSystem, int>(), py::arg("grid"), py::arg("num_components"), py::arg("coordinate_system"), + py::arg("order") = 1, "Constructor for UniformGridFieldLayout3D") - .def( - "create_field", - [](UniformGridFieldLayout<3>& self) { - return std::shared_ptr>(self.CreateFieldReal()); - }, - "Create a field with this layout") - .def("get_num_components", &UniformGridFieldLayout<3>::GetNumComponents, "Get the number of components in the field") @@ -236,7 +227,10 @@ void bind_uniform_grid_field_layout_module(py::module& m) "Get the number of cells in the grid") .def("get_num_vertices", &UniformGridFieldLayout<3>::GetNumVertices, - "Get the number of vertices in the grid"); + "Get the number of vertices in the grid") + + .def("get_order", &UniformGridFieldLayout<3>::GetOrder, + "Get the Lagrange order for the layout"); } } // namespace pcms diff --git a/src/pcms/pythonapi/test_field_copy.py b/src/pcms/pythonapi/test_field_copy.py index 3d569e9b..537341e3 100644 --- a/src/pcms/pythonapi/test_field_copy.py +++ b/src/pcms/pythonapi/test_field_copy.py @@ -20,15 +20,13 @@ def test_copy(world, dim, order, num_components): print(f" Mesh type: {type(mesh)}") print(f" Mesh object: {mesh}") - # Create layout - print(" About to create layout...") - layout = pcms.create_lagrange_layout( - mesh, - order, - num_components, - pcms.CoordinateSystem.Cartesian + # Create factory and layout + print(" About to create factory...") + factory = pcms.LagrangeFieldFactory.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian ) - print(f" Layout created successfully") + layout = factory.get_layout() + print(f" Factory and layout created successfully") print(f"Testing dim={dim}, order={order}, num_components={num_components}...") # Get number of data points @@ -40,13 +38,13 @@ def test_copy(world, dim, order, num_components): print(f" Created array of IDs from 0 to {ndata-1}") # Create original field and set data - original = layout.create_field() + original = factory.create_field_real() 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() + copied = factory.create_field_real() pcms.copy_field(original, copied) print(" Copied data to new field") diff --git a/src/pcms/pythonapi/test_omega_h_field.py b/src/pcms/pythonapi/test_omega_h_field.py index 8e741b62..95d6a0b1 100644 --- a/src/pcms/pythonapi/test_omega_h_field.py +++ b/src/pcms/pythonapi/test_omega_h_field.py @@ -23,13 +23,11 @@ def test_layout_methods(world, dim, order, num_components): False ) - # Create layout - layout = pcms.create_lagrange_layout( - mesh, - order, - num_components, - pcms.CoordinateSystem.Cartesian + # Create factory and layout + factory = pcms.LagrangeFieldFactory.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian ) + layout = factory.get_layout() # Test layout methods num_comp = layout.get_num_components() @@ -86,8 +84,8 @@ def test_layout_methods(world, dim, order, num_components): assert owned_size == num_owned * num_components assert global_size == num_global * num_components - # Create a field from the layout - field = layout.create_field() + # Create a field from the factory + field = factory.create_field_real() print(f" Created field: {type(field)}") # Test setting and getting data @@ -124,16 +122,14 @@ def test_field_evaluation(world, dim, order, num_components): False ) - # Create layout - layout = pcms.create_lagrange_layout( - mesh, - order, - num_components, - pcms.CoordinateSystem.Cartesian + # Create factory and field + factory = pcms.LagrangeFieldFactory.from_mesh( + mesh, order, num_components, pcms.CoordinateSystem.Cartesian ) + layout = factory.get_layout() # Create field - field = layout.create_field() + field = factory.create_field_real() # 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...") diff --git a/src/pcms/pythonapi/test_uniform_grid_field.py b/src/pcms/pythonapi/test_uniform_grid_field.py index b89864c6..7bb291fe 100644 --- a/src/pcms/pythonapi/test_uniform_grid_field.py +++ b/src/pcms/pythonapi/test_uniform_grid_field.py @@ -28,9 +28,12 @@ def test_uniform_grid_field_creation(): print(f"Num vertices: {layout.get_num_vertices()}") # 5x5 = 25 vertices # Create field - field = layout.create_field() + factory = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + field = factory.create_field_real() print(f"Field created: {field is not None}") - + return grid, layout, field @@ -45,8 +48,10 @@ def test_uniform_grid_field_data_operations(): layout = pcms.UniformGridFieldLayout2D( grid, 1, pcms.CoordinateSystem.Cartesian ) - ug_field = layout.create_field() - + ug_field = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ).create_field_real() + # Get the number of DOF holders (vertices) num_vertices = layout.get_num_vertices() print(f"Number of vertices: {num_vertices}") # Should be 9 @@ -77,7 +82,9 @@ def test_uniform_grid_field_mdspan_2d(): layout = pcms.UniformGridFieldLayout2D( grid, 1, pcms.CoordinateSystem.Cartesian ) - ug_field = layout.create_field() + ug_field = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ).create_field_real() num_vertices = layout.get_num_vertices() data = np.arange(num_vertices, dtype=np.float64) @@ -103,8 +110,10 @@ def test_uniform_grid_field_evaluation(): layout = pcms.UniformGridFieldLayout2D( grid, 1, pcms.CoordinateSystem.Cartesian ) - ug_field = layout.create_field() - + ug_field = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ).create_field_real() + # Set simple linear field data: f(x,y) = x + y num_vertices = layout.get_num_vertices() coords = layout.get_dof_holder_coordinates() @@ -164,9 +173,11 @@ def test_3d_uniform_grid(): ) print(f"3D Grid vertices: {layout.get_num_vertices()}") # 3x3x3 = 27 - - ug_field = layout.create_field() - + + ug_field = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ).create_field_real() + # Set and get data num_vertices = layout.get_num_vertices() data = np.ones(num_vertices) * 42.0 @@ -187,7 +198,9 @@ def test_uniform_grid_field_mdspan_3d(): layout = pcms.UniformGridFieldLayout3D( grid, 1, pcms.CoordinateSystem.Cartesian ) - ug_field = layout.create_field() + ug_field = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ).create_field_real() num_vertices = layout.get_num_vertices() data = np.arange(num_vertices, dtype=np.float64) @@ -234,11 +247,10 @@ def test_uniform_grid_workflow(world): print(f"Created mask field with {len(mask_data)} vertices") # Create Omega_h field layout with linear elements - omega_h_layout = pcms.create_lagrange_layout( - mesh, 1 - ) - omega_h_field = omega_h_layout.create_field() - + omega_h_factory = pcms.LagrangeFieldFactory.from_mesh(mesh, 1) + omega_h_layout = omega_h_factory.get_layout() + omega_h_field = omega_h_factory.create_field_real() + # 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() @@ -256,7 +268,10 @@ def test_uniform_grid_workflow(world): ug_layout = pcms.UniformGridFieldLayout2D( grid, 1, pcms.CoordinateSystem.Cartesian ) - ug_field = ug_layout.create_field() + ug_factory = pcms.LagrangeFieldFactory.from_uniform_grid( + grid, 1, pcms.CoordinateSystem.Cartesian + ) + ug_field = ug_factory.create_field_real() # Transfer from omega_h field to uniform grid field using interpolation pcms.interpolate_field(omega_h_field, ug_field) @@ -340,7 +355,7 @@ def test_uniform_grid_workflow(world): # Test deserialization print("\nTesting deserialization:") - ug_field_copy = ug_layout.create_field() + ug_field_copy = ug_factory.create_field_real() ug_field_copy.deserialize(buffer, permutation) ug_field_copy_data = ug_field_copy.get_dof_holder_data() @@ -385,14 +400,12 @@ def test_omega_h_to_omega_h_transfer_workflow(world): ) # Create Lagrange layouts and fields - src_layout = pcms.create_lagrange_layout( - src_mesh, 1 - ) - tgt_layout = pcms.create_lagrange_layout( - tgt_mesh, 1 - ) - src_field = src_layout.create_field() - tgt_field = tgt_layout.create_field() + src_factory = pcms.LagrangeFieldFactory.from_mesh(src_mesh, 1) + tgt_factory = pcms.LagrangeFieldFactory.from_mesh(tgt_mesh, 1) + src_layout = src_factory.get_layout() + tgt_layout = tgt_factory.get_layout() + src_field = src_factory.create_field_real() + tgt_field = tgt_factory.create_field_real() # Initialize source field with f(x,y) = x + 2*y src_coords = src_layout.get_dof_holder_coordinates() diff --git a/src/pcms/transfer/CMakeLists.txt b/src/pcms/transfer/CMakeLists.txt index cbd8889a..e33aa5e4 100644 --- a/src/pcms/transfer/CMakeLists.txt +++ b/src/pcms/transfer/CMakeLists.txt @@ -26,7 +26,7 @@ set(PCMS_FIELD_TRANSFER_SOURCES mls_interpolation.cpp interpolation_base.cpp) -if(PCMS_ENABLE_PETSC) +if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) list(APPEND PCMS_FIELD_TRANSFER_HEADERS calculate_load_vector.hpp calculate_mass_matrix.hpp @@ -51,16 +51,18 @@ target_sources(pcms_transfer PUBLIC BASE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/.. FILES ${PCMS_FIELD_TRANSFER_HEADERS}) add_library(pcms::transfer ALIAS pcms_transfer) -target_compile_features(pcms_transfer PUBLIC cxx_std_17) +target_compile_features(pcms_transfer PUBLIC cxx_std_20) target_link_libraries(pcms_transfer PUBLIC pcms::core Omega_h::omega_h PRIVATE Kokkos::kokkoskernels) -target_link_libraries(pcms_transfer PUBLIC meshfields::meshfields) +if(PCMS_ENABLE_MESHFIELDS) + target_link_libraries(pcms_transfer PUBLIC meshfields::meshfields) +endif() -if(PCMS_ENABLE_PETSC) +if(PCMS_ENABLE_PETSC AND PCMS_ENABLE_MESHFIELDS) target_link_libraries(pcms_transfer PRIVATE PETSc::PETSc) endif() diff --git a/src/pcms/utility/CMakeLists.txt b/src/pcms/utility/CMakeLists.txt index 8dd53372..bb4f6b42 100644 --- a/src/pcms/utility/CMakeLists.txt +++ b/src/pcms/utility/CMakeLists.txt @@ -30,7 +30,7 @@ target_include_directories( "$" "$") target_link_libraries(pcms_utility PUBLIC Kokkos::kokkos perfstubs) -target_compile_features(pcms_utility PUBLIC cxx_std_17) +target_compile_features(pcms_utility PUBLIC cxx_std_20) set_target_properties( pcms_utility PROPERTIES OUTPUT_NAME pcmsutility EXPORT_NAME diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c27a6de2..d623d49e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -398,7 +398,8 @@ if(Catch2_FOUND) test_uniform_grid_field.cpp test_interpolation_class.cpp test_mesh_geometry.cpp - test_omega_h_field2_outofbounds.cpp) + test_omega_h_field2_outofbounds.cpp + test_omega_h_lagrange_field.cpp) endif() if(PCMS_ENABLE_MESHFIELDS) diff --git a/test/field_test_utils.h b/test/field_test_utils.h new file mode 100644 index 00000000..d961e286 --- /dev/null +++ b/test/field_test_utils.h @@ -0,0 +1,187 @@ +#ifndef PCMS_TEST_FIELD_TEST_UTILS_H +#define PCMS_TEST_FIELD_TEST_UTILS_H + +#include +#include +#include "pcms/field.h" +#include "pcms/coordinate_system.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/memory_spaces.h" +#include +#include + +// Shared utilities for field tests that apply equally to MeshFields-backed +// and native OmegaH-backed field implementations. + +namespace pcms::test +{ + +// Affine test function — exactly representable on linear elements. +KOKKOS_INLINE_FUNCTION Real linear_f(Real x, Real y) +{ + return x + 2.0 * y; +} + +// Interior test points for a unit [0,1]^2 box mesh. +inline std::vector StandardEvalCoords2D() +{ + return {0.1, 0.2, 0.5, 0.5, 0.7, 0.3, 0.9, 0.1, 0.2, 0.8}; +} + +template +inline std::vector EvaluateReferenceFunction( + const std::vector& pts, Func func) +{ + using MemorySpace = typename ExecutionSpace::memory_space; + + int n = static_cast(pts.size()) / 2; + Kokkos::View pts_device("pts_device", n); + auto pts_host = + Kokkos::View(pts.data(), n, 2); + deep_copy_mismatch_layouts(pts_device, pts_host); + + Kokkos::View expected_device("expected_device", n); + Kokkos::parallel_for( + "field_test_utils_expected", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(int i) { + expected_device(i) = func(pts_device(i, 0), pts_device(i, 1)); + }); + + auto expected_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), expected_device); + return std::vector(expected_host.data(), + expected_host.data() + expected_host.extent(0)); +} + +// Set scalar DOF data by sampling func at each DOF-holder coordinate. +template +inline void SetField(FieldT& field, Func func) +{ + using MemorySpace = typename ExecutionSpace::memory_space; + + auto dof_coords = field.GetLayout().GetDOFHolderCoordinates().GetCoordinates(); + int n = static_cast(dof_coords.extent(0)); + Kokkos::View coords_device("coords_device", n); + auto coords_host = Kokkos::View( + dof_coords.data_handle(), n, 2); + deep_copy_mismatch_layouts(coords_device, coords_host); + + Kokkos::View data_device("data_device", n); + Kokkos::parallel_for( + "field_test_utils_set_field", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(int i) { + data_device(i) = func(coords_device(i, 0), coords_device(i, 1)); + }); + + auto data_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), data_device); + field.SetDOFHolderData( + Rank1View(data_host.data(), n)); +} + +// Evaluate field at points known to be outside the mesh and verify they all +// get the given fill_value. +inline void CheckFillMode(FieldT& field, Real fill_value, + const std::vector& outside_pts) +{ + int n = static_cast(outside_pts.size()) / 2; + std::vector eval(n); + + Rank2View coords_view(outside_pts.data(), n, 2); + Rank1View eval_view(eval.data(), n); + CoordinateView cv{field.GetCoordinateSystem(), coords_view}; + FieldDataView dv{eval_view, field.GetCoordinateSystem()}; + + auto hint = field.GetLocalizationHint(cv); + field.Evaluate(hint, dv); + + for (int i = 0; i < n; ++i) { + REQUIRE(eval[i] == fill_value); + } +} + +// Check that serialize followed by deserialize round-trips the data. +// Uses an identity permutation so permutation[i] = i. +inline void CheckSerializeDeserialize(FieldT& field) +{ + auto data_before = field.GetDOFHolderData(); + int n = static_cast(data_before.size()); + + std::vector buffer(n); + std::vector perm(n); + for (int i = 0; i < n; ++i) + perm[i] = i; + + Rank1View buf_view(buffer.data(), n); + Rank1View perm_view(perm.data(), n); + + field.Serialize(buf_view, perm_view); + field.Deserialize(Rank1View(buf_view), perm_view); + + auto data_after = field.GetDOFHolderData(); + REQUIRE(data_after.size() == data_before.size()); + for (int i = 0; i < n; ++i) { + REQUIRE(data_after[i] == Catch::Approx(data_before[i])); + } +} + +// Evaluate field at explicit test points and check each result against func. +template +void CheckEvaluation(FieldT& field, const std::vector& pts, + Func func, + double abs_tol = 1e-10) +{ + int n = static_cast(pts.size()) / 2; + std::vector eval(n); + + Rank2View coords_view(pts.data(), n, 2); + Rank1View eval_view(eval.data(), n); + CoordinateView cv{field.GetCoordinateSystem(), coords_view}; + FieldDataView dv{eval_view, field.GetCoordinateSystem()}; + + auto hint = field.GetLocalizationHint(cv); + field.Evaluate(hint, dv); + + auto expected = EvaluateReferenceFunction(pts, func); + + for (int i = 0; i < n; ++i) { + INFO("Point " << i << " (" << pts[2*i] << ", " << pts[2*i+1] << ")" + << " got=" << eval[i] << " expected=" << expected[i]); + REQUIRE(eval[i] == Catch::Approx(expected[i]).margin(abs_tol)); + } +} + +template +void CheckEvaluationWithFill(FieldT& field, const std::vector& pts, + const std::vector& is_inside, Func func, + Real fill_value, double abs_tol = 1e-10) +{ + int n = static_cast(pts.size()) / 2; + REQUIRE(static_cast(is_inside.size()) == n); + + std::vector eval(n); + + Rank2View coords_view(pts.data(), n, 2); + Rank1View eval_view(eval.data(), n); + CoordinateView cv{field.GetCoordinateSystem(), coords_view}; + FieldDataView dv{eval_view, field.GetCoordinateSystem()}; + + auto hint = field.GetLocalizationHint(cv); + field.Evaluate(hint, dv); + + auto expected = EvaluateReferenceFunction(pts, func); + + for (int i = 0; i < n; ++i) { + INFO("Point " << i << " (" << pts[2*i] << ", " << pts[2*i+1] << ")" + << " got=" << eval[i]); + if (is_inside[i]) { + INFO(" expected=" << expected[i]); + REQUIRE(eval[i] == Catch::Approx(expected[i]).margin(abs_tol)); + } else { + REQUIRE(eval[i] == fill_value); + } + } +} + +} // namespace pcms::test + +#endif // PCMS_TEST_FIELD_TEST_UTILS_H diff --git a/test/test_field_communication.cpp b/test/test_field_communication.cpp index 53eb4146..2c52d420 100644 --- a/test/test_field_communication.cpp +++ b/test/test_field_communication.cpp @@ -10,7 +10,7 @@ #include "pcms/adapter/meshfields/mesh_fields_adapter2.h" #include "pcms/field_communicator2.h" #include "pcms/field_communicator.h" -#include "pcms/create_field.h" +#include "pcms/lagrange_field_factory.h" #include "test_support.h" namespace ts = test_support; @@ -66,8 +66,9 @@ void client1(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, auto channel = rdv.CreateAdiosChannel("field2_chan1", params, redev::TransportType::BP4); - auto layout = pcms::CreateLagrangeLayout(mesh, order, 1, - pcms::CoordinateSystem::Cartesian); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, order, 1, pcms::CoordinateSystem::Cartesian); + auto layout = factory.GetLayout(); auto gids = layout->GetGids(); const auto n = layout->GetNumOwnedDofHolder(); Omega_h::HostWrite ids(n); @@ -76,7 +77,7 @@ void client1(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, "id gid", Kokkos::RangePolicy(0, n), [=](int i) { ids[i] = gids[i]; }); - auto field = layout->CreateFieldReal(); + auto field = factory.CreateFieldReal(); field->SetDOFHolderData(pcms::make_const_array_view(ids)); pcms::FieldLayoutCommunicator layout_comm(comm_name + "1", comm, rdv, channel, @@ -95,12 +96,13 @@ void client2(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, auto channel = rdv.CreateAdiosChannel("field2_chan2", params, redev::TransportType::BP4); - auto layout = pcms::CreateLagrangeLayout(mesh, order, 1, - pcms::CoordinateSystem::Cartesian); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, order, 1, pcms::CoordinateSystem::Cartesian); + auto layout = factory.GetLayout(); auto gids = layout->GetGids(); const auto n = layout->GetNumOwnedDofHolder(); - auto field = layout->CreateFieldReal(); + auto field = factory.CreateFieldReal(); pcms::FieldLayoutCommunicator layout_comm(comm_name + "2", comm, rdv, channel, *layout); pcms::FieldCommunicator2 field_comm(layout_comm, *field); @@ -153,15 +155,16 @@ void server(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, auto channel2 = rdv.CreateAdiosChannel("field2_chan2", params, redev::TransportType::BP4); - auto layout = pcms::CreateLagrangeLayout(mesh, order, 1, - pcms::CoordinateSystem::Cartesian); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, order, 1, pcms::CoordinateSystem::Cartesian); + auto layout = factory.GetLayout(); const auto n = layout->GetNumOwnedDofHolder(); Omega_h::HostWrite ids(n); Kokkos::parallel_for( "id 0", Kokkos::RangePolicy(0, n), [=](int i) { ids[i] = 0; }); - auto field = layout->CreateFieldReal(); + auto field = factory.CreateFieldReal(); pcms::FieldLayoutCommunicator layout_comm1(comm_name + "1", comm, rdv, channel1, *layout); pcms::FieldLayoutCommunicator layout_comm2(comm_name + "2", comm, rdv, diff --git a/test/test_field_copy.cpp b/test/test_field_copy.cpp index 015a8167..445c15ff 100644 --- a/test/test_field_copy.cpp +++ b/test/test_field_copy.cpp @@ -6,7 +6,8 @@ #include #include "pcms/adapter/meshfields/mesh_fields_adapter.h" #include "pcms/adapter/meshfields/mesh_fields_adapter2.h" -#include "pcms/create_field.h" +#include "pcms/lagrange_field_factory.h" +#include "pcms/utility/assert.h" #include using pcms::Real; @@ -19,18 +20,22 @@ void test_copy(Omega_h::CommPtr world, int dim, int order, int num_components) auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, nx, ny, nz, false); - auto layout = pcms::CreateLagrangeLayout(mesh, order, num_components, - pcms::CoordinateSystem::Cartesian); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, order, num_components, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::OmegaH); + // FIXME this indicates that we are not exposing the right API from the fields to be able + // to allocate buffers for getting/setting the field data. + auto layout = factory.GetLayout(); int ndata = layout->GetNumOwnedDofHolder() * num_components; Omega_h::HostWrite ids(ndata); Kokkos::parallel_for( Kokkos::RangePolicy(0, ndata), [=](int i) { ids[i] = i; }); - auto original = layout->CreateFieldReal(); + auto original = factory.CreateFieldReal(); original->SetDOFHolderData(pcms::make_const_array_view(ids)); - auto copied = layout->CreateFieldReal(); + auto copied = factory.CreateFieldReal(); pcms::copy_field2(*original, *copied); auto copied_array = copied->GetDOFHolderData(); @@ -48,6 +53,14 @@ void test_copy(Omega_h::CommPtr world, int dim, int order, int num_components) TEST_CASE("copy omega_h_field2 data") { auto lib = Omega_h::Library{}; + test_copy(lib.world(), 2, 0, 1); test_copy(lib.world(), 2, 1, 1); - test_copy(lib.world(), 2, 2, 1); + auto mesh = + Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 1, 100, 100, 0, false); + REQUIRE_THROWS_AS( + pcms::LagrangeFieldFactory::FromMesh( + mesh, + 2, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::OmegaH), + pcms::pcms_error); } diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp index 7a25defe..7b00702f 100644 --- a/test/test_field_evaluation.cpp +++ b/test/test_field_evaluation.cpp @@ -1,140 +1,86 @@ #include -#include -#include #include #include -#include "pcms/adapter/meshfields/mesh_fields_adapter.h" -#include "pcms/adapter/meshfields/mesh_fields_adapter2.h" -#include "pcms/create_field.h" -#include -#include +#include +#include "pcms/lagrange_field_factory.h" +#include "pcms/utility/assert.h" +#include "field_test_utils.h" +#include using pcms::Real; -TEST_CASE("evaluate linear 2d omega_h_field") +// Non-linear test function used to exercise quadratic interpolation. +KOKKOS_INLINE_FUNCTION static Real sin_f(Real x, Real y) { - auto lib = Omega_h::Library{}; - auto world = lib.world(); - auto mesh = - Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false); - auto layout = - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); - const auto nverts = mesh.nents(0); - auto mesh_coords = mesh.coords(); - auto f = KOKKOS_LAMBDA(Real x, Real y) - { - return std::sin(20 * x * y) / 2 + 0.5; - }; - Omega_h::Write test_f(nverts); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = f(x, y); - }); - Omega_h::HostWrite test_f_host(test_f); - auto field = layout->CreateFieldReal(); - field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); - - std::vector coords = { - 0.7681, 0.886, 0.5337, 0.5205, 0.8088, 0.1513, 0.13, - 0.43, 0.5484, 0.8263, 0.006119, 0.8642, 0.5889, 0.5622, - 0.9268, 0.1749, 0.2615, 0.1468, 0.9793, 0.9612, - }; - - std::vector evaluation(coords.size() / 2); - pcms::Rank1View eval_view{evaluation.data(), - evaluation.size()}; - pcms::Rank2View coords_view( - coords.data(), coords.size() / 2, 2); - pcms::FieldDataView data_view( - eval_view, field->GetCoordinateSystem()); - pcms::CoordinateView coordinate_view{ - field->GetCoordinateSystem(), coords_view}; - - auto locale = field->GetLocalizationHint(coordinate_view); - field->Evaluate(locale, data_view); - - for (int i = 0; i < coords.size() / 2; ++i) { - Real x = coords[2 * i + 0]; - Real y = coords[2 * i + 1]; + return std::sin(20 * x * y) / 2 + 0.5; +} - Real test_value = evaluation[i]; - Real reference_value = f(x, y); - Real percent_error = - 100 * std::abs(test_value - reference_value) / reference_value; +// Standard interior test points shared across all evaluation tests. +static const std::vector kEvalCoords = { + 0.7681, 0.886, 0.5337, 0.5205, 0.8088, 0.1513, + 0.13, 0.43, 0.5484, 0.8263, 0.006119, 0.8642, + 0.5889, 0.5622, 0.9268, 0.1749, 0.2615, 0.1468, + 0.9793, 0.9612, +}; - REQUIRE(percent_error < 1.0); - } +TEST_CASE("evaluate linear 2d omega_h_field") +{ + auto lib = Omega_h::Library{}; + auto mesh = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, + 100, 100, 0, false); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto field = factory.CreateFieldReal(); + + pcms::test::SetField(*field, pcms::test::linear_f); + pcms::test::CheckEvaluation( + *field, pcms::test::StandardEvalCoords2D(), pcms::test::linear_f); } -TEST_CASE("evaluate quadratic 2d omega_h_field") +#ifdef PCMS_ENABLE_MESHFIELDS +TEST_CASE("evaluate quadratic 2d meshfields_field") { auto lib = Omega_h::Library{}; - auto world = lib.world(); - auto mesh = - Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false); - auto layout = - pcms::CreateLagrangeLayout(mesh, 2, 1, pcms::CoordinateSystem::Cartesian); + auto mesh = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, + 100, 100, 0, false); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::MeshFields); + + // Quadratic DOF holders span vertices and edge midpoints; set them inline. const auto nverts = mesh.nents(0); const auto nedges = mesh.nents(1); auto mesh_coords = mesh.coords(); auto edge_verts = mesh.ask_verts_of(1); - auto f = KOKKOS_LAMBDA(Real x, Real y) - { - return std::sin(20 * x * y) / 2 + 0.5; - }; + Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = f(x, y); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto endpoints = Omega_h::gather_verts<2>(edge_verts, i); - Real x0 = mesh_coords[2 * endpoints[0] + 0]; - Real y0 = mesh_coords[2 * endpoints[0] + 1]; - Real x1 = mesh_coords[2 * endpoints[1] + 0]; - Real y1 = mesh_coords[2 * endpoints[1] + 1]; - Real cx = (x0 + x1) / 2; - Real cy = (y0 + y1) / 2; - test_f[nverts + i] = f(cx, cy); - }); + Omega_h::parallel_for(nverts, OMEGA_H_LAMBDA(int i) { + test_f[i] = sin_f(mesh_coords[2 * i], mesh_coords[2 * i + 1]); + }); + Omega_h::parallel_for(nedges, OMEGA_H_LAMBDA(int i) { + auto ep = Omega_h::gather_verts<2>(edge_verts, i); + Real cx = (mesh_coords[2 * ep[0]] + mesh_coords[2 * ep[1]]) / 2; + Real cy = (mesh_coords[2 * ep[0] + 1] + mesh_coords[2 * ep[1] + 1]) / 2; + test_f[nverts + i] = sin_f(cx, cy); + }); Omega_h::HostWrite test_f_host(test_f); - auto field = layout->CreateFieldReal(); + auto field = factory.CreateFieldReal(); field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); - std::vector coords = { - 0.7681, 0.886, 0.5337, 0.5205, 0.8088, 0.1513, 0.13, - 0.43, 0.5484, 0.8263, 0.006119, 0.8642, 0.5889, 0.5622, - 0.9268, 0.1749, 0.2615, 0.1468, 0.9793, 0.9612, - }; - - std::vector evaluation(coords.size() / 2); - pcms::Rank1View eval_view{evaluation.data(), - evaluation.size()}; - pcms::Rank2View coords_view( - coords.data(), coords.size() / 2, 2); - pcms::FieldDataView data_view( - eval_view, field->GetCoordinateSystem()); - pcms::CoordinateView coordinate_view{ - field->GetCoordinateSystem(), coords_view}; - - auto locale = field->GetLocalizationHint(coordinate_view); - field->Evaluate(locale, data_view); - - for (int i = 0; i < coords.size() / 2; ++i) { - Real x = coords[2 * i + 0]; - Real y = coords[2 * i + 1]; - - Real test_value = evaluation[i]; - Real reference_value = f(x, y); - Real percent_error = - 100 * std::abs(test_value - reference_value) / reference_value; + pcms::test::CheckEvaluation(*field, kEvalCoords, sin_f, 1.0e-2); +} +#endif - REQUIRE(percent_error < 1.0); - } +TEST_CASE("evaluate quadratic 2d omega_h_field throws") +{ + auto lib = Omega_h::Library{}; + auto mesh = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1, 1, 0, + 100, 100, 0, false); + + REQUIRE_THROWS_AS( + pcms::LagrangeFieldFactory::FromMesh( + mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::OmegaH), + pcms::pcms_error); } diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp index 413f06ca..2190e0b6 100644 --- a/test/test_field_interpolation.cpp +++ b/test/test_field_interpolation.cpp @@ -7,37 +7,30 @@ #include #include "pcms/adapter/meshfields/mesh_fields_adapter.h" #include "pcms/adapter/meshfields/mesh_fields_adapter2.h" -#include "pcms/create_field.h" +#include "pcms/lagrange_field_factory.h" +#include "pcms/utility/assert.h" +#include "field_test_utils.h" #include #include using pcms::Real; +KOKKOS_INLINE_FUNCTION static Real interpolation_linear_f(Real x, Real y) +{ + return -0.3 * x + 0.5 * y; +} + TEST_CASE("interpolate linear 2d omega_h_field") { auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false); - auto layout = - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); - const auto nverts = mesh.nents(0); - auto mesh_coords = mesh.coords(); - auto f = KOKKOS_LAMBDA(Real x, Real y) - { - return -0.3 * x + 0.5 * y; - }; - Omega_h::Write test_f(nverts); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = f(x, y); - }); - Omega_h::HostWrite test_f_host(test_f); - auto field = layout->CreateFieldReal(); - auto interpolated = layout->CreateFieldReal(); - field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto field = factory.CreateFieldReal(); + auto interpolated = factory.CreateFieldReal(); + pcms::test::SetField(*field, interpolation_linear_f); pcms::interpolate_field2(*field, *interpolated); auto interpolated_dof = interpolated->GetDOFHolderData(); @@ -51,28 +44,27 @@ TEST_CASE("interpolate linear 2d omega_h_field") } } -TEST_CASE("interpolate quadratic 2d omega_h_field") +#ifdef PCMS_ENABLE_MESHFIELDS +TEST_CASE("interpolate quadratic 2d meshfields_field") { auto lib = Omega_h::Library{}; auto world = lib.world(); auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false); - auto layout = - pcms::CreateLagrangeLayout(mesh, 2, 1, pcms::CoordinateSystem::Cartesian); + auto factory2 = pcms::LagrangeFieldFactory::FromMesh( + mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::MeshFields); + auto layout = factory2.GetLayout(); const auto nverts = mesh.nents(0); const auto nedges = mesh.nents(1); auto mesh_coords = mesh.coords(); auto edge_verts = mesh.ask_verts_of(1); - auto f = KOKKOS_LAMBDA(Real x, Real y) - { - return -0.3 * x + 0.5 * y; - }; Omega_h::Write test_f(nverts + nedges); Omega_h::parallel_for( nverts, OMEGA_H_LAMBDA(int i) { Real x = mesh_coords[2 * i + 0]; Real y = mesh_coords[2 * i + 1]; - test_f[i] = f(x, y); + test_f[i] = interpolation_linear_f(x, y); }); Omega_h::parallel_for( nedges, OMEGA_H_LAMBDA(int i) { @@ -83,12 +75,12 @@ TEST_CASE("interpolate quadratic 2d omega_h_field") Real y1 = mesh_coords[2 * endpoints[1] + 1]; Real cx = (x0 + x1) / 2; Real cy = (y0 + y1) / 2; - test_f[nverts + i] = f(cx, cy); + test_f[nverts + i] = interpolation_linear_f(cx, cy); }); Omega_h::HostWrite test_f_host(test_f); - auto field = layout->CreateFieldReal(); - auto interpolated = layout->CreateFieldReal(); + auto field = factory2.CreateFieldReal(); + auto interpolated = factory2.CreateFieldReal(); field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); // interpolate the field from one mesh to another mesh with the same @@ -105,3 +97,18 @@ TEST_CASE("interpolate quadratic 2d omega_h_field") Catch::Matchers::WithinAbs(original_dof[i], 1E-10)); } } +#endif + +TEST_CASE("interpolate quadratic 2d omega_h_field throws") +{ + auto lib = Omega_h::Library{}; + auto world = lib.world(); + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 100, 100, 0, false); + + REQUIRE_THROWS_AS( + pcms::LagrangeFieldFactory::FromMesh( + mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFieldFactory::Backend::OmegaH), + pcms::pcms_error); +} diff --git a/test/test_omega_h_field2_outofbounds.cpp b/test/test_omega_h_field2_outofbounds.cpp index ab77f05f..9d07cd0b 100644 --- a/test/test_omega_h_field2_outofbounds.cpp +++ b/test/test_omega_h_field2_outofbounds.cpp @@ -3,13 +3,18 @@ #include #include #include -#include "pcms/adapter/meshfields/mesh_fields_adapter2.h" -#include "pcms/create_field.h" +#include "pcms/lagrange_field_factory.h" +#include "field_test_utils.h" #include #include using pcms::Real; +KOKKOS_INLINE_FUNCTION static Real fill_mode_f(Real x, Real y) +{ + return x + y; +} + TEST_CASE("omega_h_field2 out of bounds FILL mode") { auto lib = Omega_h::Library{}; @@ -17,26 +22,10 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") // Create a 1x1 box mesh (coords from 0 to 1) auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 0, 10, 10, 0, false); - auto layout = - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); - const auto nverts = mesh.nents(0); - auto mesh_coords = mesh.coords(); - - // Set up a simple linear field - auto f = KOKKOS_LAMBDA(Real x, Real y) - { - return x + y; - }; - Omega_h::Write test_f(nverts); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = f(x, y); - }); - Omega_h::HostWrite test_f_host(test_f); - auto field = layout->CreateFieldReal(); - field->SetDOFHolderData(pcms::make_const_array_view(test_f_host)); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto field = factory.CreateFieldReal(); + pcms::test::SetField(*field, fill_mode_f); // Set FILL mode with fill value of -999.0 Real fill_value = -999.0; @@ -51,32 +40,7 @@ TEST_CASE("omega_h_field2 out of bounds FILL mode") -0.1, 0.5, // outside (x < 0) - should return fill_value }; - std::vector evaluation(coords.size() / 2); - pcms::Rank1View eval_view{evaluation.data(), - evaluation.size()}; - pcms::Rank2View coords_view( - coords.data(), coords.size() / 2, 2); - pcms::FieldDataView data_view( - eval_view, field->GetCoordinateSystem()); - pcms::CoordinateView coordinate_view{ - field->GetCoordinateSystem(), coords_view}; - - auto locale = field->GetLocalizationHint(coordinate_view); - field->Evaluate(locale, data_view); - - // Check results - // Point 0: (0.5, 0.5) - inside, should be close to f(0.5, 0.5) = 1.0 - REQUIRE(std::abs(evaluation[0] - 1.0) < 0.1); - - // Point 1: (1.5, 0.5) - outside, should be fill_value - REQUIRE(evaluation[1] == fill_value); - - // Point 2: (0.5, -0.1) - outside, should be fill_value - REQUIRE(evaluation[2] == fill_value); - - // Point 3: (0.3, 0.7) - inside, should be close to f(0.3, 0.7) = 1.0 - REQUIRE(std::abs(evaluation[3] - 1.0) < 0.1); - - // Point 4: (-0.1, 0.5) - outside, should be fill_value - REQUIRE(evaluation[4] == fill_value); -} \ No newline at end of file + std::vector is_inside = {true, false, false, true, false}; + pcms::test::CheckEvaluationWithFill( + *field, coords, is_inside, fill_mode_f, fill_value, 1.0e-10); +} diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp new file mode 100644 index 00000000..030a0475 --- /dev/null +++ b/test/test_omega_h_lagrange_field.cpp @@ -0,0 +1,310 @@ +#include +#include +#include +#include +#include + +#include "pcms/adapter/omega_h/omega_h_lagrange_layout.h" +#include "pcms/adapter/omega_h/omega_h_lagrange_field.h" +#include "pcms/lagrange_field_factory.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/mesh_geometry.h" +#include "field_test_utils.h" + +#include +#include + +using pcms::Real; +using pcms::LO; + +static Omega_h::Mesh MakeBox2D(Omega_h::CommPtr world, int nx = 10, int ny = 10) +{ + return Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, nx, ny, 0, + false); +} + +// ---- Layout tests ----------------------------------------------------------- + +TEST_CASE("OmegaHLagrangeLayout order-1 properties") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + pcms::OmegaHLagrangeLayout layout(mesh, 1, 2, + pcms::CoordinateSystem::Cartesian); + + REQUIRE(layout.GetOrder() == 1); + REQUIRE(layout.GetNumComponents() == 2); + REQUIRE(layout.GetNumOwnedDofHolder() == mesh.nents(0)); + REQUIRE(layout.GetNumGlobalDofHolder() == mesh.nglobal_ents(0)); + REQUIRE(layout.IsDistributed()); // always true for Omega_h mesh layouts + + // DOF holder coordinates should match vertex coordinates + auto coords = layout.GetDOFHolderCoordinates().GetCoordinates(); + auto mesh_coords = Omega_h::HostRead(mesh.coords()); + int nverts = mesh.nents(0); + REQUIRE(static_cast(coords.extent(0)) == nverts); + REQUIRE(static_cast(coords.extent(1)) == mesh.dim()); + for (int v = 0; v < nverts; ++v) { + for (int d = 0; d < mesh.dim(); ++d) { + REQUIRE(coords(v, d) == + Catch::Approx(mesh_coords[v * mesh.dim() + d])); + } + } + + // GetEntOffsets: vertices are at slot 0, all other slots = nverts + auto offsets = layout.GetEntOffsets(); + REQUIRE(offsets[0] == 0); + for (int i = 1; i < pcms::ent_offsets_len; ++i) + REQUIRE(offsets[i] == nverts); +} + +TEST_CASE("OmegaHLagrangeLayout order-0 properties") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + pcms::OmegaHLagrangeLayout layout(mesh, 0, 1, + pcms::CoordinateSystem::Cartesian); + + REQUIRE(layout.GetOrder() == 0); + REQUIRE(layout.GetNumComponents() == 1); + REQUIRE(layout.GetNumOwnedDofHolder() == mesh.nelems()); + REQUIRE(layout.GetNumGlobalDofHolder() == mesh.nglobal_ents(mesh.dim())); + + // DOF holder coordinates should match element centroids + auto coords = layout.GetDOFHolderCoordinates().GetCoordinates(); + auto centroids = + Omega_h::HostRead(pcms::get_entity_centroids(mesh, mesh.dim())); + int nelems = mesh.nelems(); + REQUIRE(static_cast(coords.extent(0)) == nelems); + for (int e = 0; e < nelems; ++e) { + for (int d = 0; d < mesh.dim(); ++d) { + REQUIRE(coords(e, d) == + Catch::Approx(centroids[e * mesh.dim() + d])); + } + } + + // GetEntOffsets: all DOFs are at entity_dim = mesh.dim() (slot 2 for 2D) + auto offsets = layout.GetEntOffsets(); + for (int i = 0; i <= mesh.dim(); ++i) + REQUIRE(offsets[i] == 0); + for (int i = mesh.dim() + 1; i < pcms::ent_offsets_len; ++i) + REQUIRE(offsets[i] == nelems); +} + +TEST_CASE("OmegaHLagrangeLayout invalid order throws") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + REQUIRE_THROWS_AS( + pcms::OmegaHLagrangeLayout(mesh, 2, 1, pcms::CoordinateSystem::Cartesian), + std::invalid_argument); + REQUIRE_THROWS_AS( + pcms::OmegaHLagrangeLayout(mesh, -1, 1, pcms::CoordinateSystem::Cartesian), + std::invalid_argument); +} + +TEST_CASE("OmegaHLagrangeLayout layout sharing") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + + auto f1 = std::make_unique>(layout); + auto f2 = std::make_unique>(layout); + + REQUIRE(&f1->GetLayout() == &f2->GetLayout()); +} + +// ---- Order-1 field tests ---------------------------------------------------- + +TEST_CASE("OmegaHLagrangeField order-1: set/get DOF data round-trip") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + int n = layout->GetNumOwnedDofHolder(); + std::vector data(n); + for (int i = 0; i < n; ++i) + data[i] = static_cast(i); + + pcms::Rank1View view(data.data(), n); + field.SetDOFHolderData(view); + + auto got = field.GetDOFHolderData(); + REQUIRE(static_cast(got.size()) == n); + for (int i = 0; i < n; ++i) + REQUIRE(got[i] == Catch::Approx(data[i])); +} + +TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world(), 20, 20); + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + pcms::test::SetField(field, pcms::test::linear_f); + pcms::test::CheckEvaluation( + field, pcms::test::StandardEvalCoords2D(), pcms::test::linear_f); +} + +// The same linear evaluation test run on the MeshFields-backed order-1 field +// ensures both backends produce identical results for the same inputs. +TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world(), 20, 20); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto field = factory.CreateFieldReal(); + + pcms::test::SetField(*field, pcms::test::linear_f); + pcms::test::CheckEvaluation( + *field, pcms::test::StandardEvalCoords2D(), pcms::test::linear_f); +} + +TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + pcms::test::SetField(field, pcms::test::linear_f); + + Real fill_value = -999.0; + field.SetOutOfBoundsMode(pcms::OutOfBoundsMode::FILL, fill_value); + + // Points clearly outside the [0,1]^2 domain + std::vector outside{-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; + pcms::test::CheckFillMode(field, fill_value, outside); +} + +TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + pcms::test::SetField(field, pcms::test::linear_f); + pcms::test::CheckSerializeDeserialize(field); +} + +// ---- Order-0 field tests ---------------------------------------------------- + +TEST_CASE("OmegaHLagrangeField order-0: set/get DOF data round-trip") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + int n = layout->GetNumOwnedDofHolder(); + std::vector data(n, 3.14); + pcms::Rank1View view(data.data(), n); + field.SetDOFHolderData(view); + + auto got = field.GetDOFHolderData(); + REQUIRE(static_cast(got.size()) == n); + for (int i = 0; i < n; ++i) + REQUIRE(got[i] == Catch::Approx(3.14)); +} + +TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world(), 10, 10); + auto layout = std::make_shared( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + // Set every element to the same value + const Real kValue = 42.0; + int nelems = mesh.nelems(); + std::vector data(nelems, kValue); + pcms::Rank1View view(data.data(), nelems); + field.SetDOFHolderData(view); + + // All evaluation points should return kValue + auto pts = pcms::test::StandardEvalCoords2D(); + int n = static_cast(pts.size()) / 2; + std::vector eval(n); + pcms::Rank2View coords_view(pts.data(), n, + 2); + pcms::Rank1View eval_view(eval.data(), n); + pcms::CoordinateView cv{field.GetCoordinateSystem(), + coords_view}; + pcms::FieldDataView dv{eval_view, + field.GetCoordinateSystem()}; + auto hint = field.GetLocalizationHint(cv); + field.Evaluate(hint, dv); + + for (int i = 0; i < n; ++i) + REQUIRE(eval[i] == Catch::Approx(kValue)); +} + +TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + int nelems = mesh.nelems(); + std::vector data(nelems, 1.0); + pcms::Rank1View view(data.data(), nelems); + field.SetDOFHolderData(view); + + Real fill_value = -1.0; + field.SetOutOfBoundsMode(pcms::OutOfBoundsMode::FILL, fill_value); + + std::vector outside{-0.5, 0.5, 1.5, 0.5}; + pcms::test::CheckFillMode(field, fill_value, outside); +} + +TEST_CASE("OmegaHLagrangeField order-0: serialize / deserialize round-trip") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + auto layout = std::make_shared( + mesh, 0, 1, pcms::CoordinateSystem::Cartesian); + pcms::OmegaHLagrangeField field(layout); + + int nelems = mesh.nelems(); + std::vector data(nelems); + for (int i = 0; i < nelems; ++i) + data[i] = static_cast(i); + pcms::Rank1View view(data.data(), nelems); + field.SetDOFHolderData(view); + + pcms::test::CheckSerializeDeserialize(field); +} + +// ---- Temporary factory lifetime safety -------------------------------------- + +TEST_CASE("OmegaHLagrangeField: field valid after layout destruction") +{ + auto lib = Omega_h::Library{}; + auto mesh = MakeBox2D(lib.world()); + + std::unique_ptr> field; + { + auto layout = std::make_shared( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + field = std::make_unique>(layout); + } // layout shared_ptr goes out of scope here; field keeps it alive + + pcms::test::SetField(*field, pcms::test::linear_f); + pcms::test::CheckEvaluation( + *field, pcms::test::StandardEvalCoords2D(), pcms::test::linear_f); +} diff --git a/test/test_proxy_coupling.cpp b/test/test_proxy_coupling.cpp index 6f8a3ffc..d66a086c 100644 --- a/test/test_proxy_coupling.cpp +++ b/test/test_proxy_coupling.cpp @@ -7,7 +7,7 @@ #include #include "test_support.h" #include "pcms/coupler2.h" -#include "pcms/create_field.h" +#include "pcms/lagrange_field_factory.h" #include #include @@ -91,12 +91,12 @@ void xgc_delta_f(MPI_Comm comm, Omega_h::Mesh& mesh) pcms::Coupler2 coupler("proxy_couple", comm, false, {}); pcms::Application2* app = coupler.AddApplication("proxy_couple_xgc_delta_f"); - auto& layout = app->AddLayout( - "gids", - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian)); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto& layout = app->AddLayout("gids", factory.GetLayout()); - auto gids_field = layout.CreateFieldReal(); - auto gids2_field = layout.CreateFieldReal(); + auto gids_field = factory.CreateFieldReal(); + auto gids2_field = factory.CreateFieldReal(); auto* gids_ptr = gids_field.get(); auto* gids2_ptr = gids2_field.get(); @@ -132,11 +132,11 @@ void xgc_total_f(MPI_Comm comm, Omega_h::Mesh& mesh) pcms::Coupler2 coupler("proxy_couple", comm, false, {}); pcms::Application2* app = coupler.AddApplication("proxy_couple_xgc_total_f"); - auto& layout = app->AddLayout( - "gids", - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian)); + auto factory = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto& layout = app->AddLayout("gids", factory.GetLayout()); - auto gids_field = layout.CreateFieldReal(); + auto gids_field = factory.CreateFieldReal(); auto* gids_ptr = gids_field.get(); @@ -175,17 +175,17 @@ void xgc_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, std::string_view cpn_file) const auto partition = std::get(cpl.GetPartition()); auto* total_f = cpl.AddApplication("proxy_couple_xgc_total_f"); auto* delta_f = cpl.AddApplication("proxy_couple_xgc_delta_f"); - auto& layout_total = total_f->AddLayout( - "gids", - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian)); + auto factory_total = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto& layout_total = total_f->AddLayout("gids", factory_total.GetLayout()); - auto& layout_delta = delta_f->AddLayout( - "gids", - pcms::CreateLagrangeLayout(mesh, 1, 1, pcms::CoordinateSystem::Cartesian)); + auto factory_delta = pcms::LagrangeFieldFactory::FromMesh( + mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto& layout_delta = delta_f->AddLayout("gids", factory_delta.GetLayout()); // TODO, fields should have a transfer policy rather than parameters - auto total_gids_field = layout_total.CreateFieldReal(); - auto delta_gids_field = layout_delta.CreateFieldReal(); - auto delta_gids2_field = layout_delta.CreateFieldReal(); + auto total_gids_field = factory_total.CreateFieldReal(); + auto delta_gids_field = factory_delta.CreateFieldReal(); + auto delta_gids2_field = factory_delta.CreateFieldReal(); auto* total_gids_ptr = total_gids_field.get(); auto* delta_gids_ptr = delta_gids_field.get(); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 13c96544..4ecdd391 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -9,27 +9,14 @@ #include "pcms/transfer_field2.h" #include "pcms/adapter/meshfields/mesh_fields_adapter_layout.h" #include "pcms/adapter/meshfields/mesh_fields_adapter2.h" -#include "pcms/create_field.h" +#include "pcms/adapter/uniform_grid/uniform_grid_binary_field.h" +#include "pcms/lagrange_field_factory.h" #include "pcms/utility/arrays.h" +#include "field_test_utils.h" #include -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, @@ -78,19 +65,71 @@ TEST_CASE("UniformGrid field creation") grid.divisions = {5, 5}; // Create field layout with 1 component (scalar field) - pcms::UniformGridFieldLayout<2> layout(grid, 1, - pcms::CoordinateSystem::Cartesian); + auto layout = std::make_shared>( + 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()); + 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.CreateFieldReal(); + auto field = std::make_unique>(layout); REQUIRE(field != nullptr); } +TEST_CASE("UniformGrid order-0 field creation and evaluation") +{ + pcms::UniformGrid<2> grid; + grid.bot_left = {0.0, 0.0}; + grid.edge_length = {10.0, 10.0}; + grid.divisions = {2, 2}; + + auto layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian, 0); + auto field = std::make_unique>(layout); + + REQUIRE(layout->GetOrder() == 0); + REQUIRE(layout->GetNumOwnedDofHolder() == 4); + REQUIRE(layout->GetNumGlobalDofHolder() == 4); + + auto coords = layout->GetDOFHolderCoordinates().GetCoordinates(); + REQUIRE(coords(0, 0) == Catch::Approx(2.5)); + REQUIRE(coords(0, 1) == Catch::Approx(2.5)); + REQUIRE(coords(3, 0) == Catch::Approx(7.5)); + REQUIRE(coords(3, 1) == Catch::Approx(7.5)); + + std::vector data = {1.0, 2.0, 3.0, 4.0}; + field->SetDOFHolderData( + pcms::Rank1View( + data.data(), data.size())); + + std::vector eval_coords = { + 1.0, 1.0, + 9.0, 1.0, + 1.0, 9.0, + 9.0, 9.0 + }; + 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); + + REQUIRE(results[0] == Catch::Approx(1.0)); + REQUIRE(results[1] == Catch::Approx(2.0)); + REQUIRE(results[2] == Catch::Approx(3.0)); + REQUIRE(results[3] == Catch::Approx(4.0)); +} + TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") { pcms::UniformGrid<2> grid; @@ -98,9 +137,9 @@ TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") 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.CreateFieldReal(); + auto layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto field = std::make_unique>(layout); // Initialize field data with vertex indices (5x5 = 25 vertices for 4x4 cells) std::vector data(25); @@ -128,9 +167,9 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant") 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.CreateFieldReal(); + auto layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto field = std::make_unique>(layout); // Set vertex values for a 2x2 cell grid (3x3 = 9 vertices) // Vertex layout: @@ -190,9 +229,9 @@ TEST_CASE("UniformGrid field serialization") 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.CreateFieldReal(); + auto layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto field = std::make_unique>(layout); // Set field data (4x4 = 16 vertices for 3x3 cells) std::vector data(16); @@ -204,39 +243,7 @@ TEST_CASE("UniformGrid field serialization") 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.CreateFieldReal(); - 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]); - } + pcms::test::CheckSerializeDeserialize(*field); } TEST_CASE("UniformGrid field copy") @@ -246,9 +253,9 @@ TEST_CASE("UniformGrid field copy") 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.CreateFieldReal(); + auto layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto field = std::make_unique>(layout); // 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), @@ -292,7 +299,7 @@ TEST_CASE("UniformGrid field copy") REQUIRE(std::abs(results[3] - 15.0) < 1e-10); // Test 2: test copy_field2 - auto field2 = layout.CreateFieldReal(); + auto field2 = std::make_unique>(layout); pcms::copy_field2(*field, *field2); auto copied_data = field2->GetDOFHolderData(); @@ -311,20 +318,12 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") 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->CreateFieldReal(); + auto omega_h_factory = + pcms::LagrangeFieldFactory::FromMesh(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto omega_h_field = omega_h_factory.CreateFieldReal(); // Initialize omega_h field with a simple function f(x,y) = x + 2*y - auto coords = omega_h_layout->GetDOFHolderCoordinates(); - int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); - std::vector omega_h_data = - CreateOmegaHFieldData(coords, num_nodes); - - auto omega_h_data_view = - pcms::Rank1View( - omega_h_data.data(), omega_h_data.size()); - omega_h_field->SetDOFHolderData(omega_h_data_view); + pcms::test::SetField(*omega_h_field, pcms::test::linear_f); // Create a uniform grid field covering the same domain [0,1] x [0,1] pcms::UniformGrid<2> grid; @@ -332,12 +331,12 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") 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.CreateFieldReal(); + auto ug_layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto ug_field = std::make_unique>(ug_layout); // Transfer from omega_h field to uniform grid field using interpolation - auto coords_interpolation = ug_layout.GetDOFHolderCoordinates(); + auto coords_interpolation = ug_layout->GetDOFHolderCoordinates(); std::vector evaluation( coords_interpolation.GetCoordinates().size() / 2); auto evaluation_view = pcms::Rank1View( @@ -353,9 +352,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") // Verify the transferred data at uniform grid vertices auto transferred_data = ug_field->GetDOFHolderData(); - auto ug_coords = ug_layout.GetDOFHolderCoordinates(); + auto ug_coords = ug_layout->GetDOFHolderCoordinates(); auto ug_coords_data = ug_coords.GetCoordinates(); - int num_ug_nodes = ug_layout.GetNumOwnedDofHolder(); // 5x5 = 25 vertices + int num_ug_nodes = ug_layout->GetNumOwnedDofHolder(); // 5x5 = 25 vertices // Check a few sample points for (int i = 0; i < num_ug_nodes; ++i) { @@ -376,19 +375,18 @@ TEST_CASE("Create binary field from uniform grid") 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {5, 5}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{5, 5}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 36); // (5+1) * (5+1) = 36 vertices - // All grid vertices should be inside the mesh pcms::Real sum = 0.0; for (size_t i = 0; i < field_data.extent(0); ++i) { + REQUIRE((field_data(i) == 0.0 || field_data(i) == 1.0)); sum += field_data(i); } REQUIRE(sum == 36.0); @@ -399,35 +397,75 @@ 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 8}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{10, 8}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 99); // (10+1) * (8+1) = 99 vertices - // Most vertices should be inside pcms::Real inside_count = 0.0; - for (size_t i = 0; i < field_data.extent(0); ++i) { + 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") + 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); + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, 8, 8, 0, false); - auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 8); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{10, 10}); auto field_data = field->GetDOFHolderData(); - REQUIRE(field_data.extent(0) == 81); // (8+1) * (8+1) = 81 vertices + for (size_t i = 0; i < field_data.extent(0); ++i) { + REQUIRE((field_data(i) == 0.0 || field_data(i) == 1.0)); + } + } + SECTION("Grid larger than mesh - vertices outside should be marked 0") + { + auto mesh = + Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.5, 0.5, 0.0, 5, 5, 0, false); + + pcms::UniformGrid<2> grid; + grid.edge_length = {1.0, 1.0}; + grid.bot_left = {0.0, 0.0}; + grid.divisions = {10, 10}; + + auto [layout, field] = pcms::CreateUniformGridBinaryField<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); } SECTION("Fine grid over coarse mesh") @@ -437,8 +475,9 @@ TEST_CASE("Create binary field from uniform grid") 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {20, 20}); + auto grid = CreateUniformGridFromMesh<2>(mesh, std::array{20, 20}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{20, 20}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 441); // (20+1) * (20+1) = 441 vertices @@ -453,50 +492,14 @@ TEST_CASE("Create binary field from uniform grid") REQUIRE(field_data(corner_idx) == 1.0); } - 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10); - auto field_data = field->GetDOFHolderData(); - - // All values should be 0 or 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 - vertices 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {10, 10}); - auto field_data = field->GetDOFHolderData(); - - REQUIRE(field_data.extent(0) == 121); // (10+1) * (10+1) = 121 vertices - - // 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") { // 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {30, 10}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{30, 10}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 341); // (30+1) * (10+1) = 341 vertices @@ -524,8 +527,9 @@ TEST_CASE("Binary field integration with grid methods") SECTION("Query field value at specific grid vertex") { - auto grid = CreateUniformGridFromMesh<2>(mesh, {8, 8}); - auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {8, 8}); + auto grid = CreateUniformGridFromMesh<2>(mesh, std::array{8, 8}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{8, 8}); auto field_data = field->GetDOFHolderData(); // Get field value for a specific vertex (middle vertex at i=4, j=4) @@ -545,8 +549,9 @@ TEST_CASE("Binary field integration with grid methods") SECTION("Count vertices by region") { - auto grid = CreateUniformGridFromMesh<2>(mesh, 10); - auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, 10); + auto grid = CreateUniformGridFromMesh<2>(mesh, std::array{10, 10}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{10, 10}); auto field_data = field->GetDOFHolderData(); // Count vertices in different quadrants @@ -603,7 +608,8 @@ 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, 50); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{50, 50}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 2601); // (50+1) * (50+1) = 2601 vertices @@ -621,7 +627,8 @@ TEST_CASE("Performance and edge cases") 10, 0, false); // Very coarse grid - auto [layout, field] = CreateUniformGridBinaryField<2>(mesh, {2, 2}); + auto [layout, field] = + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{2, 2}); auto field_data = field->GetDOFHolderData(); REQUIRE(field_data.extent(0) == 9); // (2+1) * (2+1) = 9 vertices @@ -639,64 +646,17 @@ 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 [layout, field] = CreateUniformGridBinaryField<2>(mesh, {25, 10}); - auto field_data = field->GetDOFHolderData(); - - REQUIRE(field_data.extent(0) == 286); // (25+1) * (10+1) = 286 vertices - - 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 - vertices 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 [layout, field] = - pcms::CreateUniformGridBinaryFieldFromGrid<2>(mesh, grid); + pcms::CreateUniformGridBinaryField<2>(mesh, std::array{25, 10}); auto field_data = field->GetDOFHolderData(); - REQUIRE(field_data.extent(0) == 121); // (10+1) * (10+1) = 121 vertices + REQUIRE(field_data.extent(0) == 286); // (25+1) * (10+1) = 286 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); } } @@ -709,33 +669,26 @@ 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_layout, 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->CreateFieldReal(); + auto omega_h_factory = + pcms::LagrangeFieldFactory::FromMesh(mesh, 1, 1, pcms::CoordinateSystem::Cartesian); + auto omega_h_field = omega_h_factory.CreateFieldReal(); // Initialize omega_h field with a simple function f(x,y) = x + 2*y - auto coords = omega_h_layout->GetDOFHolderCoordinates(); - int num_nodes = omega_h_layout->GetNumOwnedDofHolder(); - std::vector omega_h_data = - CreateOmegaHFieldData(coords, num_nodes); - - auto omega_h_data_view = - pcms::Rank1View( - omega_h_data.data(), omega_h_data.size()); - omega_h_field->SetDOFHolderData(omega_h_data_view); + pcms::test::SetField(*omega_h_field, pcms::test::linear_f); // Create uniform grid field layout - pcms::UniformGridFieldLayout<2> ug_layout(grid, 1, - pcms::CoordinateSystem::Cartesian); - auto ug_field = ug_layout.CreateFieldReal(); + auto ug_layout = std::make_shared>( + grid, 1, pcms::CoordinateSystem::Cartesian); + auto ug_field = std::make_unique>(ug_layout); + + auto [mask_layout, mask_field] = + pcms::CreateUniformGridBinaryField<2>(mesh, grid); // 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(); + auto ug_coords = ug_layout->GetDOFHolderCoordinates(); // Verify ug_field values directly from the field object auto ug_field_data = ug_field->GetDOFHolderData();