diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6052152..d71f9c5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,6 +5,7 @@ FIND_PACKAGE(IBAMR REQUIRED)
 
 SET(CXX_SRC
     CohesionStressRHS.cpp
+    utility_functions.cpp
     )
 ADD_EXECUTABLE(main2d main.cpp)
 TARGET_SOURCES(main2d PRIVATE ${CXX_SRC})
diff --git a/Makefile b/Makefile
index cc44d1d..bb347e7 100644
--- a/Makefile
+++ b/Makefile
@@ -13,7 +13,7 @@ include $(IBAMR_BUILD_DIR)/config/make.inc
 ## main driver is in main.cpp
 ##
 ## PDIM = 2 implies two spatial dimensions
-OBJS = CohesionRHS.o
+OBJS = CohesionRHS.o utility_functions.o
 #main.o
 PDIM = 2
 
diff --git a/utility_functions.cpp b/utility_functions.cpp
new file mode 100644
index 0000000..2fe0147
--- /dev/null
+++ b/utility_functions.cpp
@@ -0,0 +1,42 @@
+#include <ibamr/app_namespaces.h>
+
+#include <ibtk/ibtk_utilities.h>
+
+#include <Box.h>
+
+// Local includes
+#include "utility_functions.h"
+
+namespace IBAMR
+{
+double
+convolution(double alpha,
+            const CellData<NDIM, double>& a_data,
+            double beta,
+            const CellData<NDIM, double>& b_data,
+            std::function<double(double)> phi,
+            unsigned int phi_width,
+            const CellIndex<NDIM>& idx,
+            const double* const dx)
+{
+    // Box over which integration is performed
+    Box<NDIM> int_box(idx, idx);
+    int_box.grow(phi_width);
+#ifndef NDEBUG
+    // Error checking. Need ghost data for patch data.
+    TBOX_ASSERT(a_data.getGhostBox() * int_box == int_box);
+    TBOX_ASSERT(b_data.getGhostBox() * int_box == int_box);
+#endif
+    // Now compute convolution
+    double convolve = 0.0;
+    for (CellIterator<NDIM> ci(int_box); ci; ci++)
+    {
+        const CellIndex<NDIM>& n_idx = ci();
+        double phi_weight = 1.0;
+        ;
+        for (int d = 0; d < NDIM; ++d) phi_weight *= phi(static_cast<double>(n_idx(d) - idx(d))) * dx[d];
+        convolve += (alpha * a_data(n_idx) + beta * b_data(n_idx)) * phi_weight;
+    }
+    return convolve;
+}
+} // namespace IBAMR
diff --git a/utility_functions.h b/utility_functions.h
new file mode 100644
index 0000000..6f6210d
--- /dev/null
+++ b/utility_functions.h
@@ -0,0 +1,32 @@
+/////////////////////////////// INCLUDE GUARD ////////////////////////////////
+
+#ifndef included_utility_functions
+#define included_utility_functions
+
+/////////////////////////////// INCLUDES /////////////////////////////////////
+#include <ibamr/config.h>
+
+#include <CellData.h>
+#include <CellIndex.h>
+
+#include <functional>
+namespace IBAMR
+{
+/*!
+ * Computes the convolution of (alpha*a + beta*b) and the function phi at a given patch index, assuming phi has finite
+ * support.
+ *
+ * Both patch data MUST have enough ghost cells for the width of the phi function. This assumes the phi function is
+ * centered at idx and scaled by the grid spacing.
+ */
+double convolution(double alpha,
+                   const SAMRAI::pdat::CellData<NDIM, double>& a_data,
+                   double beta,
+                   const SAMRAI::pdat::CellData<NDIM, double>& b_data,
+                   std::function<double(double)> phi,
+                   unsigned int phi_width,
+                   const SAMRAI::pdat::CellIndex<NDIM>& idx,
+                   const double* const dx);
+} // namespace IBAMR
+
+#endif // included_utility_functions