|
| 1 | +// |
| 2 | +// Created by hasanm4 on 2/17/25. |
| 3 | +// |
| 4 | +#include <pcms/capi/kokkos.h> |
| 5 | +#include <pcms/capi/interpolator.h> |
| 6 | +#include <pcms/interpolator/interpolation_base.h> |
| 7 | +#include <Omega_h_file.hpp> |
| 8 | +#include <Omega_h_library.hpp> |
| 9 | +#include <Omega_h_mesh.hpp> |
| 10 | +#include <pcms/print.h> |
| 11 | + |
| 12 | +//[[nodiscard]] |
| 13 | +PcmsInterpolatorHandle pcms_create_interpolator(PcmsOmegaHMeshHandle oh_mesh, |
| 14 | + double radius) |
| 15 | +{ |
| 16 | + auto* source_mesh = reinterpret_cast<Omega_h::Mesh*>(oh_mesh.mesh_handle); |
| 17 | + auto* interpolator = new MLSMeshInterpolation(*source_mesh, radius); |
| 18 | + return {reinterpret_cast<void*>(interpolator)}; |
| 19 | +} |
| 20 | + |
| 21 | +PcmsInterpolatorHandle pcms_create_point_based_interpolator( |
| 22 | + void* source_points, int source_points_size, void* target_points, |
| 23 | + int target_points_size, double radius, int degree, int min_req_supports, |
| 24 | + double lambda, double decay_factor) |
| 25 | +{ |
| 26 | + |
| 27 | + auto source_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>( |
| 28 | + reinterpret_cast<double*>(source_points), source_points_size); |
| 29 | + auto target_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>( |
| 30 | + reinterpret_cast<double*>(target_points), target_points_size); |
| 31 | + auto* interpolator = new MLSPointCloudInterpolation( |
| 32 | + source_points_view, target_points_view, 2, radius, min_req_supports, degree, |
| 33 | + true, lambda, decay_factor); |
| 34 | + return {reinterpret_cast<void*>(interpolator)}; |
| 35 | +} |
| 36 | + |
| 37 | +Omega_h::HostRead<Omega_h::Real> read_mesh_centroids(const char* mesh_filename, |
| 38 | + int& num_elements) |
| 39 | +{ |
| 40 | + auto fname = std::string(mesh_filename); |
| 41 | + fname = fname.erase(fname.find_last_not_of(" \n\r\t") + 1); |
| 42 | + pcms::printInfo("The interpolator got dg2 mesh file: %s\n", fname.c_str()); |
| 43 | + auto mesh_lib = Omega_h::Library(nullptr, nullptr, MPI_COMM_SELF); |
| 44 | + auto mesh = Omega_h::binary::read(fname, mesh_lib.world()); |
| 45 | + auto elem_centroids = getCentroids(mesh); |
| 46 | + num_elements = mesh.nelems(); |
| 47 | + OMEGA_H_CHECK_PRINTF(num_elements * 2 == elem_centroids.size(), |
| 48 | + "Mesh element centroids size does not match the number " |
| 49 | + "of elements %d != %d\n", |
| 50 | + num_elements * 2, elem_centroids.size()); |
| 51 | + |
| 52 | + pcms::printInfo("Number of element centroids: %d\n", |
| 53 | + elem_centroids.size() / 2); |
| 54 | + OMEGA_H_CHECK_PRINTF(mesh.dim() == 2, "Mesh dimension is not 2D %d\n", |
| 55 | + mesh.dim()); |
| 56 | + |
| 57 | + return {elem_centroids}; |
| 58 | +} |
| 59 | + |
| 60 | +void write_void_int_pointer(void* pointer, int value) |
| 61 | +{ |
| 62 | + if (pointer) { |
| 63 | + int* dg2_elem_count_int = reinterpret_cast<int*>(pointer); |
| 64 | + *dg2_elem_count_int = value; |
| 65 | + } else { |
| 66 | + pcms::printError("Error: NULL pointer provided to write integer value\n"); |
| 67 | + } |
| 68 | +} |
| 69 | + |
| 70 | +PcmsInterpolatorHandle pcms_create_degas2xgcnode_interpolator( |
| 71 | + void* target_points, int target_points_size, const char* dg2_mesh_filename, |
| 72 | + double radius, void* dg2_elem_count, int degree, int min_req_supports, |
| 73 | + double lambda, double decay_factor) |
| 74 | +{ |
| 75 | + // same as above pcms_create_degas2xgc_interpolator but the target points are |
| 76 | + // provided by the user this is useful when the corresponding xgc mesh is not |
| 77 | + // available |
| 78 | + |
| 79 | + int dg2_num_elems = 0; |
| 80 | + Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host = |
| 81 | + read_mesh_centroids(dg2_mesh_filename, dg2_num_elems); |
| 82 | + write_void_int_pointer(dg2_elem_count, dg2_num_elems); |
| 83 | + |
| 84 | + return pcms_create_point_based_interpolator( |
| 85 | + (void*)dg2_elem_centroids_host.data(), dg2_elem_centroids_host.size(), |
| 86 | + target_points, target_points_size, radius, degree, min_req_supports, lambda, |
| 87 | + decay_factor); |
| 88 | +} |
| 89 | + |
| 90 | +PcmsInterpolatorHandle pcms_create_xgcnodedegas2_interpolator( |
| 91 | + const char* dg2_mesh_filename, void* source_points, int source_points_size, |
| 92 | + double radius, void* dg2_elem_count, int degree, int min_req_supports, |
| 93 | + double lambda, double decay_factor) |
| 94 | +{ |
| 95 | + int dg2_num_elems = 0; |
| 96 | + Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host = |
| 97 | + read_mesh_centroids(dg2_mesh_filename, dg2_num_elems); |
| 98 | + write_void_int_pointer(dg2_elem_count, dg2_num_elems); |
| 99 | + |
| 100 | + return pcms_create_point_based_interpolator( |
| 101 | + source_points, source_points_size, (void*)dg2_elem_centroids_host.data(), |
| 102 | + dg2_elem_centroids_host.size(), radius, degree, min_req_supports, lambda, |
| 103 | + decay_factor); |
| 104 | +} |
| 105 | + |
| 106 | +void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator) |
| 107 | +{ |
| 108 | + if (interpolator.pointer != nullptr) { |
| 109 | + delete reinterpret_cast<InterpolationBase*>(interpolator.pointer); |
| 110 | + } |
| 111 | +} |
| 112 | + |
| 113 | +void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, |
| 114 | + int input_size, void* output, int output_size) |
| 115 | +{ |
| 116 | + auto* mls_interpolator = |
| 117 | + reinterpret_cast<InterpolationBase*>(interpolator.pointer); |
| 118 | + |
| 119 | + OMEGA_H_CHECK_PRINTF( |
| 120 | + input_size == mls_interpolator->getSourceSize(), |
| 121 | + "Input array size does not match the source size %d != %zu\n", input_size, |
| 122 | + mls_interpolator->getSourceSize()); |
| 123 | + OMEGA_H_CHECK_PRINTF( |
| 124 | + output_size == mls_interpolator->getTargetSize(), |
| 125 | + "Output array size does not match the target size %d != %zu\n", output_size, |
| 126 | + mls_interpolator->getTargetSize()); |
| 127 | + |
| 128 | + pcms::Rank1View<double, pcms::HostMemorySpace> input_array( |
| 129 | + reinterpret_cast<double*>(input), input_size); |
| 130 | + pcms::Rank1View<double, pcms::HostMemorySpace> output_array( |
| 131 | + reinterpret_cast<double*>(output), output_size); |
| 132 | + |
| 133 | + mls_interpolator->eval(input_array, output_array); |
| 134 | +} |
0 commit comments