diff --git a/src/pcms/interpolator/CMakeLists.txt b/src/pcms/interpolator/CMakeLists.txt
index 1e17e6c5..b60df52c 100644
--- a/src/pcms/interpolator/CMakeLists.txt
+++ b/src/pcms/interpolator/CMakeLists.txt
@@ -12,10 +12,12 @@ set(PCMS_FIELD_TRANSFER_HEADERS
     linear_interpolant.hpp
     multidimarray.hpp
     MLS_rbf_options.hpp
+        interpolation_base.h
 )
 
 set(PCMS_FIELD_TRANSFER_SOURCES
     MLS_rbf_options.cpp
+        interpolation_base.cpp
 )
 install(FILES ${PCMS_FIELD_TRANSFER_HEADERS} DESTINATION include/pcms/interpolator)
 
diff --git a/src/pcms/interpolator/interpolation_base.cpp b/src/pcms/interpolator/interpolation_base.cpp
new file mode 100644
index 00000000..8da33ed8
--- /dev/null
+++ b/src/pcms/interpolator/interpolation_base.cpp
@@ -0,0 +1,175 @@
+//
+// Created by hasanm4 on 1/17/25.
+//
+
+#include "interpolation_base.h"
+#include "MLS_rbf_options.hpp"
+
+#include <execution>
+
+Omega_h::Reals getCentroids(Omega_h::Mesh& mesh)
+{
+  OMEGA_H_CHECK_PRINTF(
+    mesh.dim() == 2, "Only 2D meshes are supported but found %d\n", mesh.dim());
+
+  const auto& coords = mesh.coords();
+  Omega_h::Write<Omega_h::Real> centroids(mesh.nfaces() * mesh.dim(), 0.0);
+
+  auto face2node = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
+  Omega_h::parallel_for(
+    mesh.nfaces(), OMEGA_H_LAMBDA(Omega_h::LO face) {
+      auto nodes = Omega_h::gather_verts<3>(face2node, face);
+      Omega_h::Few<Omega_h::Vector<2>, 3> face_coords =
+        Omega_h::gather_vectors<3, 2>(coords, nodes);
+      Omega_h::Vector<2> centroid = Omega_h::average(face_coords);
+      centroids[2 * face + 0] = centroid[0];
+      centroids[2 * face + 1] = centroid[1];
+    });
+
+  return {centroids};
+}
+
+MLSInterpolationHandler::MLSInterpolationHandler(Omega_h::Mesh& source_mesh,
+                                                 double radius,
+                                                 uint min_req_support,
+                                                 uint degree, bool adapt_radius)
+  : source_mesh_(source_mesh),
+    target_mesh_(source_mesh),
+    radius_(radius),
+    min_req_supports_(min_req_support),
+    degree_(degree),
+    adapt_radius_(adapt_radius)
+{
+  single_mesh_ = true;
+  target_coords_ = source_mesh_.coords();
+  source_coords_ = getCentroids(source_mesh_);
+
+  OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2,
+                       "Only 2D meshes are supported but found %d\n",
+                       source_mesh_.dim());
+
+  source_field_ =
+    Omega_h::HostWrite<Omega_h::Real>(source_mesh_.nfaces(), "source field");
+
+  target_field_ =
+    Omega_h::HostWrite<Omega_h::Real>(source_mesh_.nverts(), "target field");
+
+  find_supports(min_req_supports_);
+}
+
+MLSInterpolationHandler::MLSInterpolationHandler(
+  Omega_h::Mesh& source_mesh, Omega_h::Mesh& target_mesh, const double radius,
+  uint min_req_support, uint degree, const bool adapt_radius)
+  : source_mesh_(source_mesh),
+    target_mesh_(target_mesh),
+    radius_(radius),
+    min_req_supports_(min_req_support),
+    degree_(degree),
+    adapt_radius_(adapt_radius)
+{
+  OMEGA_H_CHECK_PRINTF(source_mesh_.dim() == 2 && target_mesh_.dim() == 2,
+                       "Only 2D meshes are supported but found %d, %d\n",
+                       source_mesh_.dim(), target_mesh_.dim());
+
+  source_coords_ = source_mesh_.coords();
+  target_coords_ = target_mesh_.coords();
+
+  source_field_ =
+    Omega_h::HostWrite<Omega_h::Real>(source_mesh_.nverts(), "source field");
+  target_field_ =
+    Omega_h::HostWrite<Omega_h::Real>(target_mesh_.nverts(), "target field");
+
+  find_supports(min_req_supports_);
+}
+
+// TODO : find way to avoid this copy
+void copyHostScalarArrayView2HostWrite(
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> source,
+  Omega_h::HostWrite<Omega_h::Real>& target)
+{
+  OMEGA_H_CHECK_PRINTF(
+    source.size() == target.size(),
+    "Size mismatch in copy_data_from_ScalarArray_to_HostWrite: %zu %d\n",
+    source.size(), target.size());
+
+  for (int i = 0; i < source.size(); ++i) {
+    target[i] = source[i];
+  }
+}
+void copyHostWrite2ScalarArrayView(
+  const Omega_h::HostWrite<Omega_h::Real>& source,
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> target)
+{
+  OMEGA_H_CHECK_PRINTF(
+    source.size() == target.size(),
+    "Size mismatch in copy_data_from_HostWrite_to_ScalarArray: %d %zu\n",
+    source.size(), target.size());
+
+  for (int i = 0; i < source.size(); ++i) {
+    target[i] = source[i];
+  }
+}
+
+void MLSInterpolationHandler::eval(
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> source_field,
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> target_field)
+{
+  OMEGA_H_CHECK_PRINTF(
+    target_field.size() == target_coords_.size() / target_mesh_.dim(),
+    "Source Data and Source Points size mismatch: %zu %d\n",
+    target_field.size(), target_coords_.size() / target_mesh_.dim());
+
+  OMEGA_H_CHECK_PRINTF(
+    source_field.size() == source_coords_.size() / source_mesh_.dim(),
+    "Target Data and Target Points size mismatch: %zu %d\n",
+    source_field.size(), source_coords_.size() / source_mesh_.dim());
+
+  copyHostScalarArrayView2HostWrite(source_field, source_field_);
+
+  // TODO: make the basis function a template or pass it as a parameter
+  auto target_field_write = mls_interpolation(
+    Omega_h::Reals(source_field_), source_coords_, target_coords_, supports_, 2,
+    degree_, supports_.radii2, RadialBasisFunction::RBF_GAUSSIAN);
+
+  target_field_ = Omega_h::HostWrite<Omega_h::Real>(target_field_write);
+  copyHostWrite2ScalarArrayView(target_field_, target_field);
+}
+
+void MLSInterpolationHandler::find_supports(const uint min_req_support)
+{
+  if (single_mesh_) {
+    supports_ =
+      searchNeighbors(source_mesh_, radius_, min_req_support, adapt_radius_);
+  } else { // two mesh : vert to vert
+    supports_ = searchNeighbors(source_mesh_, target_mesh_, radius_,
+                                min_req_support, adapt_radius_);
+  }
+
+#ifndef NDEBUG
+  Omega_h::HostRead<Omega_h::Real> hostRadii2(supports_.radii2);
+  for (size_t i = 0; i < hostRadii2.size(); ++i) {
+    OMEGA_H_CHECK_PRINTF(
+      hostRadii2[i] > 1e-10,
+      "Radius squared has to be more than zero found found [%zu] = %f\n", i,
+      hostRadii2[i]);
+  }
+#endif
+}
+
+size_t MLSInterpolationHandler::getSourceSize()
+{
+  if (single_mesh_) {
+    return source_mesh_.nfaces();
+  } else {
+    return source_mesh_.nverts();
+  }
+}
+
+size_t MLSInterpolationHandler::getTargetSize()
+{
+  if (single_mesh_) {
+    return source_mesh_.nverts();
+  } else {
+    return target_mesh_.nverts();
+  }
+}
\ No newline at end of file
diff --git a/src/pcms/interpolator/interpolation_base.h b/src/pcms/interpolator/interpolation_base.h
new file mode 100644
index 00000000..deddecf2
--- /dev/null
+++ b/src/pcms/interpolator/interpolation_base.h
@@ -0,0 +1,88 @@
+//
+// Created by Fuad Hasan on 1/12/25.
+//
+#include "MLSInterpolation.hpp"
+#include "adj_search.hpp"
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <pcms/arrays.h>
+#include <string>
+
+#ifndef PCMS_INTERPOLATION_BASE_H
+#define PCMS_INTERPOLATION_BASE_H
+
+class InterpolationBase
+{
+  /**
+   * @brief Evaluate the interpolation
+   * @param source_field The field to interpolate from
+   * @param target_field The field to interpolate to
+   */
+  virtual void eval(
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> source_field,
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> target_field) = 0;
+};
+
+/**
+ * @brief Moving Least Square Radial Basis Function Interpolation
+ */
+class MLSInterpolationHandler : public InterpolationBase
+{
+
+public:
+  void eval(
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> source_field,
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> target_field) override;
+
+  /**
+   * @brief Vertex to Vertex interpolation for two given meshes
+   * @param source_mesh The source mesh
+   * @param target_mesh The target mesh
+   * @param radius The cutoff radius for the MLS interpolation
+   * @param min_req_supports The minimum number of source locations required for
+   * interpolation
+   * @param degree The degree of the polynomial used in the MLS interpolation
+   * @param adapt_radius Whether to adapt the radius based on the local density
+   */
+  MLSInterpolationHandler(Omega_h::Mesh& source_mesh,
+                          Omega_h::Mesh& target_mesh, double radius,
+                          uint min_req_supports = 10, uint degree = 3,
+                          bool adapt_radius = true);
+
+  /**
+   * @brief Centroids to Vertices interpolation for a single mesh
+   * @param source_mesh The source mesh
+   * @param radius The cutoff radius for the MLS interpolation
+   * @param adapt_radius Whether to adapt the radius based on the local density
+   */
+  MLSInterpolationHandler(Omega_h::Mesh& source_mesh, double radius,
+                          uint min_req_supports = 10, uint degree = 3,
+                          bool adapt_radius = true);
+
+  size_t getSourceSize();
+  size_t getTargetSize();
+
+private:
+  double radius_;
+  bool adapt_radius_;
+  bool single_mesh_ = false;
+  uint degree_;
+  uint min_req_supports_;
+
+  // InterpolationType interpolation_type_;
+
+  Omega_h::Mesh& source_mesh_;
+  // TODO: handle what to do with this when only 1 mesh is provided
+  Omega_h::Mesh& target_mesh_;
+  Omega_h::Reals source_coords_;
+  Omega_h::Reals target_coords_;
+
+  SupportResults supports_;
+
+  Omega_h::HostWrite<Omega_h::Real> target_field_;
+  Omega_h::HostWrite<Omega_h::Real> source_field_;
+
+  void find_supports(uint min_req_supports = 10);
+};
+
+#endif // PCMS_INTERPOLATION_BASE_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e4d8e4e7..b73a9f30 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -121,6 +121,7 @@ if(Catch2_FOUND)
 	      test_mls_basis.cpp
 	      test_rbf_interp.cpp
 	      test_linear_solver.cpp
+              test_interpolation_class.cpp
               )
   endif ()
   add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES})
diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp
new file mode 100644
index 00000000..3d0a3e7c
--- /dev/null
+++ b/test/test_interpolation_class.cpp
@@ -0,0 +1,190 @@
+//
+// Created by hasanm4 on 1/17/25.
+//
+
+#include <catch2/catch_test_macros.hpp>
+
+#include <Omega_h_build.hpp>
+#include <Omega_h_library.hpp>
+#include <pcms/interpolator/interpolation_base.h>
+
+#include <vector>
+#include <iostream>
+
+void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector) {
+  auto coords = mesh->coords();
+  auto nverts = mesh->nverts();
+  auto out = Write<Real>(coords.size());
+
+  auto f = OMEGA_H_LAMBDA(LO i) {
+    auto coord = get_vector<2>(coords, i);
+    coord = coord + translation_vector;
+    set_vector<2>(out, i, coord);
+  };
+
+  parallel_for(nverts, f);
+  mesh->set_coords(Reals(out));
+}
+
+
+bool isClose(Omega_h::HostWrite<Omega_h::Real>& array1,
+             Omega_h::HostWrite<Omega_h::Real>& array2,
+             double percent_diff = 0.1);
+
+/**
+ * @brief Create sin(x)cos(y) + 2 at each vertex of the mesh
+ */
+void createSinxCosyAtVertex(Omega_h::Mesh& mesh,
+                            Omega_h::Write<Omega_h::Real>& sinxcosy)
+{
+  OMEGA_H_CHECK_PRINTF(sinxcosy.size() == mesh.nverts(),
+                       "Given ScalarArrayView is not properly sized %d vs %d\n",
+                       mesh.nverts(), sinxcosy.size());
+
+  // get the bounding box of the mesh
+  Omega_h::BBox<2> bb = Omega_h::get_bounding_box<2>(&mesh);
+  double dx = bb.max[0] - bb.min[0];
+  double dy = bb.max[1] - bb.min[1];
+
+  const auto& coords = mesh.coords();
+
+  auto assignSinCos = OMEGA_H_LAMBDA(int node)
+  {
+    auto x = (coords[2 * node + 0] - bb.min[0]) / dx * 2.0 * M_PI;
+    auto y = (coords[2 * node + 1] - bb.min[1]) / dy * 2.0 * M_PI;
+    sinxcosy[node] = 2 + Kokkos::sin(x) * Kokkos::cos(y);
+  };
+  Omega_h::parallel_for(mesh.nverts(), assignSinCos, "assignSinCos");
+}
+
+void node2CentroidInterpolation(Omega_h::Mesh& mesh,
+                                Omega_h::Write<Omega_h::Real>& node_vals,
+                                Omega_h::Write<Omega_h::Real>& centroid_vals)
+{
+  OMEGA_H_CHECK_PRINTF(node_vals.size() == mesh.nverts(),
+                       "Given ScalarArrayView is not properly sized %d vs %d\n",
+                       mesh.nverts(), node_vals.size());
+  OMEGA_H_CHECK_PRINTF(centroid_vals.size() == mesh.nfaces(),
+                       "Given ScalarArrayView is not properly sized %d vs %d\n",
+                       mesh.nfaces(), centroid_vals.size());
+
+  const auto& face2node = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b;
+  auto coords = mesh.coords();
+  Omega_h::LO nfaces = mesh.nfaces();
+
+  auto averageSinCos = OMEGA_H_LAMBDA(Omega_h::LO face)
+  {
+    auto faceNodes = Omega_h::gather_verts<3>(face2node, face);
+    Omega_h::Real sum = node_vals[faceNodes[0]] + node_vals[faceNodes[1]] +
+                        node_vals[faceNodes[2]];
+    centroid_vals[face] = sum / 3.0;
+  };
+  Omega_h::parallel_for(nfaces, averageSinCos, "averageSinCos");
+}
+
+TEST_CASE("Test MLSInterpolationHandler")
+{
+  fprintf(stdout, "[INFO] Starting MLS Interpolation Test...\n");
+  auto lib = Library{};
+  auto world = lib.world();
+  auto source_mesh =
+    Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false);
+  printf("[INFO] Mesh created with %d vertices and %d faces\n",
+         source_mesh.nverts(), source_mesh.nfaces());
+
+  Omega_h::Write<Omega_h::Real> source_sinxcosy_node(source_mesh.nverts(),
+                                                     "source_sinxcosy_node");
+  createSinxCosyAtVertex(source_mesh, source_sinxcosy_node);
+
+  SECTION("Single Mesh")
+  {
+    fprintf(stdout, "\n-------------------- Single Mesh Interpolation Test Started --------------------\n");
+    auto mls_single = MLSInterpolationHandler(source_mesh, 0.12, 12, 3, true);
+
+    Omega_h::Write<Omega_h::Real> sinxcosy_centroid(source_mesh.nfaces(),
+                                                    "sinxcosy_centroid");
+    node2CentroidInterpolation(source_mesh, source_sinxcosy_node,
+                               sinxcosy_centroid);
+
+    Omega_h::HostWrite<double> source_data_host_write(sinxcosy_centroid);
+    Omega_h::HostWrite<double> interpolated_data_hwrite(source_mesh.nverts());
+    Omega_h::HostWrite<double> exact_values_at_nodes(source_sinxcosy_node);
+
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> sourceArrayView(
+      source_data_host_write.data(), source_data_host_write.size());
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> interpolatedArrayView(
+      interpolated_data_hwrite.data(), interpolated_data_hwrite.size());
+
+    OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(),
+                         "Source size mismatch: %zu vs %zu\n",
+                         sourceArrayView.size(), mls_single.getSourceSize());
+    OMEGA_H_CHECK_PRINTF(
+      interpolatedArrayView.size() == mls_single.getTargetSize(),
+      "Target size mismatch: %zu vs %zu\n", interpolatedArrayView.size(),
+      mls_single.getTargetSize());
+
+    mls_single.eval(sourceArrayView, interpolatedArrayView);
+
+    REQUIRE(isClose(exact_values_at_nodes, interpolated_data_hwrite, 10.0) ==
+            true);
+    fprintf(stdout, "[****] Single Mesh Interpolation Test Passed with %.2f%% tolerance!\n", 10.0);
+  }
+
+  SECTION("Double Mesh")
+  {
+    fprintf(stdout, "\n-------------------- Double Mesh Interpolation Test Started --------------------\n");
+    auto target_mesh =
+      Omega_h::build_box(world, OMEGA_H_SIMPLEX, 0.999, 0.999, 1, 17, 17, 0, false);
+    printf("[INFO] Target Mesh created with %d vertices and %d faces\n",
+           target_mesh.nverts(), target_mesh.nfaces());
+
+    // TODO: This is a way around. https://github.com/SCOREC/pcms/pull/148#discussion_r1926204199
+    translate_mesh(&target_mesh, Omega_h::Vector<2>{(1.0-0.999)/2.0, (1.0-0.999)/2.0});
+
+    auto mls_double =
+      MLSInterpolationHandler(source_mesh, target_mesh, 0.12, 12, 3, true);
+
+    Omega_h::HostWrite<double> source_data_host_write(source_sinxcosy_node);
+    Omega_h::HostWrite<double> interpolated_data_hwrite(
+      mls_double.getTargetSize());
+
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> sourceArrayView(
+      source_data_host_write.data(), source_data_host_write.size());
+    pcms::ScalarArrayView<double, pcms::HostMemorySpace> interpolatedArrayView(
+      interpolated_data_hwrite.data(), interpolated_data_hwrite.size());
+
+    mls_double.eval(sourceArrayView, interpolatedArrayView);
+
+    Omega_h::Write<Omega_h::Real> exact_target_sinxcosy_node(
+      target_mesh.nverts(), "target_sinxcosy_node");
+    createSinxCosyAtVertex(target_mesh, exact_target_sinxcosy_node);
+    Omega_h::HostWrite<double> exact_target_sinxcosy_node_hwrite(
+      exact_target_sinxcosy_node);
+
+    REQUIRE(isClose(interpolated_data_hwrite, exact_target_sinxcosy_node_hwrite,
+                    10.0) == true);
+
+    fprintf(stdout, "[INFO] Double Mesh Interpolation Test Passed with %.2f%% tolerance!\n", 10.0);
+  }
+}
+
+bool isClose(Omega_h::HostWrite<Omega_h::Real>& array1,
+             Omega_h::HostWrite<Omega_h::Real>& array2, double percent_diff)
+{
+  if (array1.size() != array2.size()) {
+    fprintf(stderr, "[ERROR] Arrays are not of the same size: %d vs %d\n",
+            array1.size(), array2.size());
+    return false;
+  }
+
+  double eps = percent_diff / 100.0;
+  for (int i = 0; i < array1.size(); i++) {
+    if (std::abs(array1[i] - array2[i]) > eps * std::abs(array2[i])) {
+      fprintf(stderr, "[ERROR] Arrays differ at index %d: %.16f vs %.16f\n", i,
+              array1[i], array2[i]);
+      return false;
+    }
+  }
+
+  return true;
+}
\ No newline at end of file