diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e573adc0d8c594c404503a915154c0b5c09ae076..d88a00fb095a9ae8aa0c2ccc4e486edfef0aed0b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
 
 ########################### Executables
 
-add_executable(numerics DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp  Vector/Vector_unit_tests.cpp  Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp  Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu  ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp DCPSE_op/DCPSE_Solver.hpp MatrixAssembler/MatrixAssembler.hpp)
+add_executable(numerics DCPSE_op/DCPSE_op_test3d.cpp DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp  Vector/Vector_unit_tests.cpp  Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp  Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu  ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp DCPSE_op/DCPSE_Solver.hpp MatrixAssembler/MatrixAssembler.hpp)
 
 
 ###########################
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 2acca4bd156428f647f2e15e3d5df76bddfe34eb..810331a0a806ad2d94dafbbee8d00bbf17bf57e9 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -108,7 +108,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
         openfpm::vector<aggregate<int>> ref_p;
 
         constexpr int x      =     0;
-        constexpr int y          =     1;
+        constexpr int y      =     1;
 
         constexpr int Polarization      =     0;
         constexpr int Velocity          =     1;
@@ -245,15 +245,15 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
 
             dV[x] = 0.5*Dy(h[y]) + zeta*Dx(delmu*Pol[x]*Pol[x]) + zeta*Dy(delmu*Pol[x]*Pol[y]) - zeta*Dx(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dx(-2*h[y]*Pol[x]*Pol[y])
                   - 0.5*nu*Dy(h[y]*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - Dx(sigma[x][x]) - Dy(sigma[x][y]) - g[x]
-                  - 0.5*nu*Dx(-gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dy(-2*gama*lambda*delmu*(Pol[x]*Pol[y])) ;
+                  - 0.5*nu*Dx(-gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dy(-2*gama*lambda*delmu*(Pol[x]*Pol[y]));
 
 
             dV[y] = 0.5*Dx(-h[y]) + zeta*Dy(delmu*Pol[y]*Pol[y]) + zeta*Dx(delmu*Pol[x]*Pol[y]) - zeta*Dy(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dy(-2*h[y]*Pol[x]*Pol[y])
                   - 0.5*nu*Dx(h[y]*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - Dx(sigma[y][x]) - Dy(sigma[y][y]) - g[y]
-                  - 0.5*nu*Dy(gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dx(-2*gama*lambda*delmu*(Pol[x]*Pol[y])) ;
+                  - 0.5*nu*Dy(gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dx(-2*gama*lambda*delmu*(Pol[x]*Pol[y]));
 
 
-/*        //Velocity Solution n iterations
+/*      //Velocity Solution n iterations
         eq_id vx,vy;
         vx.setId(0);
         vy.setId(1);
@@ -1390,7 +1390,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
             ++it2;
         }
 
-
+        domain.ghost_get<1>();
         auto Poisson = -Lap(v);
         Solver.impose(Poisson, bulk, prop_id<1>());
         Solver.impose(v, up_p, 0);
@@ -1408,6 +1408,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
 
         }
         std::cout << "Maximum Auto Error: " << worst1 << std::endl;
+        domain.ghost_get<0>();
 
         domain.write("Poisson_Periodic");
     }
diff --git a/src/DCPSE_op/DCPSE_op_test3d.cpp b/src/DCPSE_op/DCPSE_op_test3d.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8cbfd539b89a39d714b5fc4bed17251e87011213
--- /dev/null
+++ b/src/DCPSE_op/DCPSE_op_test3d.cpp
@@ -0,0 +1,115 @@
+/*
+ * DCPSE_op_test.cpp
+ *
+ *  Created on: April 9, 2020
+ *      Author: Abhinav Singh
+ *
+ */
+
+#include "config.h"
+
+#define BOOST_TEST_DYN_LINK
+
+#include "util/util_debug.hpp"
+#include <boost/test/unit_test.hpp>
+#include <iostream>
+#include "DCPSE_op.hpp"
+#include "DCPSE_Solver.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+#include "Vector/vector_dist_subset.hpp"
+
+//int vector_dist_expression_op<void,void,VECT_COPY_N_TO_N>::i = 0;
+//int vector_dist_expression_op<void,void,VECT_COPY_1_TO_N>::i = 0;
+
+//! Specify the general characteristic of system to solve
+struct equations3d2 {
+    //! dimensionaly of the equation ( 3D problem ...)
+    static const unsigned int dims = 3;
+    //! number of fields in the system
+    static const unsigned int nvar = 3;
+
+    //! boundary at X and Y
+    static const bool boundary[];
+
+    //! type of space float, double, ...
+    typedef double stype;
+
+    //! type of base particles
+    typedef vector_dist<dims, double, aggregate<double>> b_part;
+
+    //! type of SparseMatrix for the linear solver
+    typedef SparseMatrix<double, int, EIGEN_BASE> SparseMatrix_type;
+
+    //! type of Vector for the linear solver
+    typedef Vector<double> Vector_type;
+};
+
+struct equations3d {
+    //! dimensionaly of the equation ( 3D problem ...)
+    static const unsigned int dims = 3;
+    //! number of fields in the system
+    static const unsigned int nvar = 1;
+
+    //! boundary at X and Y
+    static const bool boundary[];
+
+    //! type of space float, double, ...
+    typedef double stype;
+
+    //! type of base particles
+    typedef vector_dist<dims, double, aggregate<double>> b_part;
+
+    //! type of SparseMatrix for the linear solver
+    typedef SparseMatrix<double, int, EIGEN_BASE> SparseMatrix_type;
+
+    //! type of Vector for the linear solver
+    typedef Vector<double> Vector_type;
+};
+
+const bool equations3d::boundary[] = {NON_PERIODIC, NON_PERIODIC};
+const bool equations3d2::boundary[] = {NON_PERIODIC, NON_PERIODIC};
+
+
+//template<typename T>
+//struct Debug;
+
+BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
+
+    BOOST_AUTO_TEST_CASE(dcpse_sphere) {
+        const size_t sz[3] = {31,31,31};
+        Sphere<3, double> sphere(0, 5);
+        double spacing = 0.1;
+        Ghost<3, double> ghost(spacing * 3);
+        double rCut = 2.5 * spacing;
+        //                                  P        V                 v_star           RHS            V_t       Helmholtz
+        vector_dist<2, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>, double,    double>> Particles(0, Sphere, bc, ghost);
+        auto it = Particles.getGridIterator(sz);
+        while (it.isNext()) {
+            Particles.add();
+            auto key = it.get();
+            double x = key.get(0) * it.getSpacing(0);
+            Particles.getLastPos()[0] = r*cos(theta)*cos(phi);
+            double y = key.get(1) * it.getSpacing(1);
+            Particles.getLastPos()[1] = r*sin(theta)*cos(phi);
+            double z = key.get(1) * it.getSpacing(1);
+            Particles.getLastPos()[1] = r*sin(phi);
+            ++it;
+        }
+
+        Sphere<3, double> bulk(0,5);
+
+        openfpm::vector<Sphere<3, double>> spheres;
+        spheres.add(bulk);
+        VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_sph;
+        vtk_sph.add(spheres);
+        vtk_sph.write("vtk_sph.vtk");
+
+        Particles.map();
+        Particles.ghost_get<0>();
+
+        Particles.write_frame("Sphere",i);\
+    }
+
+BOOST_AUTO_TEST_SUITE_END()
+
+