Conservative projection in C-Fortran API

src/pcms/capi/interpolator.cpp

@@ -133,3 +133,156 @@ void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input,
 
   mls_interpolator->eval(input_array, output_array);
 }
+
+// ---------------------------------------------------------------------------
+// Conservative Projection (mesh-intersection-based)
+// ---------------------------------------------------------------------------
+
+#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS)
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <pcms/field/field.h>
+#include <pcms/field/function_space/lagrange.h>
+#include <pcms/transfer/omega_h_conservative_projection.hpp>
+#include <pcms/utility/arrays.h>
+#include <cstring>
+#include <memory>
+#include <string>
+
+namespace
+{
+
+// Trims trailing whitespace from a C string (Fortran strings may be padded).
+std::string trim_filename(const char* filename)
+{
+  auto fname = std::string(filename);
+  fname.erase(fname.find_last_not_of(" \n\r\t") + 1);
+  return fname;
+}
+
+struct ConservativeProjectionContext
+{
+  // Library must be declared before meshes (Omega_h::Mesh holds a pointer to
+  // it).
+  Omega_h::Library library;
+  Omega_h::Mesh source_mesh;
+  Omega_h::Mesh target_mesh;
+  std::shared_ptr<pcms::LagrangeFunctionSpace> source_space;
+  std::shared_ptr<pcms::LagrangeFunctionSpace> target_space;
+  std::unique_ptr<pcms::OmegaHConservativeProjection> projection;
+  pcms::Function<pcms::Real> source_field;
+  pcms::Function<pcms::Real> target_field;
+
+  ConservativeProjectionContext(const char* source_mesh_name,
+                                const char* target_mesh_name, int src_order,
+                                int tgt_order)
+    : library(nullptr, nullptr, MPI_COMM_SELF),
+      source_mesh(Omega_h::binary::read(trim_filename(source_mesh_name),
+                                        library.world())),
+      target_mesh(Omega_h::binary::read(trim_filename(target_mesh_name),
+                                        library.world())),
+      source_space(pcms::LagrangeFunctionSpace::FromMesh(
+        source_mesh, src_order, 1, pcms::CoordinateSystem::Cartesian, "global",
+        pcms::LagrangeFunctionSpace::Backend::OmegaH)),
+      target_space(pcms::LagrangeFunctionSpace::FromMesh(
+        target_mesh, tgt_order, 1, pcms::CoordinateSystem::Cartesian, "global",
+        pcms::LagrangeFunctionSpace::Backend::OmegaH)),
+      projection(std::make_unique<pcms::OmegaHConservativeProjection>(
+        *source_space, *target_space)),
+      source_field(source_space->CreateFunction<pcms::Real>()),
+      target_field(target_space->CreateFunction<pcms::Real>())
+  {
+  }
+};
+
+} // namespace
+
+PcmsConservativeProjectionHandle pcms_create_conservative_projection(
+  const char* source_mesh_name, int source_order, const char* target_mesh_name,
+  int target_order)
+{
+  auto* ctx = new ConservativeProjectionContext(
+    source_mesh_name, target_mesh_name, source_order, target_order);
+  return {reinterpret_cast<void*>(ctx)};
+}
+
+int pcms_conservative_projection_get_source_size(
+  PcmsConservativeProjectionHandle projection)
+{
+  auto* ctx =
+    reinterpret_cast<ConservativeProjectionContext*>(projection.pointer);
+  return ctx->source_space->GetLayout()->GetNumOwnedDofHolder();
+}
+
+int pcms_conservative_projection_get_target_size(
+  PcmsConservativeProjectionHandle projection)
+{
+  auto* ctx =
+    reinterpret_cast<ConservativeProjectionContext*>(projection.pointer);
+  return ctx->target_space->GetLayout()->GetNumOwnedDofHolder();
+}
+
+void pcms_conservative_projection_apply(
+  PcmsConservativeProjectionHandle projection, void* source_data,
+  int source_size, void* target_data, int target_size)
+{
+  auto* ctx =
+    reinterpret_cast<ConservativeProjectionContext*>(projection.pointer);
+
+  // Copy source data into internal field
+  auto source_view = pcms::Rank2View<const pcms::Real, pcms::HostMemorySpace>(
+    reinterpret_cast<pcms::Real*>(source_data), source_size, 1);
+  ctx->source_field.SetDOFHolderDataHost(source_view);
+
+  // Apply conservative projection
+  ctx->projection->Apply(ctx->source_field, ctx->target_field);
+
+  // Copy result back to user buffer
+  auto target_view = ctx->target_field.GetDOFHolderDataHost();
+  auto flat = pcms::FlattenToRank1View(target_view);
+  std::memcpy(target_data, flat.data_handle(),
+              static_cast<std::size_t>(target_size) * sizeof(double));
+}
+
+void pcms_destroy_conservative_projection(
+  PcmsConservativeProjectionHandle projection)
+{
+  if (projection.pointer != nullptr) {
+    delete reinterpret_cast<ConservativeProjectionContext*>(projection.pointer);
+  }
+}
+
+#else // !(PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS)
+
+// Stub implementations when PETSc or MeshFields are unavailable
+
+PcmsConservativeProjectionHandle pcms_create_conservative_projection(
+  const char*, int, const char*, int)
+{
+  pcms::printError("Conservative projection requires PCMS_ENABLE_PETSC and "
+                   "PCMS_ENABLE_MESHFIELDS\n");
+  return {nullptr};
+}
+
+int pcms_conservative_projection_get_source_size(
+  PcmsConservativeProjectionHandle)
+{
+  return 0;
+}
+
+int pcms_conservative_projection_get_target_size(
+  PcmsConservativeProjectionHandle)
+{
+  return 0;
+}
+
+void pcms_conservative_projection_apply(PcmsConservativeProjectionHandle, void*,
+                                        int, void*, int)
+{
+  pcms::printError("Conservative projection requires PCMS_ENABLE_PETSC and "
+                   "PCMS_ENABLE_MESHFIELDS\n");
+}
+
+void pcms_destroy_conservative_projection(PcmsConservativeProjectionHandle) {}
+
+#endif

src/pcms/capi/interpolator.h

@@ -175,6 +175,85 @@ void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator);
 void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input,
                       int input_size, void* output, int output_size);
 
+/**
+ * @brief Holds a void pointer to ConservativeProjectionContext
+ */
+struct PcmsConservativeProjectionHandle
+{
+  void* pointer;
+};
+
+/**
+ * @brief Typedef for PcmsConservativeProjectionHandle struct
+ * @copydetails PcmsConservativeProjectionHandle
+ */
+typedef struct PcmsConservativeProjectionHandle
+  PcmsConservativeProjectionHandle;
+
+/**
+ * @brief Create a mesh-intersection-based conservative projection
+ * @param source_mesh_name C-string path to source Omega_h mesh file (.osh)
+ * @param source_order Order of source Lagrange space (0 or 1)
+ * @param target_mesh_name C-string path to target Omega_h mesh file (.osh)
+ * @param target_order Order of target Lagrange space (0 or 1)
+ * @return Handle to the created conservative projection
+ *
+ * @details Loads both meshes internally, creates LagrangeFunctionSpace objects
+ * (Backend::OmegaH), performs mesh intersection via intersectTargets(), builds
+ * quadrature data, assembles and factors the target mass matrix (PETSc KSP),
+ * and caches everything for repeated Apply calls. The meshes are owned by
+ * the returned handle and live until pcms_destroy_conservative_projection.
+ *
+ * @note Requires PCMS_ENABLE_PETSC and PCMS_ENABLE_MESHFIELDS. Returns a
+ * handle with null pointer if either is unavailable.
+ * @note Only scalar fields (1 component) are supported.
+ * @note Source and target orders can differ (e.g., P0<->P1).
+ */
+PcmsConservativeProjectionHandle pcms_create_conservative_projection(
+  const char* source_mesh_name, int source_order, const char* target_mesh_name,
+  int target_order);
+
+/**
+ * @brief Get the number of source DOF holders
+ * @param projection Handle to the conservative projection
+ * @return Number of source DOF holders
+ */
+int pcms_conservative_projection_get_source_size(
+  PcmsConservativeProjectionHandle projection);
+
+/**
+ * @brief Get the number of target DOF holders
+ * @param projection Handle to the conservative projection
+ * @return Number of target DOF holders
+ */
+int pcms_conservative_projection_get_target_size(
+  PcmsConservativeProjectionHandle projection);
+
+/**
+ * @brief Apply conservative projection
+ * @param projection Handle to the projection
+ * @param source_data Flat array of source field values (size = source_size)
+ * @param source_size Number of source DOF holders
+ * @param target_data Flat array to receive target field values (size =
+ * target_size)
+ * @param target_size Number of target DOF holders
+ *
+ * @details Copies source data into the internal field, calls
+ * OmegaHConservativeProjection::Apply(), and copies the result back.
+ * This is the cheap per-time-step call: no mesh intersection,
+ * no matrix assembly, no refactorization.
+ */
+void pcms_conservative_projection_apply(
+  PcmsConservativeProjectionHandle projection, void* source_data,
+  int source_size, void* target_data, int target_size);
+
+/**
+ * @brief Destroy conservative projection and free all internal resources
+ * @param projection Handle to the conservative projection to destroy
+ */
+void pcms_destroy_conservative_projection(
+  PcmsConservativeProjectionHandle projection);
+
 #ifdef __cplusplus
 }
 #endif

src/pcms/fortranapi/interpolator.i

@@ -29,3 +29,25 @@ void pcms_kokkos_initialize_without_args();
 void pcms_kokkos_finalize();
 
 void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size);
+
+
+// --- Conservative Projection (mesh-intersection-based) ---
+
+struct PcmsConservativeProjectionHandle {
+  void* pointer;
+};
+typedef struct PcmsConservativeProjectionHandle PcmsConservativeProjectionHandle;
+
+PcmsConservativeProjectionHandle pcms_create_conservative_projection(
+    const char* source_mesh_name, int source_order,
+    const char* target_mesh_name, int target_order);
+int pcms_conservative_projection_get_source_size(
+    PcmsConservativeProjectionHandle projection);
+int pcms_conservative_projection_get_target_size(
+    PcmsConservativeProjectionHandle projection);
+void pcms_conservative_projection_apply(
+    PcmsConservativeProjectionHandle projection,
+    void* source_data, int source_size,
+    void* target_data, int target_size);
+void pcms_destroy_conservative_projection(
+    PcmsConservativeProjectionHandle projection);