From 96012d63a3744fb61fac68b9fb66138b0a53dd6e Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Sun, 17 Aug 2025 21:16:37 -0400 Subject: [PATCH 01/29] add new rbf functions --- src/pcms/interpolator/mls_interpolation.cpp | 122 +++++++++++++++--- .../interpolator/mls_interpolation_impl.hpp | 37 +++--- 2 files changed, 128 insertions(+), 31 deletions(-) diff --git a/src/pcms/interpolator/mls_interpolation.cpp b/src/pcms/interpolator/mls_interpolation.cpp index 00140f8e..d66e24b1 100644 --- a/src/pcms/interpolator/mls_interpolation.cpp +++ b/src/pcms/interpolator/mls_interpolation.cpp @@ -31,11 +31,7 @@ struct RBF_GAUSSIAN r_sq); double r = Kokkos::sqrt(r_sq); - double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho; - double limit = 1 - ratio; - - if (limit < 0) { + if (rho_sq < r_sq) { phi = 0; } else { @@ -55,13 +51,14 @@ struct RBF_C4 { double phi; double r = Kokkos::sqrt(r_sq); - OMEGA_H_CHECK_PRINTF( - rho_sq > 0, "ERROR: rho_sq in rbf has to be positive, but got %.16f\n", - rho_sq); double rho = Kokkos::sqrt(rho_sq); double ratio = r / rho; double limit = 1 - ratio; - if (limit < 0) { + + OMEGA_H_CHECK_PRINTF( + rho_sq > 0, "ERROR: rho_sq in rbf has to be positive, but got %.16f\n", + rho_sq); + if (rho_sq < r_sq) { phi = 0; } else { @@ -86,15 +83,12 @@ struct RBF_CONST double operator()(double r_sq, double rho_sq) const { double phi; - double r = Kokkos::sqrt(r_sq); OMEGA_H_CHECK_PRINTF( rho_sq > 0, "ERROR: rho_sq in rbf has to be positive, but got %.16f\n", rho_sq); - double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho; - double limit = 1 - ratio; - if (limit < 0) { - phi = 0; + + if (rho_sq < r_sq) { + phi = 0.0; } else { phi = 1.0; @@ -113,6 +107,104 @@ struct NoOp double operator()(double, double) const { return 1.0; } }; +struct RBF_MULTIQUADRIC +{ + OMEGA_H_INLINE + double operator()(double r_sq, double rho_sq) const + { + double phi; + double r = Kokkos::sqrt(r_sq); + double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho; + if (rho_sq < r_sq) { + phi = 0.0; + } else { + phi = Kokkos::sqrt(1.0 + ratio * ratio); + } + + OMEGA_H_CHECK_PRINTF(!std::isnan(phi), + "ERROR: phi in rbf is NaN. r_sq, rho_sq = (%f, %f)\n", + r_sq, rho_sq); + + return phi; + } +}; + +struct RBF_INVMULTIQUADRIC +{ + OMEGA_H_INLINE + double operator()(double r_sq, double rho_sq) const + { + double phi; + double r = Kokkos::sqrt(r_sq); + double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho if (rho_sq < r_sq) + { + phi = 0.0; + } + else + { + phi = 1.0 / Kokkos::sqrt(1.0 + ratio * ratio); + } + + OMEGA_H_CHECK_PRINTF(!std::isnan(phi), + "ERROR: phi in rbf is NaN. r_sq, rho_sq = (%f, %f)\n", + r_sq, rho_sq); + + return phi; + } +}; + +struct RBF_THINPLATESPLINE +{ + OMEGA_H_INLINE + double operator()(double r_sq, double rho_sq) const + { + double phi; + double r = Kokkos::sqrt(r_sq); + double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho if (rho_sq < r_sq) + { + phi = 0.0; + } + else + { + phi = r * r * Kokkos::log(ratio); + } + + OMEGA_H_CHECK_PRINTF(!std::isnan(phi), + "ERROR: phi in rbf is NaN. r_sq, rho_sq = (%f, %f)\n", + r_sq, rho_sq); + + return phi; + } +}; + +struct RBF_CUBIC +{ + OMEGA_H_INLINE + double operator()(double r_sq, double rho_sq) const + { + double phi; + double r = Kokkos::sqrt(r_sq); + double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho double limit = 1 - ratio if (rho_sq < r_sq) + { + phi = 0.0; + } + else + { + phi = limit * limit * limit; + } + + OMEGA_H_CHECK_PRINTF(!std::isnan(phi), + "ERROR: phi in rbf is NaN. r_sq, rho_sq = (%f, %f)\n", + r_sq, rho_sq); + + return phi; + } +}; + Omega_h::Write mls_interpolation( const Omega_h::Reals source_values, const Omega_h::Reals source_coordinates, const Omega_h::Reals target_coordinates, const SupportResults& support, diff --git a/src/pcms/interpolator/mls_interpolation_impl.hpp b/src/pcms/interpolator/mls_interpolation_impl.hpp index 499f7a71..abab764c 100644 --- a/src/pcms/interpolator/mls_interpolation_impl.hpp +++ b/src/pcms/interpolator/mls_interpolation_impl.hpp @@ -32,7 +32,9 @@ static constexpr int MAX_DIM = 6; * eval_basis_vector are needed to evaluate the polynomial basis for any degree * and dimension For instance, polynomial basis vector for dim = 2 and degree = * 3 at the point (x,y) looks like {1, x, y, xx, xy, yy, xxx, xxy, xyy,yyy}. The - * slices can be written as [1] degree 0 [x] & [y] degree 1 + * slices can be written as + * [1] degree 0 + * [x] & [y] degree 1 * [xx] & [xy, yy] degree 2 * [xxx] & [xxy,xyy, yyy] degree 3 * @@ -55,6 +57,7 @@ namespace pcms namespace detail { + /** * @brief Computes the slice lengths of the polynomial basis * @@ -463,21 +466,25 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, calculate_scratch_shared_size(support, ntargets, basis_size, dim); team_policy tp(ntargets, Kokkos::AUTO); + // printf("num of targets: ", ntargets); int scratch_size = tp.scratch_size_max(1); // printf("Scratch Size = %d\n", scratch_size); // printf("Shared Size = %d\n", shared_size); + // printf("Shared Size to be set = %d bytes (%.2f KB)\n", shared_size, + // shared_size / 1024.0); PCMS_ALWAYS_ASSERT(scratch_size > shared_size); // calculates the interpolated values Kokkos::parallel_for( - "MLS coefficients", tp.set_scratch_size(1, Kokkos::PerTeam(scratch_size)), + "MLS coefficients", tp.set_scratch_size(1, Kokkos::PerTeam(shared_size)), KOKKOS_LAMBDA(const member_type& team) { int league_rank = team.league_rank(); int start_ptr = support.supports_ptr[league_rank]; int end_ptr = support.supports_ptr[league_rank + 1]; int nsupports = end_ptr - start_ptr; + // Logger logger(0); // local_source_point stores the coordinates of source supports of a // given target ScratchMatView local_source_points(team.team_scratch(1), nsupports, dim); @@ -485,6 +492,7 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, // rbf function values of source supports Phi(n,n) ScratchVecView phi_vector(team.team_scratch(1), nsupports); + // rbf function values of source supports Phi(n,n) // vondermonde matrix P from the vectors of basis vector of supports ScratchMatView vandermonde_matrix(team.team_scratch(1), nsupports, basis_size); @@ -505,8 +513,6 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, fill(0.0, team, target_basis_vector); fill(0.0, team, solution_coefficients); - // Logger logger(15); - /** * * the local_source_points is of the type ScratchMatView with @@ -518,22 +524,21 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, for (int j = start_ptr; j < end_ptr; ++j) { count++; auto index = support.supports_idx[j]; - for (int i = 0; i < dim; ++i) { local_source_points(count, i) = source_coordinates[index * dim + i]; } } - // logger.logMatrix(team, LogLevel::DEBUG, local_source_points, - // "Support Coordinates"); + // logger.logMatrix(team, LogLevel::DEBUG, local_source_points, + // "Support Coordinates"); double target_point[MAX_DIM] = {}; for (int i = 0; i < dim; ++i) { target_point[i] = target_coordinates[league_rank * dim + i]; } - // logger.logArray(team, LogLevel::DEBUG, target_point, dim, - // "Target points"); + // logger.logArray(team, LogLevel::DEBUG, target_point, dim, + // "Target points"); /** phi(nsupports) is the array of rbf functions evaluated at the * source supports In the actual implementation, Phi(nsupports, @@ -567,9 +572,9 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, team.team_barrier(); - // logger.log(team, LogLevel::DEBUG, "The search starts"); - // logger.logVector(team, LogLevel::DEBUG, support_values, "Support - // values"); + // logger.log(team, LogLevel::DEBUG, "The search starts"); + // logger.logVector(team, LogLevel::DEBUG, support_values, "Support + // values"); /** * step 3: normalize local source supports and target point @@ -600,8 +605,8 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, team.team_barrier(); - // logger.logMatrix(team, LogLevel::DEBUG, vandermonde_matrix, - // "vandermonde matrix"); + // logger.logMatrix(team, LogLevel::DEBUG, vandermonde_matrix, + // "vandermonde matrix"); OMEGA_H_CHECK_PRINTF( @@ -622,8 +627,8 @@ void mls_interpolation(RealConstDefaultScalarArrayView source_values, // printf("Target Point : %d \t\t Value: %5.6f\n", league_rank, // target_value); - // logger.logScalar(team, LogLevel::DEBUG, target_value, - // "interpolated value"); + // logger.logScalar(team, LogLevel::DEBUG, target_value, + // "interpolated value"); if (team.team_rank() == 0) { OMEGA_H_CHECK_PRINTF(!std::isnan(target_value), "Nan at %d\n", league_rank); From 3a149369f37a6e08f080426bdac6603b4a43df68 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 9 Oct 2025 19:22:19 -0400 Subject: [PATCH 02/29] minor changes --- src/pcms/interpolator/adj_search.hpp | 9 +++- src/pcms/interpolator/mls_interpolation.cpp | 57 +++++++++++---------- 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/src/pcms/interpolator/adj_search.hpp b/src/pcms/interpolator/adj_search.hpp index aca256fb..f7ceb7db 100644 --- a/src/pcms/interpolator/adj_search.hpp +++ b/src/pcms/interpolator/adj_search.hpp @@ -4,6 +4,7 @@ #include #include "queue_visited.hpp" +#include static constexpr int max_dim = 3; @@ -121,7 +122,7 @@ 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, @@ -129,8 +130,12 @@ inline void FindSupports::adjBasedSearch( 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, diff --git a/src/pcms/interpolator/mls_interpolation.cpp b/src/pcms/interpolator/mls_interpolation.cpp index 364c9846..de59a67f 100644 --- a/src/pcms/interpolator/mls_interpolation.cpp +++ b/src/pcms/interpolator/mls_interpolation.cpp @@ -5,12 +5,13 @@ namespace pcms { -// RBF_GAUSSIAN Functor +// 'a' is a shape parameter +// the value of 'a' is higher if the data is localized +// the value of 'a' is smaller if the data is farther + +// gaussian struct RBF_GAUSSIAN { - // 'a' is a spreading factor/decay factor - // the value of 'a' is higher if the data is localized - // the value of 'a' is smaller if the data is farther double a; @@ -42,7 +43,7 @@ struct RBF_GAUSSIAN } }; -// RBF_C4 Functor +// c4 struct RBF_C4 { @@ -74,8 +75,7 @@ struct RBF_C4 } }; -// RBF_const Functor -// +// constant/uniform struct RBF_CONST { @@ -101,25 +101,30 @@ struct RBF_CONST } }; +// no operations struct NoOp { OMEGA_H_INLINE double operator()(double, double) const { return 1.0; } }; +// multiquadric struct RBF_MULTIQUADRIC { + + double a; + RBF_MULTIQUADRIC(double a_val) : a(a_val) {} + OMEGA_H_INLINE double operator()(double r_sq, double rho_sq) const { double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho; if (rho_sq < r_sq) { phi = 0.0; } else { - phi = Kokkos::sqrt(1.0 + ratio * ratio); + phi = Kokkos::sqrt(1.0 + a * a * r * r); } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -130,21 +135,22 @@ struct RBF_MULTIQUADRIC } }; +// inverse multiquadric struct RBF_INVMULTIQUADRIC { + double a; + RBF_INVMULTIQUADRIC(double a_val) : a(a_val) {} + OMEGA_H_INLINE double operator()(double r_sq, double rho_sq) const { double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho if (rho_sq < r_sq) - { + if (rho_sq < r_sq) { phi = 0.0; - } - else - { - phi = 1.0 / Kokkos::sqrt(1.0 + ratio * ratio); + } else { + phi = 1.0 / Kokkos::sqrt(1.0 + a * a * r * r); } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -155,6 +161,7 @@ struct RBF_INVMULTIQUADRIC } }; +// thin plate spline struct RBF_THINPLATESPLINE { OMEGA_H_INLINE @@ -163,13 +170,11 @@ struct RBF_THINPLATESPLINE double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho if (rho_sq < r_sq) - { + + if (rho_sq < r_sq) { phi = 0.0; - } - else - { - phi = r * r * Kokkos::log(ratio); + } else { + phi = r * r * Kokkos::log(r); } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -180,6 +185,7 @@ struct RBF_THINPLATESPLINE } }; +// cubic struct RBF_CUBIC { OMEGA_H_INLINE @@ -188,13 +194,10 @@ struct RBF_CUBIC double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); - double ratio = r / rho double limit = 1 - ratio if (rho_sq < r_sq) - { + if (rho_sq < r_sq) { phi = 0.0; - } - else - { - phi = limit * limit * limit; + } else { + phi = r * r * r; } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), From 9e925287cd458d87d920a00abad9ce839ec698a8 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 12:42:24 -0400 Subject: [PATCH 03/29] add doxygen comment --- src/pcms/interpolator/pcms_interpolator_view_utils.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/pcms/interpolator/pcms_interpolator_view_utils.hpp b/src/pcms/interpolator/pcms_interpolator_view_utils.hpp index 1c52752f..611b3f2c 100644 --- a/src/pcms/interpolator/pcms_interpolator_view_utils.hpp +++ b/src/pcms/interpolator/pcms_interpolator_view_utils.hpp @@ -49,6 +49,13 @@ void fill(double value, member_type team, ScratchVecView vector) [=](int j) { vector(j) = value; }); } +/** + * @brief Evaluates the square root of each element in scratch view + * + * @param team The team member + * @pram array The scratch vector + * + */ KOKKOS_INLINE_FUNCTION void find_sq_root_each(member_type team, ScratchVecView& array) { From 219c35e608324543e361318c816e8ca81d960412 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 12:53:33 -0400 Subject: [PATCH 04/29] add mesh intersection files --- .../mesh_intersection/calculateLoadVector.hpp | 113 ++++++ .../mesh_intersection/calculateMassMatrix.hpp | 167 +++++++++ .../loadVectorIntegrator.hpp | 341 ++++++++++++++++++ .../massMatrixIntegrator.hpp | 81 +++++ 4 files changed, 702 insertions(+) create mode 100644 src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp create mode 100644 src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp create mode 100644 src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp create mode 100644 src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp diff --git a/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp b/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp new file mode 100644 index 00000000..6aa2219d --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp @@ -0,0 +1,113 @@ +/** + * @file calculateLoadVector.hpp + * @brief Functions for calculating load vector on finite element meshes + * + * This file contains functions for creating and computing load vectors + * for finite element calculations using Omega_h mesh structures and PETSc. + */ + +#ifndef COMPUTING_AT_SCALE_CALCULATE_LOAD_VECTOR_HPP +#define COMPUTING_AT_SCALE_CALCULATE_LOAD_VECTOR_HPP +#include +#include +#include //Omega_h::atomic_fetch_add +#include +#include +#include +#include +#include //Omega_h::binary +#include +#include //project_by_fit +#include +#include +#include //precision +#include +#include +#include //ostringstream + +#include "loadVectorIntegrator.hpp" +#include + +#include + +// detect floating point exceptions +#include + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; + +/** + * @brief Calculates the load + * + * This function constructs a mass matrix based on the provided mesh using + * a finite element approach. It creates coordinate field elements, builds + * the mass matrix using the massMatrixIntegrator, and sets up the PETSc matrix + * with appropriate values. + * + * @param mesh The Omega_h mesh to calculate the mass matrix for + * @param[out] mass_out Pointer to the resulting mass matrix + * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) + */ +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); + + /* + std::cerr << " DEBUG: printing the load vector \n"; + for (int i = 0; i < hostElmLoadVector.size(); ++i) { + std::cout << hostElmLoadVector[i] << " "; + } + std::cout << "\n"; + + target_mesh.add_tag(2, "elmLoadVector", 3, + Omega_h::read(Omega_h::Write(elmLoadVector))); + + Omega_h::vtk::write_parallel("loadVector.vtk", &target_mesh, 2); + */ + 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); +} + +#endif diff --git a/src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp b/src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp new file mode 100644 index 00000000..71588503 --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp @@ -0,0 +1,167 @@ +/** + * @file calculateMassMatrix.hpp + * @brief Functions for calculating mass matrices on finite element meshes + * @author [Cameron Smith] + * @date April, 2025 + * + * This file contains functions for creating and computing mass matrices + * for finite element calculations using Omega_h mesh structures and PETSc. + */ + +#ifndef PCMS_INTERPOLATOR_CALCULATE_MASS_MATRIX_HPP +#define PCMS_INTERPOLATOR_CALCULATE_MASS_MATRIX_HPP + +#include +#include +#include //Omega_h::atomic_fetch_add +#include +#include +#include +#include +#include //Omega_h::binary +#include +#include //project_by_fit +#include +#include +#include //precision +#include +#include +#include //ostringstream + +#include "massMatrixIntegrator.hpp" +#include + +#include + +// detect floating point exceptions +#include + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; + +/** + * @brief Sets a constant value for a field at all mesh vertices + * + * @tparam ShapeField Type of the shape field + * @param mesh The Omega_h mesh + * @param field The field to set values for + * @param val The constant value to set at all vertices + */ +template +void setFieldAtVertices(Omega_h::Mesh& mesh, ShapeField field, + const MeshField::Real val) +{ + const auto MeshDim = mesh.dim(); + auto setFieldAtVertices = KOKKOS_LAMBDA(const int& vtx) + { + field(0, 0, vtx, MeshField::Vertex) = val; + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, + setFieldAtVertices, "setFieldAtVertices"); +} + +/** + * @brief Creates a PETSc matrix based on mesh connectivity + * + * This function creates a sparse matrix with the proper sparsity pattern + * according to the mesh connectivity. The matrix size corresponds to the + * number of vertices in the mesh. + * + * @param mesh The Omega_h mesh to create the matrix from + * @param[out] A Pointer to the PETSc matrix to be created + * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) + */ +static PetscErrorCode CreateMatrix(Omega_h::Mesh& mesh, Mat* A) +{ + const auto numNodesPerTri = 3; // FIXME query the mesh + const auto matSize = numNodesPerTri * numNodesPerTri * mesh.nelems(); + auto elmVerts = Omega_h::HostRead(mesh.ask_elem_verts()); + PetscInt *oor, *ooc, cnt = 0; + PetscFunctionBeginUser; + PetscCall(MatCreate(PETSC_COMM_WORLD, A)); + PetscCall( + MatSetSizes(*A, mesh.nverts(), mesh.nverts(), PETSC_DECIDE, PETSC_DECIDE)); + PetscCall(MatSetFromOptions(*A)); + /* determine for each entry in each element stiffness matrix the global row + * and column */ + /* since the element is triangular with piecewise linear basis functions there + * are three degrees of freedom per element, one for each vertex */ + PetscCall(PetscMalloc2(matSize, &oor, matSize, &ooc)); + for (PetscInt e = 0; e < mesh.nelems(); e++) { + for (PetscInt vi = 0; vi < numNodesPerTri; vi++) { + for (PetscInt vj = 0; vj < numNodesPerTri; vj++) { + oor[cnt] = elmVerts[numNodesPerTri * e + vi]; + ooc[cnt++] = elmVerts[numNodesPerTri * e + vj]; + } + } + } + PetscCall(MatSetPreallocationCOO(*A, matSize, oor, ooc)); + PetscCall(PetscFree2(oor, ooc)); + PetscFunctionReturn(PETSC_SUCCESS); +} + +/** + * @brief Calculates the mass matrix for a given mesh + * + * This function constructs a mass matrix based on the provided mesh using + * a finite element approach. It creates coordinate field elements, builds + * the mass matrix using the massMatrixIntegrator, and sets up the PETSc matrix + * with appropriate values. + * + * @param mesh The Omega_h mesh to calculate the mass matrix for + * @param[out] mass_out Pointer to the resulting mass matrix + * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) + */ +inline PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out) +{ + PetscFunctionBeginUser; + + MeshField::OmegahMeshField omf( + mesh); + + const auto ShapeOrder = 1; + auto coordField = omf.getCoordField(); + const auto [shp, map] = + MeshField::Omegah::getTriangleElement(mesh); + MeshField::FieldElement coordFe(mesh.nelems(), coordField, shp, map); + + auto elmMassMatrix = buildMassMatrix(mesh, coordFe); + + auto host_elmMassMatrix = Kokkos::create_mirror_view(elmMassMatrix); + + /* + std::cerr << " DEBUG: printing the element mass matrix \n"; + for (int i = 0; i < host_elmMassMatrix.size(); ++i) { + std::cout << host_elmMassMatrix[i] << " "; + } + + */ + + // std::cout << "\n"; + // mesh.add_tag(2, "elmMassMatrix", 3 * 3, + // Omega_h::read(Omega_h::Write(elmMassMatrix))); + + // Omega_h::vtk::write_parallel("massMatrix.vtk", &mesh, 2); + + Mat mass; + PetscCall(CreateMatrix(mesh, &mass)); + PetscBool is_kokkos; + PetscCall( + PetscObjectBaseTypeCompare((PetscObject)mass, MATSEQAIJKOKKOS, &is_kokkos)); + PetscCall(MatZeroEntries(mass)); + PetscCall( + MatSetValuesCOO(mass, elmMassMatrix.data(), + INSERT_VALUES)); // FIXME fails here on gpu, calls into host + // implementation... AFAIK, petsc checks + // the type of the input array of values to + // decide which backend to use... + // + if (mesh.nelems() < 10) { + PetscCall(MatView(mass, PETSC_VIEWER_STDOUT_WORLD)); + } + + *mass_out = mass; + PetscFunctionReturn(PETSC_SUCCESS); +} + +#endif diff --git a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp new file mode 100644 index 00000000..69e90615 --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp @@ -0,0 +1,341 @@ +#ifndef PCMS_INTERPOLATOR_LOAD_VECTOR_INTEGRATOR_HPP +#define PCMS_INTERPOLATOR_LOAD_VECTOR_INTEGRATOR_HPP + +#include +#include +#include +#include +#include "mesh_intersection.hpp" +#include + +// computes the load vector for each element +// This routine is only applicable linear elements +// It calculates the rhs of the conservative field transfer formula from source +// to the target mesh + +static const double EPS_AREA = 1e-16; +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> global_from_barycentric( + const MeshField::Vector3& barycentric_coord, + const Omega_h::Few, 3>& verts_coord) +{ + Omega_h::Vector<2> real_coords = {0.0, 0.0}; + + for (int i = 0; i < 3; ++i) { + real_coords[0] += barycentric_coord[i] * verts_coord[i][0]; + real_coords[1] += barycentric_coord[i] * verts_coord[i][1]; + } + return real_coords; +} + +[[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<3> evaluate_barycentric( + const Omega_h::Vector<2>& point, + const r3d::Few, 3>& verts_coord) +{ + Omega_h::Few, 3> omegah_vector; + for (int i = 0; i < 3; ++i) { + omegah_vector[i][0] = verts_coord[i][0]; + omegah_vector[i][1] = verts_coord[i][1]; + } + + auto barycentric_coordinate = + Omega_h::barycentric_from_global<2, 2>(point, omegah_vector); + + return barycentric_coordinate; +} + +[[nodiscard]] OMEGA_H_INLINE double evaluate_function_value( + const Omega_h::Reals& nodal_values, const Omega_h::LOs& faces2nodes, + const Omega_h::Vector<3>& bary_coords, const int elm_id) +{ + + double value = 0; + const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, elm_id); + for (int i = 0; i < 3; ++i) { + int nid = elm_verts[i]; + value += nodal_values[nid] * bary_coords[i]; + } + + return value; +} + +template +struct IntegrationData +{ + Kokkos::View bary_coords; + Kokkos::View weights; + + IntegrationData() + { + auto ip_vec = MeshField::getIntegrationPoints(MeshField::Triangle, order); + std::size_t num_ip = ip_vec.size(); + + bary_coords = Kokkos::View("bary_coords", num_ip); + weights = Kokkos::View("weights", num_ip); + + auto bary_coords_host = Kokkos::create_mirror_view(bary_coords); + auto weights_host = Kokkos::create_mirror_view(weights); + + for (std::size_t i = 0; i < num_ip; ++i) { + bary_coords_host(i) = ip_vec[i].param; + weights_host(i) = ip_vec[i].weight; + } + + Kokkos::deep_copy(bary_coords, bary_coords_host); + Kokkos::deep_copy(weights, weights_host); + } + + int size() const { return bary_coords.extent(0); } +}; + +[[nodiscard]] OMEGA_H_INLINE int get_polygon_cycle_ccw( + const r3d::Polytope<2>& poly, int* order) +{ + const int m = poly.nverts; + if (m <= 0) + return 0; + order[0] = 0; + int prev = -1; + int curr = 0; + for (int s = 1; s < m; ++s) { + const int a = poly.verts[curr].pnbrs[0]; + const int b = poly.verts[curr].pnbrs[1]; + const int next = (a != prev ? a : b); // pick neighbor that's not 'prev' + order[s] = next; + prev = curr; + curr = next; + } + return m; +} + +Kokkos::View buildLoadVector( + Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, const Omega_h::Reals& source_values) +{ + + const auto& tgt_coords = target_mesh.coords(); + const auto& src_coords = source_mesh.coords(); + const auto& tgt_faces2nodes = + target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_faces2nodes = + source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + + IntegrationData<2> integrationPoints; + int npts = integrationPoints.size(); + + // TODO: Make it generalised; hardcoded for liner 2D + Kokkos::View elmLoadVector("elmLoadVector", + target_mesh.nelems() * 3); + int count = 0; + Kokkos::parallel_for( + "calculate load vector", target_mesh.nelems(), + KOKKOS_LAMBDA(const int& elm) { + auto tgt_elm_vert_coords = + get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); + const int start = intersection.tgt2src_offsets[elm]; + const int end = intersection.tgt2src_offsets[elm + 1]; + // printf("number of intersections : %d\n", end-start); + Omega_h::Vector<3> part_integration = {0.0, 0.0, 0.0}; + + for (int i = start; i < end; ++i) { + const int current_src_elm = intersection.tgt2src_indices[i]; + auto src_elm_vert_coords = + get_vert_coords_of_elem(src_coords, src_faces2nodes, current_src_elm); + r3d::Polytope<2> poly; + r3d::intersect_simplices(poly, tgt_elm_vert_coords, + src_elm_vert_coords); + auto nverts = poly.nverts; + int order[r3d::MaxVerts<2>::value]; + auto m = get_polygon_cycle_ccw(poly, order); + auto poly_area = r3d::measure(poly); + // printf("No of verts : %d\n", nverts); + // decompose if more than 3 vertices + // printf(" area from polytope : %f\n", + //r3d::measure(poly)); + double sum_area = 0; + for (int j = 1; j < nverts - 1; ++j) { + // build triangle from poly.verts[order[0]], poly.verts[order[j]], + // poly.verts[order[j+1]] + auto& p0 = poly.verts[order[0]].pos; + auto& p1 = poly.verts[order[j]].pos; + auto& p2 = poly.verts[order[j + 1]].pos; + + Omega_h::Few, 3> tri_coords; + + tri_coords[0] = {p0[0], p0[1]}; + tri_coords[1] = {p1[0], p1[1]}; + tri_coords[2] = {p2[0], p2[1]}; + + Omega_h::Few, 2> basis; + + basis[0] = tri_coords[1] - tri_coords[0]; + basis[1] = tri_coords[2] - tri_coords[0]; + + Omega_h::Real area = + Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + sum_area += area; + + if (area <= EPS_AREA) + continue; // drops duplicates and colinear/degenerates + // printf(" area = %f\n", area); + + for (int ip = 0; ip < npts; ++ip) { + auto bary = integrationPoints.bary_coords(ip); + auto weight = integrationPoints.weights(ip); + + // convert barycentric to real coords in triangle + auto real_coords = global_from_barycentric(bary, tri_coords); + + // evaluate shape function (barycentric wrt target for linear) + auto shape_fn = + evaluate_barycentric(real_coords, tgt_elm_vert_coords); + + // evaluate function at point (barycentric wrt source for linear) + auto src_bary = + evaluate_barycentric(real_coords, src_elm_vert_coords); + auto fval = evaluate_function_value(source_values, src_faces2nodes, + src_bary, current_src_elm); + + // integration + for (int k = 0; k < 3; ++k) { + part_integration[k] += shape_fn[k] * fval * weight * 2 * area; + } + } + } + // printf(" sum_area = %.16e\n", sum_area); + // printf(" poly_area = %.16e\n", poly_area); + // OMEGA_H_CHECK(nearly_equal(sum_area, poly_area)); + } + + for (int j = 0; j < 3; ++j) { + elmLoadVector(elm * 3 + j) = part_integration[j]; + } + }); + + return elmLoadVector; +} + +struct Errors +{ + double proj_err; + double cons_err; +}; + +Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& target_values, + const Omega_h::Reals& source_values) +{ + + const auto& tgt_coords = target_mesh.coords(); + const auto& src_coords = source_mesh.coords(); + const auto& tgt_faces2nodes = + target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_faces2nodes = + source_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + + IntegrationData<2> integrationPoints; + int npts = integrationPoints.size(); + + constexpr double EPS_DEN = 1e-30; + + Kokkos::View accum("accum", 4); + Kokkos::deep_copy(accum, 0.0); + + Kokkos::parallel_for( + "evaluate relative errors", target_mesh.nelems(), + KOKKOS_LAMBDA(const int& elm) { + auto tgt_elm_vert_coords = + get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); + const int start = intersection.tgt2src_offsets[elm]; + const int end = intersection.tgt2src_offsets[elm + 1]; + + double N2 = 0.0, D2 = 0.0, C = 0.0, QD = 0.0; + + for (int i = start; i < end; ++i) { + const int current_src_elm = intersection.tgt2src_indices[i]; + auto src_elm_vert_coords = + get_vert_coords_of_elem(src_coords, src_faces2nodes, current_src_elm); + r3d::Polytope<2> poly; + r3d::intersect_simplices(poly, tgt_elm_vert_coords, + src_elm_vert_coords); + auto nverts = poly.nverts; + int order[r3d::MaxVerts<2>::value]; + auto m = get_polygon_cycle_ccw(poly, order); + auto poly_area = r3d::measure(poly); + + double sum_area = 0.0; + for (int j = 1; j < nverts - 1; ++j) { + // build triangle from poly.verts[order[0]], poly.verts[order[j]], + // poly.verts[order[j+1]] + auto& p0 = poly.verts[order[0]].pos; + auto& p1 = poly.verts[order[j]].pos; + auto& p2 = poly.verts[order[j + 1]].pos; + + Omega_h::Few, 3> tri_coords; + + tri_coords[0] = {p0[0], p0[1]}; + tri_coords[1] = {p1[0], p1[1]}; + tri_coords[2] = {p2[0], p2[1]}; + + Omega_h::Few, 2> basis; + + basis[0] = tri_coords[1] - tri_coords[0]; + basis[1] = tri_coords[2] - tri_coords[0]; + + Omega_h::Real area = + Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); + sum_area += area; + + const double EPS_AREA = Kokkos::fmax(abs_tol, rel_tol * area); + if (area <= EPS_AREA) + continue; // drops duplicates and colinear/degenerates + // printf(" area = %f\n", area); + + for (int ip = 0; ip < npts; ++ip) { + auto bary = integrationPoints.bary_coords(ip); + auto weight = integrationPoints.weights(ip); + + // convert barycentric to real coords in triangle + auto real_coords = global_from_barycentric(bary, tri_coords); + auto tgt_bary = + evaluate_barycentric(real_coords, tgt_elm_vert_coords); + + // evaluate shape function (barycentric wrt target for linear) + auto tgtVal = evaluate_function_value( + target_values, tgt_faces2nodes, tgt_bary, elm); + + // evaluate function at point (barycentric wrt source for linear) + auto src_bary = + evaluate_barycentric(real_coords, src_elm_vert_coords); + auto srcVal = evaluate_function_value( + source_values, src_faces2nodes, src_bary, current_src_elm); + + // integration + auto diff = srcVal - tgtVal; + auto w = 2 * weight * area; + N2 += diff * diff * w; + D2 += srcVal * srcVal * w; + C += diff * w; + QD += srcVal * w; + } + } + } + + Kokkos::atomic_add(&accum(0), N2); + Kokkos::atomic_add(&accum(1), D2); + Kokkos::atomic_add(&accum(2), C); + Kokkos::atomic_add(&accum(3), QD); + }); + + auto h_accum = Kokkos::create_mirror(accum); + Kokkos::deep_copy(h_accum, accum); + const double proj_err = + Kokkos::sqrt(h_accum(0)) / Kokkos::max(Kokkos::sqrt(h_accum(1)), EPS_DEN); + const double cons_err = + Kokkos::fabs(h_accum(2)) / Kokkos::max(Kokkos::fabs(h_accum(3)), EPS_DEN); + + return Errors{.proj_err = proj_err, .cons_err = cons_err}; +} + +#endif diff --git a/src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp new file mode 100644 index 00000000..449504a1 --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp @@ -0,0 +1,81 @@ +#ifndef PCMS_INTERPOLATOR_MASS_MATRIX_INTEGRATOR_H +#define PCMS_INTERPOLATOR_MASS_MATRIX_INTEGRATOR_H + +#include +#include +#include +#include +#include +#include +#include +#include + +// computes the mass matrix for each element +template +class MassMatrixIntegrator : public MeshField::Integrator +{ +public: + MassMatrixIntegrator(Omega_h::Mesh mesh_in, FieldElement& fe_in, + int order = 2) + : mesh(mesh_in), + fe(fe_in), + subMatrixSize(3 * 3), // FIXME remove hard coded size + elmMassMatrix("elmMassMatrix", mesh_in.nelems() * 3 * 3), + Integrator(order) + { + Kokkos::deep_copy(elmMassMatrix, 0); + assert(mesh.dim() == 2); // TODO support 1d,2d,3d + assert(mesh.family() == OMEGA_H_SIMPLEX); + } + void atPoints(Kokkos::View p, + Kokkos::View w, + Kokkos::View dV) + { + // std::cerr << "MassMatrixIntegrator::atPoints(...)\n"; + const size_t numPtsPerElem = p.extent(0) / mesh.nelems(); + // std::cerr << " Number points per Elem : " << numPtsPerElem << "\n"; + assert(numPtsPerElem >= 1); + const size_t ptDim = p.extent(1); + assert(ptDim == fe.MeshEntDim + 1); + Kokkos::parallel_for( + "eval", mesh.nelems(), KOKKOS_CLASS_LAMBDA(const int& elm) { + const auto first = elm * numPtsPerElem; + const auto last = first + numPtsPerElem; + for (auto pt = first; pt < last; pt++) { + // FIXME better way to fill? pass kokkos::subview to getValues? + Kokkos::Array + localCoord; + for (auto i = 0; i < localCoord.size(); i++) { + localCoord[i] = p(pt, i); + } + const auto N = fe.shapeFn.getValues(localCoord); + const auto wPt = w(pt); + const auto dVPt = dV(pt); + // printf("Shape Functions: %f, %f, %f \n", N[0], N[1], N[2]); + // printf("wPt, dVPt: %f, %f \n", wPt, dVPt); + for (auto i = 0; i < N.size(); i++) { + for (auto j = 0; j < N.size(); j++) { + elmMassMatrix(elm * subMatrixSize + i * 3 + j) += + N[i] * N[j] * wPt * dVPt; + } + } + } + }); + } + Omega_h::Mesh& mesh; + FieldElement& fe; + const int subMatrixSize; + Kokkos::View + elmMassMatrix; // numNodes^2 entries per element +}; + +template +Kokkos::View buildMassMatrix(Omega_h::Mesh& mesh, + FieldElement& coordFe) +{ + MassMatrixIntegrator mmi(mesh, coordFe); + mmi.process(coordFe); + return mmi.elmMassMatrix; +} + +#endif From e3e7fc1c5e4e4974eef6236cb9a91e89e4a601d5 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 15:35:26 -0400 Subject: [PATCH 05/29] add doxygen comments --- .../loadVectorIntegrator.hpp | 213 ++++++++++++++++-- 1 file changed, 195 insertions(+), 18 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp index 69e90615..660d66a9 100644 --- a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp @@ -8,12 +8,35 @@ #include "mesh_intersection.hpp" #include -// computes the load vector for each element -// This routine is only applicable linear elements -// It calculates the rhs of the conservative field transfer formula from source -// to the target mesh +/** + * @brief Computes the load vector for each target element in the conservative + * field transfer. + * + * This routine is used for constructing the right-hand side (RHS) of the + * conservative field transfer formulation, projecting field quantities from the + * source mesh to the target mesh. + * + * The underlying algorithm computes contributions to the load vector + * using geometric intersection data between source and target elements. + * + * @note Currently this method works for a two-dimensional linear triangles. + */ + +/** + * @brief Converts barycentric coordinates to global (physical) coordinates. + * + * Given barycentric coordinates within a 2D triangle and the coordinates of + * the triangle's vertices, this function computes the corresponding global + * position. + * + * @param barycentric_coord The barycentric coordinates \f$(\lambda_1, + * \lambda_2, \lambda_3)\f$ of the point. + * @param verts_coord The coordinates of the triangle's three vertices in global + * space. + * @return The 2D global coordinates corresponding to the given barycentric + * position. + */ -static const double EPS_AREA = 1e-16; [[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<2> global_from_barycentric( const MeshField::Vector3& barycentric_coord, const Omega_h::Few, 3>& verts_coord) @@ -27,6 +50,24 @@ static const double EPS_AREA = 1e-16; return real_coords; } +/** + * @brief Computes the barycentric coordinates of a 2D point with respect to a + * triangle. + * + * Given a point in global (x, y) coordinates and the coordinates of the three + * vertices of a triangle, this function evaluates the barycentric coordinates + * \f$(\lambda_1, \lambda_2, \lambda_3)\f$ of the point with respect to that + * triangle. + * + * @param point The 2D global coordinates of the point to evaluate (in + * Omega_h::Vector<2> format). + * @param verts_coord The vertex coordinates of the triangle (in r3d::Vector<2> + * format). + * @return A vector of three barycentric coordinates corresponding to the input + * point. + * + */ + [[nodiscard]] OMEGA_H_INLINE Omega_h::Vector<3> evaluate_barycentric( const Omega_h::Vector<2>& point, const r3d::Few, 3>& verts_coord) @@ -43,6 +84,26 @@ static const double EPS_AREA = 1e-16; return barycentric_coordinate; } +/** + * @brief Evaluates the value of a linear function at a given point using + * barycentric coordinates. + * + * This function computes the interpolated value of a nodal scalar field over a + * triangle, using barycentric coordinates within the specified element. + * + * @param nodal_values The global array of nodal field values. + * @param faces2nodes The element-to-node connectivity array. + * @param bary_coords The barycentric coordinates of the evaluation point within + * the triangle. + * @param elm_id The ID of the triangle element being evaluated. + * @return The interpolated function value at the given point. + * + * @note This function assumes linear (3-node) triangular element. + */ + +[[nodiscard]] OMEGA_H_INLINE double evaluate_function_value( + const Omega_h::Reals& nodal_values, const Omega_h::LOs& faces2nodes, + const Omega_h::Vector<3>& bary_coords, const int elm_id); [[nodiscard]] OMEGA_H_INLINE double evaluate_function_value( const Omega_h::Reals& nodal_values, const Omega_h::LOs& faces2nodes, const Omega_h::Vector<3>& bary_coords, const int elm_id) @@ -58,12 +119,33 @@ static const double EPS_AREA = 1e-16; return value; } +/** + * @brief Provides barycentric integration points and weights for a triangle + * element. + * + * This templated struct stores the barycentric coordinates + * and quadrature weights for performing numerical integration over a reference + * triangle. It is used for integrating functions over elements in the + * conservative field transfer. + * + * @tparam order The quadrature order (number of integration points and + * polynomial accuracy). + */ template struct IntegrationData { + // Barycentric coordinates of integration points Kokkos::View bary_coords; + + // Quadrature weights associated with each integration point Kokkos::View weights; + /** + * @brief Constructs the integration data for a given quadrature order + * + * Initializes barycentric coordinates and weights using + * MeshField's predefined triangle quadrature rules. + */ IntegrationData() { auto ip_vec = MeshField::getIntegrationPoints(MeshField::Triangle, order); @@ -84,9 +166,38 @@ struct IntegrationData Kokkos::deep_copy(weights, weights_host); } + /** + * @brief Returns the number of integration points + * + * @return Number of integration points for the selected order. + */ + int size() const { return bary_coords.extent(0); } }; +/** + * @brief Computes the counter-clockwise (CCW) vertex ordering of a 2D polygon. + * + * Given a 2D polygon stored in an `r3d::Polytope<2>` struct where each vertex + * has exactly two neighbors. This function reconstructs the CCW traversal order + * of the polygon's vertices and stores the result in the provided `order` + * array. + * + * The traversal begins at vertex 0 and proceeds by selecting the neighbor that + * is not previously visited vertex, thereby completing a full loop around the + * polygon. + * + * @param poly The polygon to process, represented as an `r3d::Polytope<2>`. + * Each vertex includes a list of two neighbor indices + * (`pnbrs[2]`) forming the cycle. + * + * @param[out] order An output array to store the CCW vertex order. Must be + * preallocated to hold at least `poly.nverts` size. + * @return The number of vertices in the polygon (i.e., `poly.nverts`) + * + * @see r3d::Polytope + */ + [[nodiscard]] OMEGA_H_INLINE int get_polygon_cycle_ccw( const r3d::Polytope<2>& poly, int* order) { @@ -107,6 +218,45 @@ struct IntegrationData return m; } +/** + * @brief Computes the per-element RHS load vectors for conservative field + * projection from source to target mesh. + * + * This function computes local (element-wise) right-hand side (RHS) + * contributions for the Galerkin projection of a scalar field from the source + * mesh to the target mesh. It integrates over the polygonal intersection + * regions between each target element and its intersecting source elements + * using barycentric quadrature. + * + * The output is a flat array containing unassembled load vector contributions + * at the nodes of each target triangle. + * + * @param target_mesh The target mesh object receiving the projected scalar + * field. + * @param source_mesh The source mesh object 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 Scalar field values defined at the nodes of the source + * mesh. + * + * @return A Kokkos view containing per-element load vectors. + * Each triangle contributes 3 values (one per node), so the view has + * size 3 × (number of target elements). + * + * @note + * - This function assumes 2D linear triangular elements. + * - Degenerate or near-zero-area intersection polygons are skipped. + * - Each polygon is triangulated using a fan structure and integrated using + * barycentric quadrature rules. + * - The returned vector must be assembled into a global RHS vector in a later + * step. + * + * @see evaluate_barycentric, evaluate_function_value, global_from_barycentric + * @see IntersectionResults + */ + Kokkos::View buildLoadVector( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values) @@ -133,7 +283,6 @@ Kokkos::View buildLoadVector( get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, elm); const int start = intersection.tgt2src_offsets[elm]; const int end = intersection.tgt2src_offsets[elm + 1]; - // printf("number of intersections : %d\n", end-start); Omega_h::Vector<3> part_integration = {0.0, 0.0, 0.0}; for (int i = start; i < end; ++i) { @@ -147,10 +296,6 @@ Kokkos::View buildLoadVector( int order[r3d::MaxVerts<2>::value]; auto m = get_polygon_cycle_ccw(poly, order); auto poly_area = r3d::measure(poly); - // printf("No of verts : %d\n", nverts); - // decompose if more than 3 vertices - // printf(" area from polytope : %f\n", - //r3d::measure(poly)); double sum_area = 0; for (int j = 1; j < nverts - 1; ++j) { // build triangle from poly.verts[order[0]], poly.verts[order[j]], @@ -174,9 +319,9 @@ Kokkos::View buildLoadVector( Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); sum_area += area; + const double EPS_AREA = abs_tol + rel_tol * poly_area; if (area <= EPS_AREA) continue; // drops duplicates and colinear/degenerates - // printf(" area = %f\n", area); for (int ip = 0; ip < npts; ++ip) { auto bary = integrationPoints.bary_coords(ip); @@ -201,9 +346,6 @@ Kokkos::View buildLoadVector( } } } - // printf(" sum_area = %.16e\n", sum_area); - // printf(" poly_area = %.16e\n", poly_area); - // OMEGA_H_CHECK(nearly_equal(sum_area, poly_area)); } for (int j = 0; j < 3; ++j) { @@ -214,10 +356,46 @@ Kokkos::View buildLoadVector( return elmLoadVector; } +/** + * @brief Stores scalar error metrics for supermesh-based field projection + * diagnostics. + * + * This structure holds two exactly computable error measures that quantify + * the accuracy and conservation properties of field transfer between + * nonconforming meshes. Both metrics are evaluated on the *supermesh*—the + * polygonal intersection mesh constructed between the source and target meshes. + * + * - **`proj_err`** — The L2 projection error between the projected source and + * target fields, integrated over the supermesh. By definition of the supermesh, + * both fields are representable in the same function superspace, allowing exact + * evaluation (up to roundoff) of + * \f$ g = P_S^D(q_D) - P_S^T(q_T) = q_D - q_T \f$. + * The Galerkin projection being optimal in the L2 norm, this value measures + * the exact difference between the donor and recipient fields. + * + * - **`cons_err`** — The relative conservation error, representing the + * imbalance between the total scalar quantities of the projected source and + * projected target fields over the supermesh domain. + * + * The use of the supermesh provides a common integration space that contains + * the basis functions of both meshes. As a result, both projection and + * conservation errors are exactly computable without additional search or + * remapping operations. Only the evaluation of the parent basis functions is + * required. + * + * @note + * - Errors are computed on the supermesh, not on the original target mesh. + * - The projection error corresponds to the L2 norm difference between source + * and target representations. + * - The conservation error measures the total integral imbalance over the + * supermesh. + * + * @see evaluate_proj_and_cons_errors, IntersectionResults + */ struct Errors { - double proj_err; - double cons_err; + double proj_err; ///< L2 projection error computed on the supermesh. + double cons_err; ///< Relative conservation error over the supermesh. }; Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, @@ -287,10 +465,9 @@ Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); sum_area += area; - const double EPS_AREA = Kokkos::fmax(abs_tol, rel_tol * area); + const double EPS_AREA = abs_tol + rel_tol * ploy_area; if (area <= EPS_AREA) continue; // drops duplicates and colinear/degenerates - // printf(" area = %f\n", area); for (int ip = 0; ip < npts; ++ip) { auto bary = integrationPoints.bary_coords(ip); From 53b6498e0d150e8b9f38addc9f056476c4b08ca7 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 15:51:30 -0400 Subject: [PATCH 06/29] add doxygen comments --- .../loadVectorIntegrator.hpp | 101 +++++++++++------- 1 file changed, 65 insertions(+), 36 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp index 660d66a9..730114ce 100644 --- a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp @@ -356,48 +356,77 @@ Kokkos::View buildLoadVector( return elmLoadVector; } -/** - * @brief Stores scalar error metrics for supermesh-based field projection - * diagnostics. - * - * This structure holds two exactly computable error measures that quantify - * the accuracy and conservation properties of field transfer between - * nonconforming meshes. Both metrics are evaluated on the *supermesh*—the - * polygonal intersection mesh constructed between the source and target meshes. - * - * - **`proj_err`** — The L2 projection error between the projected source and - * target fields, integrated over the supermesh. By definition of the supermesh, - * both fields are representable in the same function superspace, allowing exact - * evaluation (up to roundoff) of - * \f$ g = P_S^D(q_D) - P_S^T(q_T) = q_D - q_T \f$. - * The Galerkin projection being optimal in the L2 norm, this value measures - * the exact difference between the donor and recipient fields. - * - * - **`cons_err`** — The relative conservation error, representing the - * imbalance between the total scalar quantities of the projected source and - * projected target fields over the supermesh domain. - * - * The use of the supermesh provides a common integration space that contains - * the basis functions of both meshes. As a result, both projection and - * conservation errors are exactly computable without additional search or - * remapping operations. Only the evaluation of the parent basis functions is - * required. - * - * @note - * - Errors are computed on the supermesh, not on the original target mesh. - * - The projection error corresponds to the L2 norm difference between source - * and target representations. - * - The conservation error measures the total integral imbalance over the - * supermesh. - * - * @see evaluate_proj_and_cons_errors, IntersectionResults - */ +/// Holds projection and conservation error metrics returned by +/// evaluate_pro_and_cons_errors(). struct Errors { double proj_err; ///< L2 projection error computed on the supermesh. double cons_err; ///< Relative conservation error over the supermesh. }; +/** + * @brief Computes projection and conservation errors over the supermesh for + * scalar field transfer. + * + * This function quantifies the accuracy and conservation properties of + * conservative field transfer between nonconforming meshes using + * supermesh-based integration. It returns a struct containing two error + * metrics: + * + * - **Projection Error (`proj_err`)** — Measures the L2 norm of the difference + * between the projected source field and the target field over the + * supermesh.This reflects how accurately the field has been projected. + * + * - **Conservation Error (`cons_err`)** — Relative difference in the integrated + * field values between source and target representations. This captures + * conservation loss across the transfer. + * + * ### Mathematical Definitions: + * \f[ + * \text{proj\_err} = \frac{ \| q_D - q_T \|_{L_2(\Omega_S)} } + * { \| q_D \|_{L_2(\Omega_S)} }, \quad + * \text{cons\_err} = \frac{ \left| \int_{\Omega_S} q_D - \int_{\Omega_S} q_T + * \right| } { \left| \int_{\Omega_S} q_D \right| } + * \f] + * + * where: + * - \f$q_D\f$ is the scalar field defined on the source mesh (mesh from where + * the field is defined), + * - \f$q_T\f$ is the projected field on the target mesh (mesh to where the + * field is projected), + * - \f$\Omega_S\f$ is the supermesh formed by polygonal intersections of source + * and target elements. + * + * Integration is performed by triangulating each intersection region and + * applying barycentric quadrature. Degenerate or near-zero-area triangles are + * skipped based on area tolerance. + * + * + * @param target_mesh The target mesh object receiving the projected scalar + * field. + * @param source_mesh The source mesh object 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 target_values Nodal scalar field values evaluated on the target mesh + * using galerkin projection. + * @param source_values Scalar field values defined at the nodes of the source + * mesh. + * + * + * @return A struct containing: + * - `proj_err`: Exact L2 projection error over the supermesh. + * - `cons_err`: Relative conservation error over the supermesh. + * + * @note + * - Assumes 2D linear (P1) triangular elements. + * - Ideal for validating conservative transfer schemes or testing projection + * fidelity. + * + * @see IntersectionResults, buildLoadVector + */ + Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, From 3dea928abaac5e91e042c15fb16b1176f0d62506 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 16:07:24 -0400 Subject: [PATCH 07/29] add doxygen comments --- .../mesh_intersection/calculateLoadVector.hpp | 64 +++++++++++-------- 1 file changed, 38 insertions(+), 26 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp b/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp index 6aa2219d..bf7cdab8 100644 --- a/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp @@ -1,13 +1,18 @@ /** * @file calculateLoadVector.hpp - * @brief Functions for calculating load vector on finite element meshes + * @brief Utilities for assembling load vectors in finite element field transfer + * routines. * - * This file contains functions for creating and computing load vectors - * for finite element calculations using Omega_h mesh structures and PETSc. + * This header defines functions for constructing global load vectors for finite + * element meshes required for conservative field transfer methods between + * source and target meshes. It leverages Omega_h mesh data structures and PETSc + * for vector assembly and storage. + * + * Typical usage involves computing element-wise load vector contributions + * and assembling them into a global PETSc vector in COO format. */ - -#ifndef COMPUTING_AT_SCALE_CALCULATE_LOAD_VECTOR_HPP -#define COMPUTING_AT_SCALE_CALCULATE_LOAD_VECTOR_HPP +#ifndef PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP +#define PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP #include #include #include //Omega_h::atomic_fetch_add @@ -37,17 +42,36 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; /** - * @brief Calculates the load + * @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 + * @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. * - * This function constructs a mass matrix based on the provided mesh using - * a finite element approach. It creates coordinate field elements, builds - * the mass matrix using the massMatrixIntegrator, and sets up the PETSc matrix - * with appropriate values. + * @return PetscErrorCode Returns PETSC_SUCCESS if successful, or an appropriate + * PETSc error code otherwise. * - * @param mesh The Omega_h mesh to calculate the mass matrix for - * @param[out] mass_out Pointer to the resulting mass matrix - * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) + * @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 */ + inline PetscErrorCode calculateLoadVector( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values, @@ -72,18 +96,6 @@ inline PetscErrorCode calculateLoadVector( auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector); Kokkos::deep_copy(hostElmLoadVector, elmLoadVector); - /* - std::cerr << " DEBUG: printing the load vector \n"; - for (int i = 0; i < hostElmLoadVector.size(); ++i) { - std::cout << hostElmLoadVector[i] << " "; - } - std::cout << "\n"; - - target_mesh.add_tag(2, "elmLoadVector", 3, - Omega_h::read(Omega_h::Write(elmLoadVector))); - - Omega_h::vtk::write_parallel("loadVector.vtk", &target_mesh, 2); - */ PetscInt idx = 0; for (PetscInt e = 0; e < target_mesh.nelems(); ++e) { for (PetscInt vi = 0; vi < numNodesPerTri; ++vi) { From 6b5c25ee2b07b5bb41eb65a4151009bfe3178994 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 16:29:44 -0400 Subject: [PATCH 08/29] add doxygen comments --- .../conservativeProjectionSolver.hpp | 204 ++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp diff --git a/src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp b/src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp new file mode 100644 index 00000000..1068929e --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp @@ -0,0 +1,204 @@ +#ifndef PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP +#ifndef PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP + +#include +#include +#include +#include +#include +#include +#include +#include // Must be included before other PETSc headers +#include +#include +#include +#include + +#include "calculateLoadVector.hpp" +#include "calculateMassMatrix.hpp" + +/** + * @brief Solves a linear system Ax = b using PETSc's KSP solvers + * + * Uses PETSc's Krylov Subspace solvers to find x in Ax = b. + * The solver can be configured through PETSc runtime options. + * + * @param A The system matrix + * @param b The right-hand side vector + * @return Vec Solution vector x + */ +static Vec solveLinearSystem(Mat A, Vec b) +{ + PetscInt m, n; + PetscErrorCode ierr; + + ierr = MatGetSize(A, &m, &n); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Vec x; + ierr = VecCreateSeqKokkos(PETSC_COMM_SELF, n, &x); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + KSP ksp; + ierr = KSPCreate(PETSC_COMM_SELF, &ksp); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = KSPSetOperators(ksp, A, A); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = KSPSetComputeSingularValues(ksp, PETSC_TRUE); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = KSPSetFromOptions(ksp); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = KSPSetUp(ksp); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = KSPSolve(ksp, b, x); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + /// compute and print condition number estimate + PetscReal smax = 0.0, smin = 0.0; + ierr = KSPComputeExtremeSingularValues(ksp, &smax, &smin); + if (!ierr && smin > 0.0) { + PetscPrintf(PETSC_COMM_WORLD, + "Estimated condition number of matrix A: %.6e\n", smax / smin); + } else { + PetscPrintf(PETSC_COMM_WORLD, + "Condition number estimate unavailable (smin <= 0 or error)\n"); + } + + ierr = KSPDestroy(&ksp); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + return x; +} + +/** + * @brief Solves a conservative projection problem to transfer scalar field + * values onto a target mesh. + * + * This function assembles and solves a linear system of the form: + * \f[ + * M \cdot x = f + * \f] + * where: + * - \f$M\f$ is the mass matrix on the target mesh (based on P1 finite + * elements), + * - \f$f\f$ is the load vector computed on the supermesh, + * - \f$x\f$ is the unknown nodal field on the target mesh (solution). + * + * The method computes the conservative field transfer between two nonconforming + * meshes using mesh intersections (supermesh). + * + * ### Algorithm Steps: + * 1. Compute and assemble mass matrix and load vector + * 2. Solve the linear system using PETSc. + * 3. Return the solution as a nodal field on the target mesh. + * + * @param target_mesh The Omega_h mesh where the field is projected. + * @param source_mesh The Omega_h mesh containing the original field data. + * @param intersection Precomputed intersection information between source and + * target meshes. + * @param source_values Nodal scalar field values on the source mesh. + * + * @return A vector of nodal values on the target mesh after projection + * (Omega_h::Reals). + * + * + */ +Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& source_values) +{ + + if ((PetscInt)source_values.size() != + source_mesh.coords().size() / source_mesh.dim()) { + std::cerr << "ERROR: source_values size (" << source_values.size() + << ") doesn't match expected size (" + << source_mesh.coords().size() / source_mesh.dim() << ")" + << std::endl; + throw std::runtime_error("source_values length mismatch"); + } + + Mat mass; + PetscErrorCode ierr = calculateMassMatrix(target_mesh, &mass); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Vec vec; + ierr = calculateLoadVector(target_mesh, source_mesh, intersection, + source_values, &vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Vec x = solveLinearSystem(mass, vec); + + Omega_h::Write solution_vector( + target_mesh.nverts(), 0, "stores the solution coefficients"); + + PetscScalar* array; + + ierr = VecGetArray(x, &array); + + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + auto solution_host = Omega_h::HostWrite(target_mesh.nverts()); + + for (PetscInt i = 0; i < target_mesh.nverts(); ++i) { + solution_host[i] = array[i]; + } + + solution_vector = Omega_h::Write(solution_host); + ierr = VecRestoreArray(x, &array); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = VecDestroy(&x); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = MatDestroy(&mass); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = VecDestroy(&vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + return Omega_h::read(solution_vector); +} + +Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& source_values) +{ + Vec vec; + PetscErrorCode ierr; + ierr = calculateLoadVector(target_mesh, source_mesh, intersection, + source_values, &vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Omega_h::Write rhsvector(target_mesh.nverts(), 0, + "stores the rhs vector"); + { + PetscErrorCode ierr; + PetscScalar* array; + ierr = VecGetArray(vec, &array); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + auto rhsvec_host = Omega_h::HostWrite(target_mesh.nverts()); + + for (PetscInt i = 0; i < target_mesh.nverts(); ++i) { + rhsvec_host[i] = array[i]; + } + + rhsvector = Omega_h::Write(rhsvec_host); + ierr = VecRestoreArray(vec, &array); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + } + + ierr = VecDestroy(&vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + return Omega_h::read(rhsvector); +} + +#endif // PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP From 2ef3545ab69d59f2914bce745412e3275c0aa713 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 16:30:06 -0400 Subject: [PATCH 09/29] minor fixes --- .../mesh_intersection/mesh_intersection.hpp | 219 ++++++++++++++++++ 1 file changed, 219 insertions(+) create mode 100644 src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp diff --git a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp new file mode 100644 index 00000000..1fc7c0fa --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp @@ -0,0 +1,219 @@ +#ifndef PCMS_INTERPOLATOR_MESH_INTERSECTION_HPP +#define PCMS_INTERPOLATOR_MESH_INTERSECTION_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +constexpr static double abs_tol = 1e-18; +constexpr static double rel_tol = 1e-12; + +struct IntersectionResults +{ + Omega_h::LOs tgt2src_offsets; + Omega_h::LOs tgt2src_indices; +}; + +class FindIntersections +{ +private: + Omega_h::Mesh& source_mesh_; + Omega_h::Mesh& target_mesh_; + +public: + FindIntersections(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh) + : source_mesh_(source_mesh), target_mesh_(target_mesh) {}; + + void adjBasedIntersectSearch(const Omega_h::LOs& tgt2src_offsets, + Omega_h::Write& nIntersections, + Omega_h::Write& tgt2src_indices, + bool is_count_only); +}; + +[[nodiscard]] OMEGA_H_INLINE r3d::Few, 3> +get_vert_coords_of_elem(const Omega_h::Reals& coords, + const Omega_h::LOs& faces2nodes, const int id) +{ + + const auto elm_verts = Omega_h::gather_verts<3>(faces2nodes, id); + + const Omega_h::Matrix<2, 3> elm_vert_coords = + Omega_h::gather_vectors<3, 2>(coords, elm_verts); + + r3d::Few, 3> r3d_vector; + for (int i = 0; i < 3; ++i) { + + r3d_vector[i][0] = elm_vert_coords[i][0]; + r3d_vector[i][1] = elm_vert_coords[i][1]; + } + + return r3d_vector; +} + +inline Kokkos::View compute_centroid(Omega_h::Mesh& mesh) +{ + + const auto& mesh_coords = mesh.coords(); + const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& nfaces = mesh.nfaces(); + + Kokkos::View cell_centroids("centroid coordinates", + nfaces); + Omega_h::parallel_for( + "calculate the centroid in each tri element", nfaces, + OMEGA_H_LAMBDA(const Omega_h::LO id) { + const auto current_el_verts = Omega_h::gather_verts<3>(faces2nodes, id); + const Omega_h::Few, 3> current_el_vert_coords = + Omega_h::gather_vectors<3, 2>(mesh_coords, current_el_verts); + auto centroid = Omega_h::average(current_el_vert_coords); + cell_centroids(id, 0) = centroid[0]; + cell_centroids(id, 1) = centroid[1]; + }); + + return cell_centroids; +} + +void FindIntersections::adjBasedIntersectSearch( + const Omega_h::LOs& tgt2src_offsets, + Omega_h::Write& nIntersections, + Omega_h::Write& tgt2src_indices, bool is_count_only) +{ + + const auto& tgt_coords = target_mesh_.coords(); + const auto& src_coords = source_mesh_.coords(); + const auto& tgt_faces2nodes = + target_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_faces2nodes = + source_mesh_.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& src_elem_areas = measure_elements_real(&source_mesh_); + const auto& tgt_elem_areas = measure_elements_real(&target_mesh_); + const auto& t2t = + source_mesh_.ask_dual(); // gives connected element neighbors + const auto& t2tt = t2t.a2ab; + const auto& tt2t = t2t.ab2b; + + auto centroids = compute_centroid(target_mesh_); + + pcms::GridPointSearch search_cell(source_mesh_, 20, 20); + auto results = search_cell(centroids); + + auto nfaces_target = target_mesh_.nfaces(); + Omega_h::parallel_for( + nfaces_target, + OMEGA_H_LAMBDA(const Omega_h::LO id) { + Queue queue; + Track visited; + + auto current_cell_id = results(id).tri_id; + auto current_tgt_elm_area = tgt_elem_areas[id]; + + OMEGA_H_CHECK_PRINTF(current_cell_id >= 0, + "ERROR: source cell id not found for given target " + "centroid %d (%f, %f)\n", + id, centroids(id, 0), centroids(id, 1)); + + auto tgt_elm_vert_coords = + get_vert_coords_of_elem(tgt_coords, tgt_faces2nodes, id); + + Omega_h::LO start_counter; + if (!is_count_only) { + start_counter = tgt2src_offsets[id]; + } + + int count = 0; + + count++; + visited.push_back(current_cell_id); + queue.push_back(current_cell_id); + + if (!is_count_only) { + int idx_count = count - 1; + tgt2src_indices[start_counter + idx_count] = current_cell_id; + } + + while (!queue.isEmpty()) { + Omega_h::LO currentElm = queue.front(); + queue.pop_front(); + auto start = t2tt[currentElm]; + auto end = t2tt[currentElm + 1]; + + for (int i = start; i < end; ++i) { + auto neighborElmId = tt2t[i]; + + if (visited.notVisited(neighborElmId)) { + visited.push_back(neighborElmId); + auto elm_vert_coords = get_vert_coords_of_elem( + src_coords, src_faces2nodes, neighborElmId); + r3d::Polytope<2> intersection; + r3d::intersect_simplices(intersection, tgt_elm_vert_coords, + elm_vert_coords); + auto intersected_area = r3d::measure(intersection); + auto current_src_elm_area = src_elem_areas[neighborElmId]; + auto scale = + Kokkos::fmax(current_tgt_elm_area, current_src_elm_area); + auto eps = Kokkos::fmax(abs_tol, rel_tol * scale); + if (intersection.nverts >= 3 && intersected_area >= eps) { + count++; + + OMEGA_H_CHECK_PRINTF( + count < 500, "WARNING: count exceeds 500 for target %d", id); + + queue.push_back(neighborElmId); + + if (!is_count_only) { + Omega_h::LO idx_count = count - 1; + tgt2src_indices[start_counter + idx_count] = neighborElmId; + + } // end of tgt2src_indices check + + } // end of intersection with bbox check + + } // end of not visited check + + } // end of loop over adj elements to the current element + + } // end of while loop + + nIntersections[id] = count; + }, // end of lambda + "count the number of intersections for each target element"); +} + +IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, + Omega_h::Mesh& target_mesh) +{ + FindIntersections intersect(source_mesh, target_mesh); + + auto nfaces_target = target_mesh.nfaces(); + + Omega_h::Write nIntersections( + nfaces_target, 0, "number of intersections in each target vertex"); + + Omega_h::Write tgt2src_indices; + + intersect.adjBasedIntersectSearch(Omega_h::LOs(), nIntersections, + tgt2src_indices, true); + + Kokkos::fence(); + auto tgt2src_offsets = Omega_h::offset_scan(Omega_h::Read(nIntersections), + "offsets for intersections"); + auto ntotal_intersections = tgt2src_offsets.last(); + + Kokkos::fence(); + + tgt2src_indices = Omega_h::Write( + ntotal_intersections, 0, + "indices of the source elements that intersect the given target element"); + + intersect.adjBasedIntersectSearch(tgt2src_offsets, nIntersections, + tgt2src_indices, false); + return {.tgt2src_offsets = tgt2src_offsets, + .tgt2src_indices = Omega_h::read(tgt2src_indices)}; +} + +#endif From bdf44ebc1360f92e2fa19022b3bb9d6b3d2b9b56 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 18:20:11 -0400 Subject: [PATCH 10/29] update doxygen comments --- .../mesh_intersection/mesh_intersection.hpp | 62 ++++++++++++++++++- 1 file changed, 59 insertions(+), 3 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp index 1fc7c0fa..c4d35c34 100644 --- a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp +++ b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp @@ -10,9 +10,20 @@ #include #include -constexpr static double abs_tol = 1e-18; -constexpr static double rel_tol = 1e-12; - +constexpr static double abs_tol = 1e-18; /// abs tolerance +constexpr static double rel_tol = 1e-12; /// rel tolerance + +/** + * @brief Stores results of mesh element intersections for conservative + * transfer. + * + * Contains mappings from each target element to the list of source elements + * that intersect with it. Used to guide integration over overlapping regions. + * + * - `tgt2src_offsets[i]` is the offset into `tgt2src_indices` where source + * elements for target element `i` begin. + * - `tgt2src_indices` contains flattened indices of source elements per target. + */ struct IntersectionResults { Omega_h::LOs tgt2src_offsets; @@ -77,6 +88,31 @@ inline Kokkos::View compute_centroid(Omega_h::Mesh& mesh) return cell_centroids; } +/** + * @brief Performs adjacency-based intersection search between target and source + * elements. + * + * For each target element, starting from the source element that contains its + * centroid, a queue-based BFS traversal is used over the adjacency graph of + * source elements. If an element intersects the target triangle (based on area + * tolerance), it is included. + * + * @param tgt2src_offsets Offsets array (only used when writing indices). + * @param[out] nIntersections Number of intersecting source elements per target + * element. + * @param[out] tgt2src_indices Indices of intersecting source elements. + * @param is_count_only If true, only counts intersections; if false, also fills + * tgt2src_indices. + * + * @note This method assumes 2D linear triangles and uses + * `r3d::intersect_simplices` for geometric intersection. + * + * @see r3d::intersect_simplices, intersectTargets + */ +void FindIntersections::adjBasedIntersectSearch( + const Omega_h::LOs& tgt2src_offsets, + Omega_h::Write& nIntersections, + Omega_h::Write& tgt2src_indices, bool is_count_only); void FindIntersections::adjBasedIntersectSearch( const Omega_h::LOs& tgt2src_offsets, @@ -184,6 +220,26 @@ void FindIntersections::adjBasedIntersectSearch( "count the number of intersections for each target element"); } +/** + * @brief Computes source-target element intersections for conservative + * projection. + * + * For each target element in the target mesh, this function identifies source + * elements from the source mesh that geometrically intersect with it using an + * adjacency-based breadth-first search strategy. The result is returned as a + * compact mapping. + * + * @param source_mesh The source Omega_h mesh. + * @param target_mesh The target Omega_h mesh. + * @return An IntersectionResults struct containing target-to-source mapping + * data. + * + * @note The intersection test is done using 2D polygon intersection routines + * from r3d. Only valid (non-degenerate) polygonal intersections are included. + * + * @see FindIntersections::adjBasedIntersectSearch + */ + IntersectionResults intersectTargets(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh) { From 93a08d5787e67b483fe61660e790200333b8c11e Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 18:52:52 -0400 Subject: [PATCH 11/29] update comments --- ...adVector.hpp => calculate_load_vector.hpp} | 19 +++++++-------- ...ssMatrix.hpp => calculate_mass_matrix.hpp} | 14 ----------- ...hpp => conservative_projection_solver.hpp} | 23 ++++++++++++++++--- ...egrator.hpp => load_vector_integrator.hpp} | 20 ++++++++++++++++ ...egrator.hpp => mass_matrix_integrator.hpp} | 0 5 files changed, 50 insertions(+), 26 deletions(-) rename src/pcms/interpolator/mesh_intersection/{calculateLoadVector.hpp => calculate_load_vector.hpp} (87%) rename src/pcms/interpolator/mesh_intersection/{calculateMassMatrix.hpp => calculate_mass_matrix.hpp} (93%) rename src/pcms/interpolator/mesh_intersection/{conservativeProjectionSolver.hpp => conservative_projection_solver.hpp} (89%) rename src/pcms/interpolator/mesh_intersection/{loadVectorIntegrator.hpp => load_vector_integrator.hpp} (96%) rename src/pcms/interpolator/mesh_intersection/{massMatrixIntegrator.hpp => mass_matrix_integrator.hpp} (100%) diff --git a/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp similarity index 87% rename from src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp rename to src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp index bf7cdab8..8077e85f 100644 --- a/src/pcms/interpolator/mesh_intersection/calculateLoadVector.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp @@ -1,16 +1,16 @@ /** - * @file calculateLoadVector.hpp - * @brief Utilities for assembling load vectors in finite element field transfer - * routines. + * @file calculate_load_vector.hpp + * @brief Routines for assembling global load vector in conservative field + * projection. * - * This header defines functions for constructing global load vectors for finite - * element meshes required for conservative field transfer methods between - * source and target meshes. It leverages Omega_h mesh data structures and PETSc - * for vector assembly and storage. + * Provides functionality to compute and assemble the global load vector + * used in Galerkin-based conservative field transfer between non-conforming + * meshes. * - * Typical usage involves computing element-wise load vector contributions - * and assembling them into a global PETSc vector in COO format. + * @created by Abhiyan Paudel + * @date August 2025 */ + #ifndef PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP #define PCMS_INTERPOLATOR_CALCULATE_LOAD_VECTOR_HPP #include @@ -55,6 +55,7 @@ using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; * 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. diff --git a/src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp similarity index 93% rename from src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp rename to src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp index 71588503..96a1202e 100644 --- a/src/pcms/interpolator/mesh_intersection/calculateMassMatrix.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp @@ -129,20 +129,6 @@ inline PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out) auto host_elmMassMatrix = Kokkos::create_mirror_view(elmMassMatrix); - /* - std::cerr << " DEBUG: printing the element mass matrix \n"; - for (int i = 0; i < host_elmMassMatrix.size(); ++i) { - std::cout << host_elmMassMatrix[i] << " "; - } - - */ - - // std::cout << "\n"; - // mesh.add_tag(2, "elmMassMatrix", 3 * 3, - // Omega_h::read(Omega_h::Write(elmMassMatrix))); - - // Omega_h::vtk::write_parallel("massMatrix.vtk", &mesh, 2); - Mat mass; PetscCall(CreateMatrix(mesh, &mass)); PetscBool is_kokkos; diff --git a/src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp similarity index 89% rename from src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp rename to src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp index 1068929e..0c167d09 100644 --- a/src/pcms/interpolator/mesh_intersection/conservativeProjectionSolver.hpp +++ b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp @@ -1,3 +1,20 @@ +/** + * @file conservative_projection_solver.hpp + * @brief Solves the conservative projection of scalar fields between + * non-matching meshes. + * + * Provides the main interface to perform Galerkin projection of scalar fields + * from a source mesh to a target mesh using conservative transfer using a + * supermesh generated from mesh intersections. + * + * The solver computes the right-hand side (load vector), assembles the mass + * matrix, and solves the resulting linear system to obtain projected nodal + * values. + * + * @created by Abhiyan Paudel + * @date August 2025 + */ + #ifndef PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP #ifndef PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP @@ -76,8 +93,8 @@ static Vec solveLinearSystem(Mat A, Vec b) } /** - * @brief Solves a conservative projection problem to transfer scalar field - * values onto a target mesh. + * @brief Solves a conservative galerkin projection problem to transfer scalar + * field values onto a target mesh. * * This function assembles and solves a linear system of the form: * \f[ @@ -89,7 +106,7 @@ static Vec solveLinearSystem(Mat A, Vec b) * - \f$f\f$ is the load vector computed on the supermesh, * - \f$x\f$ is the unknown nodal field on the target mesh (solution). * - * The method computes the conservative field transfer between two nonconforming + * The method computes the conservative field transfer between two non-matching * meshes using mesh intersections (supermesh). * * ### Algorithm Steps: diff --git a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp similarity index 96% rename from src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp rename to src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp index 730114ce..2fd1a859 100644 --- a/src/pcms/interpolator/mesh_intersection/loadVectorIntegrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp @@ -1,3 +1,23 @@ +/** + * @file load_vector_integrator.hpp + * @brief Functions for computing load vectors in conservative field projection. + * + * This file implements routines to compute the element-wise load vector + * (right-hand side) contributions used in Galerkin projection of scalar + * fields from a source mesh to a target mesh. + * + * The integration is performed over polygonal intersections (supermesh) between + * source and target elements using barycentric quadrature. The resulting values + * represent unassembled local contributions that can later be combined into a + * global load vector. + * + * @note + * - Assumes 2D linear triangular meshes. + * - Intersection data is provided via the `IntersectionResults` structure. + * + * @created by Abhiyan Paudel + * @date August 2025 + */ #ifndef PCMS_INTERPOLATOR_LOAD_VECTOR_INTEGRATOR_HPP #define PCMS_INTERPOLATOR_LOAD_VECTOR_INTEGRATOR_HPP diff --git a/src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp b/src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp similarity index 100% rename from src/pcms/interpolator/mesh_intersection/massMatrixIntegrator.hpp rename to src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp From aabe4023952fe8d6e8be04e991bd1d341ac28927 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 21:58:57 -0400 Subject: [PATCH 12/29] update CMakeLists --- .../mesh_intersection/CMakeLists.txt | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 src/pcms/interpolator/mesh_intersection/CMakeLists.txt diff --git a/src/pcms/interpolator/mesh_intersection/CMakeLists.txt b/src/pcms/interpolator/mesh_intersection/CMakeLists.txt new file mode 100644 index 00000000..e0bee537 --- /dev/null +++ b/src/pcms/interpolator/mesh_intersection/CMakeLists.txt @@ -0,0 +1,45 @@ +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 + $ + $ +) + +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 ${PETSC_INCLUDES}) + target_link_libraries(pcms_mesh_intersection PUBLIC ${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 + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) + +install(EXPORT meshIntersectionTargets + FILE pcms_mesh_intersection-targets.cmake + NAMESPACE pcms:: + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms) From a3330efa4fbc71d7779c48288757c81a0d771282 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 21:59:40 -0400 Subject: [PATCH 13/29] update CMakeLists --- src/pcms/interpolator/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pcms/interpolator/CMakeLists.txt b/src/pcms/interpolator/CMakeLists.txt index 5319a4ea..c8cf349a 100644 --- a/src/pcms/interpolator/CMakeLists.txt +++ b/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) @@ -21,7 +23,7 @@ set_target_properties(pcms_interpolator PROPERTIES OUTPUT_NAME pcmsinterpolator EXPORT_NAME interpolator) add_library(pcms::interpolator ALIAS pcms_interpolator) -target_link_libraries(pcms_interpolator PUBLIC pcms::core Kokkos::kokkoskernels) +target_link_libraries(pcms_interpolator PUBLIC pcms::core Kokkos::kokkoskernels pcms_mesh_intersection) install( TARGETS pcms_interpolator From 7ab665b27010a95303ac464c8880cec0a690bae0 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 22:00:03 -0400 Subject: [PATCH 14/29] update topmost CMakeLists --- CMakeLists.txt | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 02065c60..700119a9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,6 +61,30 @@ endif() find_package(ADIOS2 CONFIG 2.5 REQUIRED) find_package(Kokkos CONFIG 4.5 REQUIRED) +option(PCMS_ENABLE_MESHFIELDS "Enable MeshFields support" OFF) +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() + set(PETSC_INCLUDE_DIRS "${PETSC_DIR}/include" "${PETSC_DIR}/${PETSC_ARCH}/include") + set(PETSC_LIB_DIR "${PETSC_DIR}/${PETSC_ARCH}/lib") + find_library(PETSC_LIBRARY NAMES petsc HINTS ${PETSC_LIB_DIR} REQUIRED) + + 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) include(CTest) From 7a7ea7b550771211561ace0a103ce016d49270a3 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 22:06:03 -0400 Subject: [PATCH 15/29] update CMakeLists --- src/pcms/interpolator/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/interpolator/CMakeLists.txt b/src/pcms/interpolator/CMakeLists.txt index c8cf349a..d7bd0bb1 100644 --- a/src/pcms/interpolator/CMakeLists.txt +++ b/src/pcms/interpolator/CMakeLists.txt @@ -23,7 +23,7 @@ set_target_properties(pcms_interpolator PROPERTIES OUTPUT_NAME pcmsinterpolator EXPORT_NAME interpolator) add_library(pcms::interpolator ALIAS pcms_interpolator) -target_link_libraries(pcms_interpolator PUBLIC pcms::core Kokkos::kokkoskernels pcms_mesh_intersection) +target_link_libraries(pcms_interpolator PUBLIC pcms::core Kokkos::kokkoskernels) install( TARGETS pcms_interpolator From 7cce586941718e254995f18744c4dd44239a6793 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 22:15:03 -0400 Subject: [PATCH 16/29] add pcms namespace --- .../interpolator/mesh_intersection/calculate_load_vector.hpp | 5 ++++- .../interpolator/mesh_intersection/calculate_mass_matrix.hpp | 5 ++++- .../mesh_intersection/conservative_projection_solver.hpp | 4 ++++ .../mesh_intersection/load_vector_integrator.hpp | 4 +++- .../mesh_intersection/mass_matrix_integrator.hpp | 4 +++- 5 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp index 8077e85f..56c4a954 100644 --- a/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp @@ -73,6 +73,9 @@ using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; * @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, @@ -122,5 +125,5 @@ inline PetscErrorCode calculateLoadVector( *loadVec_out = vec; PetscFunctionReturn(PETSC_SUCCESS); } - +} // namespace pcms #endif diff --git a/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp index 96a1202e..29f7477d 100644 --- a/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp @@ -112,6 +112,9 @@ static PetscErrorCode CreateMatrix(Omega_h::Mesh& mesh, Mat* A) * @param[out] mass_out Pointer to the resulting mass matrix * @return PetscErrorCode PETSc error code (PETSC_SUCCESS if successful) */ + +namespace pcms +{ inline PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out) { PetscFunctionBeginUser; @@ -149,5 +152,5 @@ inline PetscErrorCode calculateMassMatrix(Omega_h::Mesh& mesh, Mat* mass_out) *mass_out = mass; PetscFunctionReturn(PETSC_SUCCESS); } - +} // namespace pcms #endif diff --git a/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp index 0c167d09..643e9393 100644 --- a/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp +++ b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp @@ -125,6 +125,9 @@ static Vec solveLinearSystem(Mat A, Vec b) * * */ + +namespace pcms +{ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, @@ -217,5 +220,6 @@ Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, return Omega_h::read(rhsvector); } +} // namespace pcms #endif // PCMS_INTERPOLATOR_GALERKIN_PROJECTION_SOLVER_HPP diff --git a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp index 2fd1a859..7a141dc8 100644 --- a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp @@ -277,6 +277,8 @@ struct IntegrationData * @see IntersectionResults */ +namespace pcms +{ Kokkos::View buildLoadVector( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values) @@ -375,7 +377,7 @@ Kokkos::View buildLoadVector( return elmLoadVector; } - +} // namespace pcms /// Holds projection and conservation error metrics returned by /// evaluate_pro_and_cons_errors(). struct Errors diff --git a/src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp b/src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp index 449504a1..d1fde063 100644 --- a/src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/mass_matrix_integrator.hpp @@ -10,6 +10,8 @@ #include #include +namespace pcms +{ // computes the mass matrix for each element template class MassMatrixIntegrator : public MeshField::Integrator @@ -77,5 +79,5 @@ Kokkos::View buildMassMatrix(Omega_h::Mesh& mesh, mmi.process(coordFe); return mmi.elmMassMatrix; } - +} // namespace pcms #endif From 9984c7bafc84409587b52325d531d28f3dab5cc2 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:14:25 -0400 Subject: [PATCH 17/29] add load vector test --- test/test_load_vector.cpp | 146 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 test/test_load_vector.cpp diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp new file mode 100644 index 00000000..358fa382 --- /dev/null +++ b/test/test_load_vector.cpp @@ -0,0 +1,146 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +TEST_CASE("Load vector computation on intersected regions", "[load_vector]") +{ + + Omega_h::Library lib; + + // 2D coordinates (x,y) : 4 vertices + Omega_h::Reals coords({ + 0.0, 0.0, // v0 + 1.0, 0.0, // v1 + 1.0, 1.0, // v2 + 0.0, 1.0 // v3 + }); + + // Target Mesh with two traingles + // Two triangles, CCW + // T0: (v0,v1,v3) = (0,1,3) + // T1: (v1,v2,v3) = (1,2,3) + Omega_h::LOs ev2v_source({0, 1, 3, 1, 2, 3}); + + Omega_h::Mesh target_mesh(&lib); + Omega_h::build_from_elems_and_coords(&target_mesh, OMEGA_H_SIMPLEX, 2, + ev2v_source, coords); + + // Source Mesh with two traingles + // Two triangles, CCW + // T0: (v0,v1,v3) = (0,1,2) + // T1: (v1,v2,v3) = (0,2,3) + Omega_h::LOs ev2v_target({0, 1, 2, 0, 2, 3}); + + Omega_h::Mesh source_mesh(&lib); + Omega_h::build_from_elems_and_coords(&source_mesh, OMEGA_H_SIMPLEX, 2, + ev2v_target, coords); + + REQUIRE(source_mesh.dim() == 2); + REQUIRE(target_mesh.dim() == 2); + + int num_tgt_elems = target_mesh.nelems(); + + auto intersection = intersectTargets(source_mesh, target_mesh); + SECTION("check localization routine for coincident cases") + { + Kokkos::View point("point", 3); + point(0, 0) = 0.0; + point(0, 1) = 0.0; + point(1, 0) = 1.0; + point(0, 1) = 0.0; + point(2, 0) = 1.0; + point(2, 1) = 1.0; + + pcms::GridPointSearch search_cell(target_mesh, 20, 20); + auto d_results = search_cell(point); + + auto h_results = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), d_results); + + REQUIRE(h_results.extent(0) == 3); + + for (int i = 0; i < h_results.extent(0); ++i) { + auto r = h_results(i); + std::cout << "tri_id = " << r.tri_id << "\n"; + std::cout << "dim = " << static_cast(r.dimensionality) << "\n"; + // std::cout << "parametric = " << r.parametric_coords << "\n"; + INFO("tri_id : " << r.tri_id); + REQUIRE(r.tri_id >= 0); + } + } + + SECTION("Basic shape and consistency of result") + { + + // Fill dummy values at each vertex of source + Omega_h::Write values(source_mesh.nverts(), 1.0); + + auto load_vector = + buildLoadVector(target_mesh, source_mesh, intersection, values); + + auto load_vector_host = Kokkos::create_mirror(load_vector); + Kokkos::deep_copy(load_vector_host, load_vector); + + REQUIRE(load_vector_host.extent(0) == num_tgt_elems * 3); + + for (int i = 0; i < load_vector_host.extent(0); ++i) + REQUIRE(load_vector_host(i) >= + 0.0); // Since function is 1.0 and everything is positive + } + + SECTION("Integration is zero if source values are zero") + { + Omega_h::Write zero_field(source_mesh.nverts(), 0.0); + + auto load_vector = + buildLoadVector(target_mesh, source_mesh, intersection, zero_field); + + auto load_vector_host = Kokkos::create_mirror_view(load_vector); + Kokkos::deep_copy(load_vector_host, load_vector); + + for (int i = 0; i < load_vector_host.extent(0); ++i) { + REQUIRE(load_vector_host(i) == Catch::Approx(0.0)); + } + } + + SECTION("load vector computed after the intersection of simple target and " + "source elements") + { + // the source elements are triangle1 (0,0), (1,0) & (1,1) and triangle2 + // (0,0), (1,1) & (0,1) the target elements are traingle1 (0,0), (1,0) & + // (0,1) and triangle2 (1,0), (1,1) & (0,1) + + Omega_h::Write constant_field(source_mesh.nverts(), 2.0); + + auto load_vector = + buildLoadVector(target_mesh, source_mesh, intersection, constant_field); + + auto load_vector_host = Kokkos::create_mirror_view(load_vector); + Kokkos::deep_copy(load_vector_host, load_vector); + + double expected_load_vector[6] = {0.333333, 0.333333, 0.333333, + 0.333333, 0.333333, 0.333333}; + double tolerance = 1e-6; + for (int i = 0; i < load_vector_host.extent(0); ++i) { + + CAPTURE(i, expected_load_vector[i], load_vector_host[i], tolerance); + CHECK_THAT(expected_load_vector[i], + Catch::Matchers::WithinAbs(load_vector_host[i], tolerance)); + } + } +} + +int main(int argc, char* argv[]) +{ + Kokkos::ScopeGuard guard(argc, argv); + Omega_h::Library lib(&argc, &argv); + return Catch::Session().run(argc, argv); +} From e67922eee44a06cb2a038448b71f533ba0f51a0e Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:15:32 -0400 Subject: [PATCH 18/29] add namespace --- .../interpolator/mesh_intersection/calculate_load_vector.hpp | 2 +- .../interpolator/mesh_intersection/calculate_mass_matrix.hpp | 2 +- .../mesh_intersection/conservative_projection_solver.hpp | 4 ++-- .../interpolator/mesh_intersection/load_vector_integrator.hpp | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp index 56c4a954..13ee05ac 100644 --- a/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_load_vector.hpp @@ -30,7 +30,7 @@ #include #include //ostringstream -#include "loadVectorIntegrator.hpp" +#include #include #include diff --git a/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp index 29f7477d..67a31205 100644 --- a/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp +++ b/src/pcms/interpolator/mesh_intersection/calculate_mass_matrix.hpp @@ -28,7 +28,7 @@ #include #include //ostringstream -#include "massMatrixIntegrator.hpp" +#include #include #include diff --git a/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp index 643e9393..37216796 100644 --- a/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp +++ b/src/pcms/interpolator/mesh_intersection/conservative_projection_solver.hpp @@ -31,8 +31,8 @@ #include #include -#include "calculateLoadVector.hpp" -#include "calculateMassMatrix.hpp" +#include +#include /** * @brief Solves a linear system Ax = b using PETSc's KSP solvers diff --git a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp index 7a141dc8..10456878 100644 --- a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp @@ -25,7 +25,7 @@ #include #include #include -#include "mesh_intersection.hpp" +#include #include /** From a9a2c579fd4678b83263b461d69883553be5410e Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:49:54 -0400 Subject: [PATCH 19/29] remove main --- test/test_load_vector.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp index 358fa382..6cf388f0 100644 --- a/test/test_load_vector.cpp +++ b/test/test_load_vector.cpp @@ -137,10 +137,3 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") } } } - -int main(int argc, char* argv[]) -{ - Kokkos::ScopeGuard guard(argc, argv); - Omega_h::Library lib(&argc, &argv); - return Catch::Session().run(argc, argv); -} From 1da8be902b42c2b10b749688e728ecb9f3f3066e Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:50:21 -0400 Subject: [PATCH 20/29] correct the name of petsc variables --- test/CMakeLists.txt | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index ee05be82..50a1e2b1 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -353,13 +353,32 @@ if(Catch2_FOUND) test_spline_interpolator.cpp test_svd_serial.cpp) endif() + + if (PCMS_ENABLE_MESHFIELDS) + list(APPEND PCMS_UNIT_TEST_SOURCES test_load_vector.cpp) + endif() + add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) - target_link_libraries(unit_tests PUBLIC Catch2::Catch2 pcms::core - pcms_interpolator) - target_include_directories(unit_tests PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) + + if (PCMS_ENABLE_MESHFIELDS) + target_link_libraries(unit_tests PUBLIC meshfields::meshfields) + endif() + + if (PCMS_ENABLE_PETSC) + target_include_directories(unit_tests PUBLIC ${PCMS_PETSC_INCLUDES}) + target_link_libraries(unit_tests PUBLIC ${PCMS_PETSC_LIBRARIES}) + endif() - include(Catch) - catch_discover_tests(unit_tests) + target_link_libraries(unit_tests PUBLIC + Catch2::Catch2 + pcms::core + pcms_interpolator + pcms_mesh_intersection + ) + + target_include_directories(unit_tests PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) + include(Catch) + Catch_discover_tests(unit_tests) else() message(WARNING "Catch2 not found. Disabling Unit Tests") endif() From 7c0fa98d4c10db84f226be980c1c1ad0eb198ac2 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:51:05 -0400 Subject: [PATCH 21/29] correct the petsc variables name --- src/pcms/interpolator/mesh_intersection/CMakeLists.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/CMakeLists.txt b/src/pcms/interpolator/mesh_intersection/CMakeLists.txt index e0bee537..e4bca7fd 100644 --- a/src/pcms/interpolator/mesh_intersection/CMakeLists.txt +++ b/src/pcms/interpolator/mesh_intersection/CMakeLists.txt @@ -24,8 +24,8 @@ if (PCMS_ENABLE_MESHFIELDS) endif() if (PCMS_ENABLE_PETSC) - target_include_directories(pcms_mesh_intersection PUBLIC ${PETSC_INCLUDES}) - target_link_libraries(pcms_mesh_intersection PUBLIC ${PETSC_LIBRARIES}) + 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) @@ -36,7 +36,6 @@ install(FILES ${MESH_INTERSECTION_HEADERS} install(TARGETS pcms_mesh_intersection EXPORT meshIntersectionTargets INCLUDES DESTINATION include/pcms/interpolator/mesh_intersection - ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) install(EXPORT meshIntersectionTargets From e72a6aa8904b6c7e3e248a85397505ad26e83172 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 21 Oct 2025 23:52:01 -0400 Subject: [PATCH 22/29] add PCMS_ENABLE_PETSC and PCMS_ENABLE_MESHFIELDS macros --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 700119a9..7115db51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -82,6 +82,7 @@ if (PCMS_ENABLE_PETSC) endif() set(PCMS_PETSC_INCLUDES ${PETSC_INCLUDE_DIRS} CACHE INTERNAL "") + set(PCMS_PETSC_LIBRARIES ${PETSC_LIBRARY} CACHE INTERNAL "") endif() From f90f795a7447b2f7e1ba22008e8faa4ba1753b31 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 13:50:20 -0400 Subject: [PATCH 23/29] fix memory space issue --- test/test_load_vector.cpp | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp index 6cf388f0..f9c560e0 100644 --- a/test/test_load_vector.cpp +++ b/test/test_load_vector.cpp @@ -51,19 +51,23 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") auto intersection = intersectTargets(source_mesh, target_mesh); SECTION("check localization routine for coincident cases") { - Kokkos::View point("point", 3); - point(0, 0) = 0.0; - point(0, 1) = 0.0; - point(1, 0) = 1.0; - point(0, 1) = 0.0; - point(2, 0) = 1.0; - point(2, 1) = 1.0; - + Kokkos::View points("test_points", 3); + auto points_h = Kokkos::create_mirror_view(points); + points_h(0, 0) = 0.0; + points_h(0, 1) = 0.0; + points_h(1, 0) = 1.0; + points_h(0, 1) = 0.0; + points_h(2, 0) = 1.0; + points_h(2, 1) = 1.0; + + Kokkos::deep_copy(points, points_h); pcms::GridPointSearch search_cell(target_mesh, 20, 20); - auto d_results = search_cell(point); - auto h_results = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), d_results); + auto d_results = search_cell(points); + + auto h_results = Kokkos::create_mirror_view(d_results); + + Kokkos::deep_copy(h_results, d_results); REQUIRE(h_results.extent(0) == 3); @@ -84,7 +88,7 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write values(source_mesh.nverts(), 1.0); auto load_vector = - buildLoadVector(target_mesh, source_mesh, intersection, values); + pcms::buildLoadVector(target_mesh, source_mesh, intersection, values); auto load_vector_host = Kokkos::create_mirror(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); @@ -101,7 +105,7 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write zero_field(source_mesh.nverts(), 0.0); auto load_vector = - buildLoadVector(target_mesh, source_mesh, intersection, zero_field); + pcms::buildLoadVector(target_mesh, source_mesh, intersection, zero_field); auto load_vector_host = Kokkos::create_mirror_view(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); @@ -120,8 +124,8 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write constant_field(source_mesh.nverts(), 2.0); - auto load_vector = - buildLoadVector(target_mesh, source_mesh, intersection, constant_field); + auto load_vector = pcms::buildLoadVector(target_mesh, source_mesh, + intersection, constant_field); auto load_vector_host = Kokkos::create_mirror_view(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); From cba1703bb3b8b94ce7ab1b5dbc1004f94e96f573 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 13:51:07 -0400 Subject: [PATCH 24/29] fix petsc variable name --- CMakeLists.txt | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7115db51..9347de39 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -73,9 +73,14 @@ 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") - find_library(PETSC_LIBRARY NAMES petsc HINTS ${PETSC_LIB_DIR} REQUIRED) + + 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}") From 7bfe489ad04b186e976f843f763cf3a3c476f1d6 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 13:51:42 -0400 Subject: [PATCH 25/29] fix header name --- .../interpolator/mesh_intersection/load_vector_integrator.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp index 10456878..e8ec96de 100644 --- a/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp +++ b/src/pcms/interpolator/mesh_intersection/load_vector_integrator.hpp @@ -516,7 +516,7 @@ Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, Kokkos::fabs(Omega_h::triangle_area_from_basis(basis)); sum_area += area; - const double EPS_AREA = abs_tol + rel_tol * ploy_area; + const double EPS_AREA = abs_tol + rel_tol * poly_area; if (area <= EPS_AREA) continue; // drops duplicates and colinear/degenerates From dd105ff50927f6caca56bd0a667015e0fb6d5161 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 13:52:12 -0400 Subject: [PATCH 26/29] fix header name --- src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp index c4d35c34..ca5f9eae 100644 --- a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp +++ b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp @@ -109,10 +109,6 @@ inline Kokkos::View compute_centroid(Omega_h::Mesh& mesh) * * @see r3d::intersect_simplices, intersectTargets */ -void FindIntersections::adjBasedIntersectSearch( - const Omega_h::LOs& tgt2src_offsets, - Omega_h::Write& nIntersections, - Omega_h::Write& tgt2src_indices, bool is_count_only); void FindIntersections::adjBasedIntersectSearch( const Omega_h::LOs& tgt2src_offsets, From 92dd379de69ee6c54498cfab1772ccd956b35148 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 14:41:01 -0400 Subject: [PATCH 27/29] add more rbfs and their use cases --- src/pcms/interpolator/mls_interpolation.cpp | 20 ++- src/pcms/interpolator/mls_interpolation.hpp | 144 +++++++++++++------- 2 files changed, 114 insertions(+), 50 deletions(-) diff --git a/src/pcms/interpolator/mls_interpolation.cpp b/src/pcms/interpolator/mls_interpolation.cpp index de59a67f..7cbfa795 100644 --- a/src/pcms/interpolator/mls_interpolation.cpp +++ b/src/pcms/interpolator/mls_interpolation.cpp @@ -121,10 +121,11 @@ struct RBF_MULTIQUADRIC double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho; if (rho_sq < r_sq) { phi = 0.0; } else { - phi = Kokkos::sqrt(1.0 + a * a * r * r); + phi = Kokkos::sqrt(1.0 + a * a * ratio * ratio); } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -147,10 +148,11 @@ struct RBF_INVMULTIQUADRIC double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho; if (rho_sq < r_sq) { phi = 0.0; } else { - phi = 1.0 / Kokkos::sqrt(1.0 + a * a * r * r); + phi = 1.0 / Kokkos::sqrt(1.0 + a * a * ratio * ratio); } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -164,17 +166,22 @@ struct RBF_INVMULTIQUADRIC // thin plate spline struct RBF_THINPLATESPLINE { + double a; + RBF_THINPLATESPLINE(double a_val) : a(a_val) {} + OMEGA_H_INLINE double operator()(double r_sq, double rho_sq) const { double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho; if (rho_sq < r_sq) { phi = 0.0; } else { - phi = r * r * Kokkos::log(r); + phi = + a * a * ratio * ratio * Kokkos::log(a * ratio); // a shouldnot be zero } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), @@ -188,16 +195,21 @@ struct RBF_THINPLATESPLINE // cubic struct RBF_CUBIC { + + double a; + RBF_INVMULTIQUADRIC(double a_val) : a(a_val) {} + OMEGA_H_INLINE double operator()(double r_sq, double rho_sq) const { double phi; double r = Kokkos::sqrt(r_sq); double rho = Kokkos::sqrt(rho_sq); + double ratio = r / rho; if (rho_sq < r_sq) { phi = 0.0; } else { - phi = r * r * r; + phi = a * a * a * ratio * ratio * ratio; } OMEGA_H_CHECK_PRINTF(!std::isnan(phi), diff --git a/src/pcms/interpolator/mls_interpolation.hpp b/src/pcms/interpolator/mls_interpolation.hpp index dfb387b1..b6b22245 100644 --- a/src/pcms/interpolator/mls_interpolation.hpp +++ b/src/pcms/interpolator/mls_interpolation.hpp @@ -5,40 +5,79 @@ namespace pcms { - /** - * @brief Enumeration of supported radial basis functions (RBF) for MLS - *interpolation. + * @brief Enumeration of supported radial basis functions (RBFs) for MLS + * interpolation. + * + * Specifies the kernel type used for computing weights in Moving Least Squares + * (MLS) interpolation. Each RBF defines a different spatial weighting behavior, + * influencing locality, smoothness, and conditioning of the interpolant. + * + * ### Values: + * + * - **RBF_GAUSSIAN (0)** + * \f$ \phi(r) = \exp\left(-\left(\frac{a r}{r_c}\right)^2\right) \f$ + * Infinitely smooth and compactly supported. + * `a` is the shape (decay) parameter — smaller values give wider influence. + * Good for smooth interpolants; may cause ill-conditioning for large `a`. + * + * - **RBF_C4 (1)** + * Compactly supported, \f$ C^4\f$-continuous polynomial kernel. + * Does not require a shape parameter. + * Efficient, stable, and good for local approximations. + * + * - **RBF_CONST (2)** + * Uniform weights within support: + * \f$ \phi(r) = 1 \text{ for } r < r_c \f$, otherwise 0. + * Equivalent to averaging. * - * This enum specifies the type of radial basis function used to weight source - *points in the Moving Least Squares (MLS) interpolation process. + * - **NO_OP (3)** + * Disables internal RBF weighting logic. + * Always returns \f$ \phi = 1 \f$ regardless of distance. + * Useful for methods like SPR patching * - * Members: - * - RBF_GAUSSIAN: - * A Gaussian RBF: `exp(- a^2 * r^2)`, smooth and commonly used. Good for - *localized support. `a` is a spreading/decay factor + * - **RBF_MULTIQUADRIC (4)** + * \f$ \phi(r) = \sqrt{1 + \left(\frac{a r}{r_c}\right)^2} \f$ + * Smooth and locally supported. + * `a` controls the width; higher values increase conditioning issues. * - * - RBF_C4: - * A compactly supported C4-continuous function. Useful for bounded support - *and efficiency. + * - **RBF_INVMULTIQUADRIC (5)** + * \f$ \phi(r) = \frac{1}{\sqrt{1 + \left(\frac{a r}{r_c}\right)^2}} \f$ + * Compactly supported and generally more stable than MQ. + * `a` adjusts the decay; larger `a` leads to faster fall-off. * - * - RBF_CONST: - * Constant basis function. Effectively uses uniform weights. + * - **RBF_THINPLATESPLINE (6)** + * \f$ \phi(r) = \left(\frac{a r}{r_c}\right)^2 \log\left(\frac{a + * r}{r_c}\right) \f$ Classic 2D surface fitting kernel. No explicit shape + * parameter is required, but `a` can scale smoothness. * - * - NO_OP: - * No operation. Disables RBF weighting — typically used when weights - * are externally defined or not needed. + * - **RBF_CUBIC (7)** + * \f$ \phi(r) = \left(\frac{a r}{r_c}\right)^3 \f$ + * Smooth with compact support. + * Simple and effective for local interpolation. * - * @note These are intended to be passed into function `mls_interpolation` to - *control weighting behavior. + * --- + * + * @note + * - The shape parameter `a` (also called decay or spread factor) is required + * for Gaussian, MQ, IMQ, ThinPlateSpline, and Cubic kernels. + * - All kernels are compactly supported within the cutoff radius \f$ r_c \f$. + * - A good starting value for `a` is 3–5. Avoid `a = 0` (leads to constant + * behavior). + * - In cases of ill-conditioning, use regularization (`lambda > 0`). + * + * @see mls_interpolation() */ enum class RadialBasisFunction : LO { RBF_GAUSSIAN = 0, RBF_C4, RBF_CONST, - NO_OP - + NO_OP, + RBF_MULTIQUADRIC, + RBF_INVMULTIQUADRIC, + RBF_THINPLATESPLINE, + RBF_CUBIC }; /** @@ -47,36 +86,49 @@ enum class RadialBasisFunction : LO * This function computes interpolated values at a set of target coordinates * using Moving Least Squares (MLS) based on the provided source values and * coordinates. It supports different radial basis functions (RBFs), polynomial - * degrees, dimension, optional regularization and tolerance + * degrees, spatial dimensions, optional regularization, and solver tolerance. * * @param source_values A flat array of source data values. Length should - * be `num_sources`. - * @param source_coordinates A flat array of source point coordinates. Length - * should be `num_sources * dim`. - * @param target_coordinates A flat array of target point coordinates. Length - * should be `num_targets * dim`. - * @param support A data structure holding neighbor information for - * each target (in - * CSR format). - * @param dim Dimension - * @param degree Polynomial degree - * @param bf The radial basis function used for weighting - * (e.g., Gaussian, C4). - * @param lambda Optional regularization parameter (default is - * 0.0). Helps with stability in ill-conditioned systems. - * @param tol Optional solver tolerance (default is 1e-6). - * Small singular values below this are discarded. - * - * @return A Write array containing the interpolated values at each target - * point. + * match the number of source points + * (`num_sources`). + * @param source_coordinates A flat array of source point coordinates of size + * `num_sources * dim`. + * @param target_coordinates A flat array of target point coordinates of size + * `num_targets * dim`. + * @param support A structure containing neighbor information for + * each target point in Compressed Sparse Row (CSR) + * format. + * @param dim Spatial dimension of the coordinate data (e.g., 2 + * or 3). + * @param degree Degree of the polynomial basis used in MLS. + * @param bf The radial basis function (RBF) used for + * weighting (e.g., Gaussian, C4). + * @param lambda Tikhonov regularization parameter (default: 0.0). + * A small positive value improves numerical + * stability for ill-conditioned local systems. + * - Well-conditioned systems: 1e-8 to 1e-6 + * - Mildly ill-conditioned: 1e-5 to 1e-3 + * - Highly ill-conditioned: 1e-2 to 1e-1 + * @param tol Optional solver tolerance (default: 1e-6). + * Singular values below this threshold are ignored. + * @param decay_factor Controls the spread of the RBF influence width + * (default: 5.0). + * - 1–3 → wide influence (smooth interpolation) + * - 3–8 → moderate locality (balanced, + * recommended) + * - 8–20 → highly local behavior (sharp detail, + * possible conditioning issues) + * + * @return A Write array containing interpolated values at each target + * point. Length matches the number of target points. * * @note - * - All input arrays are expected to reside in device memory (e.g., Kokkos - * device views). - * - Ensure consistency in dimensions: coordinate arrays must be sized as - * `num_points * dim`. - * - The result array length will match the number of target points. + * - All input arrays must reside in device memory (e.g., Kokkos device views). + * - Coordinate arrays must be correctly sized as `num_points * dim`. + * - Choose `lambda` and `decay_factor` carefully to balance accuracy and + * stability: a larger decay factor often requires a slightly higher `lambda`. */ + Omega_h::Write mls_interpolation( const Omega_h::Reals source_values, const Omega_h::Reals source_coordinates, const Omega_h::Reals target_coordinates, const SupportResults& support, From e0cf3572718f14e4a5efa2e8bbbc8a0d830e22d8 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 14:53:44 -0400 Subject: [PATCH 28/29] fix constructor --- src/pcms/interpolator/mls_interpolation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/interpolator/mls_interpolation.cpp b/src/pcms/interpolator/mls_interpolation.cpp index 7cbfa795..bc331143 100644 --- a/src/pcms/interpolator/mls_interpolation.cpp +++ b/src/pcms/interpolator/mls_interpolation.cpp @@ -197,7 +197,7 @@ struct RBF_CUBIC { double a; - RBF_INVMULTIQUADRIC(double a_val) : a(a_val) {} + RBF_CUBIC(double a_val) : a(a_val) {} OMEGA_H_INLINE double operator()(double r_sq, double rho_sq) const From ba5d9ab14bfa8a3b5d7cb0cff50a0a7d0d29af2c Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 22 Oct 2025 15:13:17 -0400 Subject: [PATCH 29/29] remove semicolon --- src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp index ca5f9eae..fab31873 100644 --- a/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp +++ b/src/pcms/interpolator/mesh_intersection/mesh_intersection.hpp @@ -38,7 +38,9 @@ class FindIntersections public: FindIntersections(Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh) - : source_mesh_(source_mesh), target_mesh_(target_mesh) {}; + : source_mesh_(source_mesh), target_mesh_(target_mesh) + { + } void adjBasedIntersectSearch(const Omega_h::LOs& tgt2src_offsets, Omega_h::Write& nIntersections,