Mesh intersection routines

CMakeLists.txt

@@ -61,7 +61,36 @@ endif()
 find_package(ADIOS2 CONFIG 2.5 REQUIRED)
 find_package(Kokkos CONFIG 4.5 REQUIRED)
 
-find_package(meshfields REQUIRED)
+option(PCMS_ENABLE_MESHFIELDS "Enable MeshFields support" ON)
+option(PCMS_ENABLE_PETSC "Enable PETSc support" OFF)
+
+if (PCMS_ENABLE_MESHFIELDS)
+	find_package(meshfields REQUIRED)
+	message(STATUS "Found MeshFields: ${meshfields_DIR} (found version ${meshfields_VERSION})")
+endif()
+
+if (PCMS_ENABLE_PETSC)
+	if (NOT PETSC_DIR OR NOT PETSC_ARCH)
+		message(FATAL_ERROR "PETSC_DIR and PETSC_ARCH must be provided when PETSc is not CMake-enabled")
+	endif()
+	message(STATUS "PETSC_DIR = ${PETSC_DIR}")
+	message(STATUS "PETSC_ARCH = ${PETSC_ARCH}")
+
+	set(PETSC_INCLUDE_DIRS "${PETSC_DIR}/include" "${PETSC_DIR}/${PETSC_ARCH}/include")
+	set(PETSC_LIB_DIR "${PETSC_DIR}/${PETSC_ARCH}/lib")
+
+	message(STATUS "PETSC_LIB_DIR = ${PETSC_LIB_DIR}")
+	find_library(PETSC_LIBRARY petsc HINTS ${PETSC_LIB_DIR})
+
+	if (NOT PETSC_LIBRARY)
+		message(FATAL_ERROR "Couldnot find libpetsc in ${PETSC_LIB_DIR}")
+	endif()
+
+	set(PCMS_PETSC_INCLUDES ${PETSC_INCLUDE_DIRS} CACHE INTERNAL "")
+
+	set(PCMS_PETSC_LIBRARIES ${PETSC_LIBRARY} CACHE INTERNAL "")
+
+endif()
 
 add_subdirectory(src)
 

src/pcms/interpolator/CMakeLists.txt

@@ -13,6 +13,8 @@ set(PCMS_FIELD_TRANSFER_HEADERS
     spline_interpolator.hpp)
 
 set(PCMS_FIELD_TRANSFER_SOURCES mls_interpolation.cpp)
+
+add_subdirectory(mesh_intersection)
 install(FILES ${PCMS_FIELD_TRANSFER_HEADERS}
         DESTINATION include/pcms/interpolator)
 

src/pcms/interpolator/adj_search.hpp

@@ -4,6 +4,7 @@
 #include <pcms/point_search.h>
 
 #include "queue_visited.hpp"
+#include <Kokkos_MathematicalFunctions.hpp>
 
 static constexpr int max_dim = 3;
 
@@ -119,16 +120,20 @@ inline void FindSupports::adjBasedSearch(
 
   pcms::GridPointSearch search_cell(source_mesh, 10, 10);
   auto results = search_cell(target_points);
-  checkTargetPoints(results);
+  // checkTargetPoints(results);
 
   Omega_h::parallel_for(
     nvertices_target,
     OMEGA_H_LAMBDA(const Omega_h::LO id) {
       Queue queue;
       Track visited;
       Omega_h::Real cutoffDistance = radii2[id];
-
       Omega_h::LO source_cell_id = results(id).tri_id;
+      // printf("target = %d source cell id : %d\n",id, source_cell_id);
+      //  if the point lies outside the mesh, it returns the negative cell id
+      //  making the negative cell id to positive
+      if (source_cell_id < 0)
+        source_cell_id = Kokkos::abs(source_cell_id);
       OMEGA_H_CHECK_PRINTF(
         source_cell_id >= 0,
         "ERROR: Source cell id not found for target %d (%f,%f)\n", id,

src/pcms/interpolator/mesh_intersection/CMakeLists.txt

@@ -0,0 +1,44 @@
+file(GLOB_RECURSE MESH_INTERSECTION_HEADERS
+    "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp")
+
+file(GLOB_RECURSE MESH_INTERSECTION_SOURCES
+    "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp")
+
+add_library(pcms_mesh_intersection
+    ${MESH_INTERSECTION_HEADERS}
+    ${MESH_INTERSECTION_SOURCES}
+)
+
+target_include_directories(pcms_mesh_intersection PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+    $<INSTALL_INTERFACE:include/pcms/interpolator/mesh_intersection>
+)
+
+target_link_libraries(pcms_mesh_intersection PUBLIC
+    pcms::core
+    Omega_h::omega_h
+)
+
+if (PCMS_ENABLE_MESHFIELDS)
+    target_link_libraries(pcms_mesh_intersection PUBLIC meshfields::meshfields)
+endif()
+
+if (PCMS_ENABLE_PETSC)
+	target_include_directories(pcms_mesh_intersection PUBLIC ${PCMS_PETSC_INCLUDES})
+	target_link_libraries(pcms_mesh_intersection PUBLIC ${PCMS_PETSC_LIBRARIES})
+endif()
+
+target_compile_definitions(pcms_mesh_intersection PUBLIC R3D_USE_CUDA)
+
+install(FILES ${MESH_INTERSECTION_HEADERS}
+        DESTINATION include/pcms/interpolator/mesh_intersection)
+
+install(TARGETS pcms_mesh_intersection
+        EXPORT meshIntersectionTargets
+        INCLUDES DESTINATION include/pcms/interpolator/mesh_intersection
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+
+install(EXPORT meshIntersectionTargets
+        FILE pcms_mesh_intersection-targets.cmake
+        NAMESPACE pcms::
+        DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms)

src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp

@@ -0,0 +1,129 @@
+/**
+ * @file calculate_load_vector.hpp
+ * @brief Routines for assembling global load vector in conservative field
+ * projection.
+ *
+ * Provides functionality to compute and assemble the global load vector
+ * used in Galerkin-based conservative field transfer between non-conforming
+ * meshes.
+ *
+ * @created by Abhiyan Paudel
+ * @date August 2025
+ */
+
+#ifndef PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP
+#define PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP
+#include <Omega_h_adapt.hpp>
+#include <Omega_h_array_ops.hpp>
+#include <Omega_h_atomics.hpp> //Omega_h::atomic_fetch_add
+#include <Omega_h_build.hpp>
+#include <Omega_h_class.hpp>
+#include <Omega_h_compare.hpp>
+#include <Omega_h_dbg.hpp>
+#include <Omega_h_file.hpp> //Omega_h::binary
+#include <Omega_h_for.hpp>
+#include <Omega_h_recover.hpp> //project_by_fit
+#include <Omega_h_shape.hpp>
+#include <Omega_h_timer.hpp>
+#include <iomanip> //precision
+#include <iostream>
+#include <petscvec_kokkos.hpp>
+#include <sstream> //ostringstream
+
+#include <pcms/interpolator/mesh_interscetion/load_vector_integrator.hpp>
+#include <MeshField.hpp>
+
+#include <petscmat.h>
+
+// detect floating point exceptions
+#include <fenv.h>
+
+using ExecutionSpace = Kokkos::DefaultExecutionSpace;
+using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space;
+
+/**
+ * @brief Assembles the global load vector.
+ *
+ * This function computes the unassembled local load vector contributions for
+ * each triangular element in the target mesh using `buildLoadVector()` and then
+ * assembles them into a global PETSc vector in COO format.
+ *
+ *
+ * @param target_mesh The target Omega_h mesh to which the scalar field is being
+ * projected.
+ * @param source_mesh The source Omega_h mesh containing the original scalar
+ * field values.
+ * @param intersection Precomputed intersection data for each target element.
+ *                     Includes the number and indices of intersecting source
+ * elements.
+ * @param source_values Nodal scalar field values defined on the source mesh.
+ * @param[out] loadVec_out Pointer to a PETSc Vec where the assembled load
+ * vector will be stored.
+ *
+ * @return PetscErrorCode Returns PETSC_SUCCESS if successful, or an appropriate
+ * PETSc error code otherwise.
+ *
+ * @note
+ * - Works for 2D linear triangular elements.
+ * - Uses COO-style preallocation and insertion into the PETSc vector.
+ * - Internally calls `buildLoadVector()` to compute per-element contributions.
+ * - The resulting vector is used as the right-hand side (RHS) in a projection
+ * solve.
+ *
+ * @see buildLoadVector,IntersectionResults
+ */
+
+namespace pcms
+{
+
+inline PetscErrorCode calculateLoadVector(
+  Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh,
+  const IntersectionResults& intersection, const Omega_h::Reals& source_values,
+  Vec* loadVec_out)
+{
+
+  PetscFunctionBeginUser;
+  const int numNodesPerTri = 3;
+
+  const int nnz = target_mesh.nelems() * numNodesPerTri;
+
+  // Allocate COO indices and values
+  PetscInt* coo_i;
+  PetscScalar* coo_vals;
+  PetscCall(PetscMalloc2(nnz, &coo_i, nnz, &coo_vals));
+
+  // Fill COO global indices and values
+  auto elmVerts = Omega_h::HostRead(target_mesh.ask_elem_verts());
+  auto elmLoadVector =
+    buildLoadVector(target_mesh, source_mesh, intersection, source_values);
+
+  auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector);
+  Kokkos::deep_copy(hostElmLoadVector, elmLoadVector);
+
+  PetscInt idx = 0;
+  for (PetscInt e = 0; e < target_mesh.nelems(); ++e) {
+    for (PetscInt vi = 0; vi < numNodesPerTri; ++vi) {
+      coo_i[idx] = elmVerts[numNodesPerTri * e + vi];
+      coo_vals[idx] = hostElmLoadVector(numNodesPerTri * e + vi);
+      ++idx;
+    }
+  }
+
+  // create vector with preallocated COO structure
+  Vec vec;
+  PetscCall(VecCreate(PETSC_COMM_WORLD, &vec));
+  PetscCall(VecSetSizes(vec, target_mesh.nverts(), PETSC_DECIDE));
+  PetscCall(VecSetFromOptions(vec));
+  PetscCall(VecSetPreallocationCOO(vec, nnz, coo_i));
+  PetscCall(VecSetValuesCOO(vec, coo_vals, ADD_VALUES));
+  PetscCall(PetscFree2(coo_i, coo_vals));
+
+  if (target_mesh.nelems() < 10) {
+    PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_WORLD));
+  }
+
+  *loadVec_out = vec;
+  PetscFunctionReturn(PETSC_SUCCESS);
+}
+} // namespace pcms
+#endif