diff --git a/src/Integration/AdamsBashforth2/AB2Kernels.hpp b/src/Integration/AdamsBashforth2/AB2Kernels.hpp
new file mode 100644
index 00000000..578cf07e
--- /dev/null
+++ b/src/Integration/AdamsBashforth2/AB2Kernels.hpp
@@ -0,0 +1,63 @@
+
+
+#ifndef __AB2Kernels_hpp__
+#define __AB2Kernels_hpp__
+
+#include "KokkosTypes.hpp"
+#include "types.hpp"
+
+namespace pFlow::AB2Kernels
+{
+inline
+bool intAllActive(
+    const word& name,
+	real dt, 
+    rangeU32 activeRng,
+	const deviceViewType1D<realx3>& y, 
+	const deviceViewType1D<realx3>& dy,
+	const deviceViewType1D<realx3>& dy1
+)
+{
+	Kokkos::parallel_for(
+		name,
+		deviceRPolicyStatic (activeRng.start(), activeRng.end()),
+		LAMBDA_HD(uint32 i){
+			y[i] +=  dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
+			dy1[i] = dy[i];
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+inline
+bool intScattered
+(
+	const word& name,
+	real dt, 
+    const pFlagTypeDevice& activePoints,
+	const deviceViewType1D<realx3>& y, 
+	const deviceViewType1D<realx3>& dy,
+	const deviceViewType1D<realx3>& dy1
+)
+{
+
+	Kokkos::parallel_for(
+		name,
+		deviceRPolicyStatic (activePoints.activeRange().start(), activePoints.activeRange().end()),
+		LAMBDA_HD(uint32 i){
+			if( activePoints(i))
+			{
+				y[i] +=  dt*(static_cast<real>(1.5) * dy[i] - static_cast<real>(0.5) * dy1[i]);
+				dy1[i] = dy[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+}
+
+#endif
\ No newline at end of file
diff --git a/src/Integration/boundaries/boundaryIntegration.cpp b/src/Integration/boundaries/boundaryIntegration.cpp
new file mode 100644
index 00000000..e13a5624
--- /dev/null
+++ b/src/Integration/boundaries/boundaryIntegration.cpp
@@ -0,0 +1,55 @@
+#include "boundaryIntegration.hpp"
+#include "pointStructure.hpp"
+
+
+pFlow::boundaryIntegration::boundaryIntegration(
+    const boundaryBase &boundary, 
+    const pointStructure &pStruct,  
+    const word &method,
+    integration& intgrtn
+)
+:
+    generalBoundary(boundary, pStruct, "", method),
+    integration_(intgrtn)
+{}
+
+pFlow::uniquePtr<pFlow::boundaryIntegration> pFlow::boundaryIntegration::create(
+    const boundaryBase &boundary, 
+    const pointStructure &pStruct,  
+    const word &method,
+    integration& intgrtn
+)
+{
+
+    word bType = angleBracketsNames2( 
+		"boundaryIntegration",
+		boundary.type(),
+        method);
+
+	word altBType{"boundaryIntegration<none>"};
+
+	if( boundaryBasevCtorSelector_.search(bType) )
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(bType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[bType](boundary, pStruct, method, intgrtn);
+	}
+	else if(boundaryBasevCtorSelector_.search(altBType))
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(altBType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[altBType](boundary, pStruct, method, intgrtn);
+	}
+	else
+	{
+		printKeys( 
+			fatalError << "Ctor Selector "<< bType<<
+			" and "<< altBType << " do not exist. \n"
+			<<"Avaiable ones are: \n",
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+	}
+	
+	return nullptr;
+}
diff --git a/src/Integration/boundaries/boundaryIntegration.hpp b/src/Integration/boundaries/boundaryIntegration.hpp
new file mode 100644
index 00000000..b3a20805
--- /dev/null
+++ b/src/Integration/boundaries/boundaryIntegration.hpp
@@ -0,0 +1,91 @@
+
+#ifndef __boundaryIntegration_hpp__
+#define __boundaryIntegration_hpp__
+
+
+#include "generalBoundary.hpp"
+#include "virtualConstructor.hpp"
+#include "pointFields.hpp"
+
+namespace pFlow
+{
+
+class integration;
+
+class boundaryIntegration
+:
+    public generalBoundary
+{
+private:
+
+	const integration& 		integration_;
+
+public:
+
+    TypeInfo("boundaryIntegration<none>");
+
+    boundaryIntegration(
+        const boundaryBase& boundary,
+        const pointStructure& pStruct, 
+        const word& method,
+        integration& intgrtn);
+
+    ~boundaryIntegration()override = default;
+
+    create_vCtor(
+        boundaryIntegration,
+        boundaryBase,
+        (
+            const boundaryBase& boundary,
+            const pointStructure& pStruct, 
+            const word& method,
+			integration& intgrtn
+        ),
+        (boundary, pStruct, method, intgrtn)
+    );
+    
+    add_vCtor(
+        boundaryIntegration,
+        boundaryIntegration,
+        boundaryBase
+    );
+
+    bool hearChanges(
+        real t,
+        real dt,
+        uint32 iter,
+        const message &msg,
+        const anyList &varList) override
+    {
+        return true;
+    }
+
+	const integration& Integration()const
+	{
+		return integration_;
+	}
+
+	virtual
+	bool correct(real dt, const realx3PointField_D& y, const realx3PointField_D& dy)
+	{
+		return true;
+	}
+
+    virtual
+    bool correctPStruct(real dt, const realx3PointField_D& vel)
+    {
+        return true;
+    }
+
+    static
+    uniquePtr<boundaryIntegration> create(
+        const boundaryBase& boundary,
+        const pointStructure& pStruct, 
+        const word& method,
+		integration& intgrtn);
+
+};
+
+}
+
+#endif
\ No newline at end of file
diff --git a/src/Integration/boundaries/boundaryIntegrationList.cpp b/src/Integration/boundaries/boundaryIntegrationList.cpp
new file mode 100644
index 00000000..c7e2ad3b
--- /dev/null
+++ b/src/Integration/boundaries/boundaryIntegrationList.cpp
@@ -0,0 +1,41 @@
+#include "boundaryIntegrationList.hpp"
+#include "integration.hpp"
+
+pFlow::boundaryIntegrationList::boundaryIntegrationList(
+    const pointStructure &pStruct, 
+    const word &method, 
+    integration &intgrtn
+)
+:
+    ListPtr<boundaryIntegration>(6),
+	boundaries_(pStruct.boundaries())
+{
+    
+    for(uint32 i=0; i<6; i++)
+	{
+		this->set(
+			i,
+			boundaryIntegration::create(
+				boundaries_[i],
+				pStruct, 
+				method,
+                intgrtn
+            )
+        );
+	}	
+
+}
+
+bool pFlow::boundaryIntegrationList::correct(real dt, realx3PointField_D &y, realx3PointField_D &dy)
+{
+	for(auto& bndry:*this)
+	{
+		if(!bndry->correct(dt, y, dy))
+		{
+			fatalErrorInFunction;
+			return false;
+		}
+		
+    }
+    return true;
+}
diff --git a/src/Integration/boundaries/boundaryIntegrationList.hpp b/src/Integration/boundaries/boundaryIntegrationList.hpp
new file mode 100644
index 00000000..a4e262a0
--- /dev/null
+++ b/src/Integration/boundaries/boundaryIntegrationList.hpp
@@ -0,0 +1,45 @@
+
+#ifndef __boundaryIntegrationList_hpp__
+#define __boundaryIntegrationList_hpp__
+
+
+#include "boundaryList.hpp"
+#include "ListPtr.hpp"
+#include "boundaryIntegration.hpp"
+
+
+namespace pFlow
+{
+
+class integration;
+
+class boundaryIntegrationList
+:
+    public ListPtr<boundaryIntegration>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundaryIntegrationList(
+        const pointStructure& pStruct, 
+        const word& method,
+		integration& intgrtn
+    );
+
+    ~boundaryIntegrationList()=default;
+
+    bool correct(
+		real dt, 
+		realx3PointField_D& y, 
+		realx3PointField_D& dy);
+};
+
+
+
+}
+
+
+#endif //__boundaryIntegrationList_hpp__
\ No newline at end of file
diff --git a/src/Interaction/CMakeLists.txt b/src/Interaction/CMakeLists.txt
index c0265808..3b985647 100644
--- a/src/Interaction/CMakeLists.txt
+++ b/src/Interaction/CMakeLists.txt
@@ -7,8 +7,7 @@ contactSearch/methods/cellBased/NBS/NBS.cpp
 contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp
 
 contactSearch/boundaries/boundaryContactSearch/boundaryContactSearch.cpp
-contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
-contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
+
 contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearchKernels.cpp
 contactSearch/boundaries/periodicBoundaryContactSearch/ppwBndryContactSearch.cpp
 contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp
@@ -28,6 +27,8 @@ if(pFlow_Build_MPI)
     list(APPEND SourceFiles 
     contactSearch/boundaries/processorBoundaryContactSearch/processorBoundaryContactSearch.cpp
     sphereInteraction/boundaries/processorBoundarySphereInteraction/processorBoundarySphereInteractions.cpp
+    contactSearch/boundaries/twoPartContactSearch/twoPartContactSearchKernels.cpp
+    contactSearch/boundaries/twoPartContactSearch/twoPartContactSearch.cpp
     )
 endif()
 
diff --git a/src/Particles/SphereParticles/boundarySphereParticles.cpp b/src/Particles/SphereParticles/boundarySphereParticles.cpp
new file mode 100644
index 00000000..8eb5eaa0
--- /dev/null
+++ b/src/Particles/SphereParticles/boundarySphereParticles.cpp
@@ -0,0 +1,64 @@
+#include "boundarySphereParticles.hpp"
+#include "boundaryBase.hpp"
+#include "sphereParticles.hpp"
+
+
+pFlow::boundarySphereParticles::boundarySphereParticles(
+    const boundaryBase &boundary, 
+    sphereParticles &prtcls
+)
+:
+    generalBoundary(boundary, prtcls.pStruct(), "", ""),
+    particles_(prtcls)
+{
+
+}
+
+pFlow::sphereParticles &pFlow::boundarySphereParticles::Particles()
+{
+    return particles_;
+}
+
+const pFlow::sphereParticles &pFlow::boundarySphereParticles::Particles() const
+{
+    return particles_;
+}
+
+pFlow::uniquePtr<pFlow::boundarySphereParticles> pFlow::boundarySphereParticles::create(
+    const boundaryBase &boundary, 
+    sphereParticles &prtcls
+)
+{
+	
+	word bType = angleBracketsNames2( 
+		"boundarySphereParticles",
+		pFlowProcessors().localRunTypeName(), 
+		boundary.type());
+
+	word altBType{"boundarySphereParticles<none>"};
+
+	if( boundaryBasevCtorSelector_.search(bType) )
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(bType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[bType](boundary, prtcls);
+	}
+	else if(boundaryBasevCtorSelector_.search(altBType))
+	{
+		pOutput.space(4)<<"Creating boundary "<< Green_Text(altBType)<<
+		" for "<<boundary.name()<<endl;
+		return boundaryBasevCtorSelector_[altBType](boundary, prtcls);
+	}
+	else
+	{
+		printKeys( 
+			fatalError << "Ctor Selector "<< bType<<
+			" and "<< altBType << " do not exist. \n"
+			<<"Avaiable ones are: \n",
+			boundaryBasevCtorSelector_
+		);
+		fatalExit;
+	}
+	
+	return nullptr;
+}
diff --git a/src/Particles/SphereParticles/boundarySphereParticles.hpp b/src/Particles/SphereParticles/boundarySphereParticles.hpp
new file mode 100644
index 00000000..f4662307
--- /dev/null
+++ b/src/Particles/SphereParticles/boundarySphereParticles.hpp
@@ -0,0 +1,80 @@
+
+
+#ifndef __boundarySphereParticles_hpp__
+#define __boundarySphereParticles_hpp__
+
+#include "generalBoundary.hpp"
+#include "virtualConstructor.hpp"
+#include "timeInfo.hpp"
+
+namespace pFlow
+{
+
+class sphereParticles;
+
+class boundarySphereParticles
+    : public generalBoundary
+{
+private:
+
+    sphereParticles&        particles_;
+
+public:
+
+    /// type info
+    TypeInfo("boundarySphereParticles<none>");
+
+    boundarySphereParticles(
+        const boundaryBase &boundary,
+        sphereParticles&  prtcls
+    );
+
+    create_vCtor(
+        boundarySphereParticles,
+        boundaryBase,
+        (
+            const boundaryBase &boundary,
+            sphereParticles&  prtcls
+        ),
+        (boundary, prtcls)
+    );
+
+    add_vCtor(
+        boundarySphereParticles,
+        boundarySphereParticles,
+        boundaryBase
+    );
+
+    sphereParticles& Particles();
+
+    const sphereParticles& Particles()const;
+
+    bool hearChanges(
+        real t,
+        real dt,
+        uint32 iter,
+        const message &msg,
+        const anyList &varList) override
+    {
+        return true;
+    }
+
+    virtual
+    bool acceleration(const timeInfo& ti, const realx3& g)
+    {
+        return true;
+    }
+
+    static
+    uniquePtr<boundarySphereParticles> create(
+        const boundaryBase &boundary,
+        sphereParticles&  prtcls);
+    
+};
+
+
+}
+
+
+
+#endif
diff --git a/src/Particles/SphereParticles/boundarySphereParticlesList.cpp b/src/Particles/SphereParticles/boundarySphereParticlesList.cpp
new file mode 100644
index 00000000..5cbebc32
--- /dev/null
+++ b/src/Particles/SphereParticles/boundarySphereParticlesList.cpp
@@ -0,0 +1,19 @@
+#include "boundarySphereParticlesList.hpp"
+
+pFlow::boundarySphereParticlesList::boundarySphereParticlesList(
+    const boundaryList &bndrs, 
+    sphereParticles &prtcls
+)
+:
+    ListPtr(bndrs.size()), 
+    boundaries_(bndrs)
+{
+    for(auto i=0; i<boundaries_.size(); i++)
+    {
+        this->set
+        (
+            i, 
+            boundarySphereParticles::create(boundaries_[i], prtcls)
+        );
+    }
+}
\ No newline at end of file
diff --git a/src/Particles/SphereParticles/boundarySphereParticlesList.hpp b/src/Particles/SphereParticles/boundarySphereParticlesList.hpp
new file mode 100644
index 00000000..cd6770df
--- /dev/null
+++ b/src/Particles/SphereParticles/boundarySphereParticlesList.hpp
@@ -0,0 +1,36 @@
+
+
+#ifndef __boundarySphereParticlesList_hpp__
+#define __boundarySphereParticlesList_hpp__
+
+#include "ListPtr.hpp"
+#include "boundaryList.hpp"
+#include "boundarySphereParticles.hpp"
+
+namespace pFlow
+{
+
+class boundarySphereParticlesList
+:
+    public ListPtr<boundarySphereParticles>
+{
+private:
+
+    const boundaryList&         boundaries_;
+
+public:
+
+    boundarySphereParticlesList(
+        const boundaryList& bndrs,
+        sphereParticles& prtcls
+    );
+
+    ~boundarySphereParticlesList()=default;
+
+};
+
+}
+
+
+
+#endif 
\ No newline at end of file
diff --git a/src/phasicFlow/globals/boundaryConfigs.hpp b/src/phasicFlow/globals/boundaryConfigs.hpp
new file mode 100644
index 00000000..11574957
--- /dev/null
+++ b/src/phasicFlow/globals/boundaryConfigs.hpp
@@ -0,0 +1,10 @@
+#ifndef __boundaryConfigs_hpp__
+#define __boundaryConfigs_hpp__
+
+
+#ifndef BoundaryModel1
+//#define BoundaryModel1
+#endif
+
+
+#endif //__boundaryConfigs_hpp__
\ No newline at end of file
diff --git a/src/phasicFlow/processors/MPITimer.cpp b/src/phasicFlow/processors/MPITimer.cpp
new file mode 100644
index 00000000..cfa55631
--- /dev/null
+++ b/src/phasicFlow/processors/MPITimer.cpp
@@ -0,0 +1,88 @@
+#include "MPITimer.hpp"
+
+#ifdef pFlow_Build_MPI
+    #include "pFlowProcessors.hpp"
+    #include "procCommunication.hpp"
+#endif
+
+pFlow::real pFlow::MPITimer::totalTimeMax() const
+{
+    return accTimersTotalMax();
+}
+
+std::vector<pFlow::real> pFlow::MPITimer::totalTimeAllToAll() const
+{
+#ifdef pFlow_Build_MPI
+    MPI::procCommunication comm(pFlowProcessors());
+    auto [totTime, succs] = comm.collectAllToAll(totalTime());
+    return totTime;
+#else
+    return std::vector<real>(1, totalTime());
+#endif
+}
+
+std::vector<pFlow::real> pFlow::MPITimer::totalTimeAllToMaster() const
+{
+#ifdef pFlow_Build_MPI
+	MPI::procCommunication comm(pFlowProcessors());
+	auto [totTime, succs] = comm.collectAllToMaster(totalTime());
+	return totTime;
+#else
+		return std::vector<real>(1, totalTime());
+#endif
+}
+
+pFlow::real pFlow::MPITimer::averageTimeMax() const
+{
+    return Timer::averageTimeMax();
+}
+
+std::vector<pFlow::real> pFlow::MPITimer::averageTimeAllToAll() const
+{
+#ifdef pFlow_Build_MPI
+    MPI::procCommunication comm(pFlowProcessors());
+    auto [totTime, succs] = comm.collectAllToAll(averageTime());
+    return totTime;
+#else
+    return std::vector<real>(1, averageTime());
+#endif
+}
+
+std::vector<pFlow::real> pFlow::MPITimer::averageTimeAllAtoMaster() const
+{
+#ifdef pFlow_Build_MPI
+	MPI::procCommunication comm(pFlowProcessors());
+	auto [totTime, succs] = comm.collectAllToMaster(averageTime());
+	return totTime;
+#else
+		return std::vector<real>(1, averageTime());
+#endif
+}
+
+bool pFlow::MPITimer::write(iOstream &os) const
+{
+	const auto ts = totalTimeAllToAll();
+	auto maxloc = std::distance(ts.begin(), std::max_element(ts.begin(), ts.end()));
+	os<<'(';
+	for(auto i=0; i<ts.size(); i++)
+	{
+		if(maxloc == i)
+			os<<Red_Text(ts[i]);
+		else
+			os<<ts[i];
+		
+		if(i != ts.size()-1)
+			os<<' ';
+		else
+			os<<')'<<endl;
+	}
+
+    return true;
+}
+
+static pFlow::MPITimer ComputationTimer__{"ComputationTimer"};
+
+pFlow::MPITimer &pFlow::ComputationTimer()
+{
+    return ComputationTimer__;
+}
diff --git a/src/phasicFlow/processors/MPITimer.hpp b/src/phasicFlow/processors/MPITimer.hpp
new file mode 100644
index 00000000..a3d1f306
--- /dev/null
+++ b/src/phasicFlow/processors/MPITimer.hpp
@@ -0,0 +1,55 @@
+
+
+#ifndef __MPITimer_hpp__
+#define __MPITimer_hpp__
+
+#include "Timer.hpp"
+
+
+namespace pFlow
+{
+
+class MPITimer
+:
+    public Timer
+{
+private:
+
+    // hiding methods
+    using Timer::accTimersTotal;
+    using Timer::averageTime;
+
+public:
+
+    TypeInfo("MPITimer");
+
+    explicit MPITimer(const word& name)
+	: 
+    Timer(name)
+	{}
+
+	
+    real totalTimeMax()const;
+
+    std::vector<real> totalTimeAllToAll()const;
+
+    std::vector<real> totalTimeAllToMaster()const;
+
+    real averageTimeMax()const;
+
+    std::vector<real> averageTimeAllToAll()const;
+
+    std::vector<real> averageTimeAllAtoMaster()const;
+
+    // call this from all processors in pFlowProcessors
+    bool write(iOstream& os)const;
+
+};
+
+MPITimer& ComputationTimer();
+
+}
+
+
+
+#endif
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/boundaries/boundariesMask.hpp b/src/phasicFlow/structuredData/boundaries/boundariesMask.hpp
new file mode 100644
index 00000000..13875f2b
--- /dev/null
+++ b/src/phasicFlow/structuredData/boundaries/boundariesMask.hpp
@@ -0,0 +1,50 @@
+
+#ifndef __boundariesMask_hpp__
+#define __boundariesMask_hpp__
+
+#include <array>
+
+namespace pFlow
+{
+
+template<size_t N>
+class boundariesMask
+:
+    public std::array<bool,N> 
+{
+public:
+    
+    boundariesMask()=default;
+
+    boundariesMask(bool val)
+    {
+        this->fill(val);
+    }
+
+    boundariesMask(std::initializer_list<bool> l)
+    :
+        std::array<bool,N>(l)
+    {}
+
+    inline
+    bool allElements(bool val)
+    {
+        return std::all_of(this->begin(), this->end(), [val](bool v) { return v==val;} );
+    }
+
+    inline
+    bool anyElement(bool val)
+    {
+        return std::any_of(this->begin(), this->end(), [val](bool v) { return v==val;} );
+    }
+    
+    inline
+    bool noElement(bool val)
+    {
+        return std::none_of(this->begin(), this->end(), [val](bool v) { return v==val;} );
+    }
+};
+
+}
+
+#endif //__boundariesMask_hpp__
\ No newline at end of file