diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.C b/src/Integration/AdamsBashforth2/AdamsBashforth2.C
index baf76a54..7739e219 100644
--- a/src/Integration/AdamsBashforth2/AdamsBashforth2.C
+++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.C
@@ -74,6 +74,13 @@ bool pFlow::AdamsBashforth2::correct
 	return true;
 }
 
+bool pFlow::AdamsBashforth2::setInitialVals(
+	const int32IndexContainer& newIndices,
+	const realx3Vector& y)
+{
+	return true;
+}
+
 bool pFlow::AdamsBashforth2::intAll(
 	real dt, 
 	realx3Vector_D& y, 
diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.H b/src/Integration/AdamsBashforth2/AdamsBashforth2.H
index 22dd792c..aa191f4c 100644
--- a/src/Integration/AdamsBashforth2/AdamsBashforth2.H
+++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.H
@@ -68,6 +68,20 @@ public:
 
 		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
 
+		bool setInitialVals(
+			const int32IndexContainer& newIndices,
+			const realx3Vector& y) override;
+
+		bool needSetInitialVals()const override
+		{
+			return false;
+		}
+
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth2>(*this);
+		}
+
 		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
 
 		template<typename activeFunctor>
diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.C b/src/Integration/AdamsBashforth3/AdamsBashforth3.C
index c772ad13..fa8e9927 100644
--- a/src/Integration/AdamsBashforth3/AdamsBashforth3.C
+++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.C
@@ -60,7 +60,7 @@ pFlow::AdamsBashforth3::AdamsBashforth3
 			AB3History({zero3,zero3})))
 
 {
-
+	
 }
 
 bool pFlow::AdamsBashforth3::predict
@@ -94,6 +94,13 @@ bool pFlow::AdamsBashforth3::correct
 	return true;
 }
 
+bool pFlow::AdamsBashforth3::setInitialVals(
+	const int32IndexContainer& newIndices,
+	const realx3Vector& y)
+{
+	return true;
+}
+
 bool pFlow::AdamsBashforth3::intAll(
 	real dt, 
 	realx3Vector_D& y, 
diff --git a/src/Integration/AdamsBashforth3/AdamsBashforth3.H b/src/Integration/AdamsBashforth3/AdamsBashforth3.H
index 6b9d9e94..d9baecb7 100644
--- a/src/Integration/AdamsBashforth3/AdamsBashforth3.H
+++ b/src/Integration/AdamsBashforth3/AdamsBashforth3.H
@@ -117,6 +117,20 @@ public:
 			realx3Vector_D & y,
 			realx3Vector_D& dy) override;
 
+		bool setInitialVals(
+			const int32IndexContainer& newIndices,
+			const realx3Vector& y) override;
+
+		bool needSetInitialVals()const override
+		{
+			return false;
+		}
+
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth3>(*this);
+		}
+
 		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.C b/src/Integration/AdamsBashforth4/AdamsBashforth4.C
index 13d25435..313eb58c 100644
--- a/src/Integration/AdamsBashforth4/AdamsBashforth4.C
+++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.C
@@ -76,6 +76,13 @@ bool pFlow::AdamsBashforth4::correct
 	return true;
 }
 
+bool pFlow::AdamsBashforth4::setInitialVals(
+	const int32IndexContainer& newIndices,
+	const realx3Vector& y)
+{
+	return true;
+}
+
 bool pFlow::AdamsBashforth4::intAll(
 	real dt, 
 	realx3Vector_D& y, 
diff --git a/src/Integration/AdamsBashforth4/AdamsBashforth4.H b/src/Integration/AdamsBashforth4/AdamsBashforth4.H
index 00b7c752..7cb31e88 100644
--- a/src/Integration/AdamsBashforth4/AdamsBashforth4.H
+++ b/src/Integration/AdamsBashforth4/AdamsBashforth4.H
@@ -122,6 +122,20 @@ public:
 			realx3Vector_D & y,
 			realx3Vector_D& dy) override;
 
+		bool setInitialVals(
+			const int32IndexContainer& newIndices,
+			const realx3Vector& y) override;
+
+		bool needSetInitialVals()const override
+		{
+			return false;
+		}
+
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth4>(*this);
+		}
+
 		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
diff --git a/src/Integration/AdamsBashforth5/AdamsBashforth5.C b/src/Integration/AdamsBashforth5/AdamsBashforth5.C
index 19971468..23f7fc09 100644
--- a/src/Integration/AdamsBashforth5/AdamsBashforth5.C
+++ b/src/Integration/AdamsBashforth5/AdamsBashforth5.C
@@ -76,6 +76,13 @@ bool pFlow::AdamsBashforth5::correct
 	return true;
 }
 
+bool pFlow::AdamsBashforth5::setInitialVals(
+	const int32IndexContainer& newIndices,
+	const realx3Vector& y)
+{
+	return true;
+}
+
 bool pFlow::AdamsBashforth5::intAll(
 	real dt, 
 	realx3Vector_D& y, 
diff --git a/src/Integration/AdamsBashforth5/AdamsBashforth5.H b/src/Integration/AdamsBashforth5/AdamsBashforth5.H
index c389e8c4..aee15244 100644
--- a/src/Integration/AdamsBashforth5/AdamsBashforth5.H
+++ b/src/Integration/AdamsBashforth5/AdamsBashforth5.H
@@ -120,6 +120,20 @@ public:
 			realx3Vector_D & y,
 			realx3Vector_D& dy) override;
 
+		bool setInitialVals(
+			const int32IndexContainer& newIndices,
+			const realx3Vector& y) override;
+
+		bool needSetInitialVals()const override
+		{
+			return false;
+		}
+
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsBashforth5>(*this);
+		}
+
 		bool intAll(real dt, 
 			realx3Vector_D& y, 
 			realx3Vector_D& dy, 
diff --git a/src/Integration/AdamsMoulton3/AdamsMoulton3.C b/src/Integration/AdamsMoulton3/AdamsMoulton3.C
new file mode 100644
index 00000000..aa6d45d9
--- /dev/null
+++ b/src/Integration/AdamsMoulton3/AdamsMoulton3.C
@@ -0,0 +1,175 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#include "AdamsMoulton3.H"
+
+//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0};
+
+pFlow::AdamsMoulton3::AdamsMoulton3
+(
+	const word& baseName,
+	repository& owner,
+	const pointStructure& pStruct,
+	const word& method
+)
+:
+	integration(baseName, owner, pStruct, method),
+	y0_(
+		owner.emplaceObject<realx3PointField_D>(
+			objectFile(
+				groupNames(baseName,"y0"),
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_ALWAYS),
+			pStruct,
+			zero3,
+			false
+			)
+		),
+	dy0_(
+		owner.emplaceObject<realx3PointField_D>(
+			objectFile(
+				groupNames(baseName,"dy0"),
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_ALWAYS),
+			pStruct,
+			zero3
+			)
+		),
+	dy1_(
+		owner.emplaceObject<realx3PointField_D>(
+			objectFile(
+				groupNames(baseName,"dy1"),
+				"",
+				objectFile::READ_IF_PRESENT,
+				objectFile::WRITE_ALWAYS),
+			pStruct,
+			zero3
+			)
+		)
+{
+
+}
+
+bool pFlow::AdamsMoulton3::predict
+(
+	real dt,
+	realx3Vector_D& y,
+	realx3Vector_D& dy
+)
+{
+
+	if(this->pStruct().allActive())
+	{
+		return predictAll(dt, y, dy, this->pStruct().activeRange());
+	}
+	else
+	{
+		return predictRange(dt, y, dy, this->pStruct().activePointsMaskD());
+	}
+
+	return true;
+}
+
+bool pFlow::AdamsMoulton3::correct
+(
+	real dt,
+	realx3Vector_D& y,
+	realx3Vector_D& dy
+)
+{
+	if(this->pStruct().allActive())
+	{
+		return intAll(dt, y, dy, this->pStruct().activeRange());
+	}
+	else
+	{
+		return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
+	}
+
+	return true;
+}
+
+bool pFlow::AdamsMoulton3::setInitialVals(
+	const int32IndexContainer& newIndices,
+	const realx3Vector& y)
+{
+	y0_.insertSetElement(newIndices, y);
+	//output<< "y0_  "<< y0_<<endl<<endl;
+	return true;
+}
+
+bool pFlow::AdamsMoulton3::predictAll(
+	real dt, 
+	realx3Vector_D& y, 
+	realx3Vector_D& dy, 
+	range activeRng)
+{
+
+	auto d_dy = dy.deviceVectorAll();
+	auto d_y  = y.deviceVectorAll();
+	auto d_y0 = y0_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceVectorAll();
+	auto d_dy1= dy1_.deviceVectorAll();
+
+	Kokkos::parallel_for(
+		"AdamsMoulton3::predict",
+		rpIntegration (activeRng.first, activeRng.second),
+		LAMBDA_HD(int32 i){
+			d_dy0[i] = d_dy[i];
+			d_y[i] = d_y0[i] + dt*(static_cast<real>(3.0 / 2.0) * d_dy[i] - static_cast<real>(1.0 / 2.0) * d_dy1[i]);
+			
+		});
+	Kokkos::fence();
+
+	return true;	
+}
+
+bool pFlow::AdamsMoulton3::intAll(
+	real dt, 
+	realx3Vector_D& y, 
+	realx3Vector_D& dy, 
+	range activeRng)
+{
+
+	auto d_dy = dy.deviceVectorAll();
+	auto d_y  = y.deviceVectorAll();
+
+	auto d_dy0 = dy0_.deviceVectorAll();
+	auto d_y0  = y0_.deviceVectorAll();
+	auto d_dy1 = dy1_.deviceVectorAll();
+
+	Kokkos::parallel_for(
+		"AdamsMoulton3::correct",
+		rpIntegration (activeRng.first, activeRng.second),
+		LAMBDA_HD(int32 i){
+			auto corrct_y = d_y0[i] + dt*(
+				  static_cast<real>(5.0/12.0)*d_dy[i] 
+				+ static_cast<real>(8.0/12.0)*d_dy0[i] 
+				- static_cast<real>(1.0/12.0)*d_dy1[i]);
+			d_y[i]  = corrct_y;
+			d_y0[i] = corrct_y;
+			d_dy1[i]= d_dy0[i];
+		});
+	Kokkos::fence();
+	
+	return true;	
+}
\ No newline at end of file
diff --git a/src/Integration/AdamsMoulton3/AdamsMoulton3.H b/src/Integration/AdamsMoulton3/AdamsMoulton3.H
new file mode 100644
index 00000000..c6891966
--- /dev/null
+++ b/src/Integration/AdamsMoulton3/AdamsMoulton3.H
@@ -0,0 +1,179 @@
+/*------------------------------- phasicFlow ---------------------------------
+      O        C enter of
+     O O       E ngineering and
+    O   O      M ultiscale modeling of
+   OOOOOOO     F luid flow       
+------------------------------------------------------------------------------
+  Copyright (C): www.cemf.ir
+  email: hamid.r.norouzi AT gmail.com
+------------------------------------------------------------------------------  
+Licence:
+  This file is part of phasicFlow code. It is a free software for simulating 
+  granular and multiphase flows. You can redistribute it and/or modify it under
+  the terms of GNU General Public License v3 or any other later versions. 
+ 
+  phasicFlow is distributed to help others in their research in the field of 
+  granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
+  implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+
+-----------------------------------------------------------------------------*/
+
+#ifndef __AdamsMoulton3_H__
+#define __AdamsMoulton3_H__
+
+
+#include "integration.H"
+#include "pointFields.H"
+
+namespace pFlow
+{
+
+
+class AdamsMoulton3
+:
+	public integration
+{
+protected:
+
+	realx3PointField_D& y0_;
+
+	realx3PointField_D& dy0_;
+
+	realx3PointField_D& dy1_;
+
+	using rpIntegration = Kokkos::RangePolicy<
+			DefaultExecutionSpace,
+			Kokkos::Schedule<Kokkos::Static>,
+			Kokkos::IndexType<int32>
+			>;
+public:
+
+	// type info
+	TypeName("AdamsMoulton3");
+
+	//// - Constructors
+		AdamsMoulton3(
+			const word& baseName,
+			repository& owner,
+			const pointStructure& pStruct,
+			const word& method);
+
+		virtual ~AdamsMoulton3()=default;
+
+		// - add a virtual constructor 
+		add_vCtor(
+			integration,
+			AdamsMoulton3,
+			word);
+
+
+	//// - Methods
+		bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+
+		bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override;
+
+		bool setInitialVals(
+			const int32IndexContainer& newIndices,
+			const realx3Vector& y) override;
+
+		bool needSetInitialVals()const override
+		{
+			return true;
+		}
+
+		uniquePtr<integration> clone()const override
+		{
+			return makeUnique<AdamsMoulton3>(*this);
+		}
+
+		bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+
+		template<typename activeFunctor>
+		bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP);
+
+		bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng);
+
+		template<typename activeFunctor>
+		bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP );
+
+};
+
+
+template<typename activeFunctor>
+bool AdamsMoulton3::predictRange(
+	real dt, 
+	realx3Vector_D& y,
+	realx3Vector_D& dy,
+	activeFunctor activeP )
+{
+	auto d_dy = dy.deviceVectorAll();
+	auto d_y  = y.deviceVectorAll();
+	auto d_y0 = y0_.deviceVectorAll();
+	auto d_dy0 = dy0_.deviceVectorAll();
+	auto d_dy1= dy1_.deviceVectorAll();
+	
+	auto activeRng = activeP.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsMoulton3::predictRange",
+		rpIntegration (activeRng.first, activeRng.second),
+		LAMBDA_HD(int32 i){
+			if(activeP(i))
+			{
+				d_dy0[i] = d_dy[i];
+				d_y[i] = d_y0[i] + 
+				dt*
+				(
+					static_cast<real>(3.0 / 2.0) * d_dy[i] 
+					- static_cast<real>(1.0 / 2.0) * d_dy1[i]
+				);
+			}
+		});
+	Kokkos::fence();
+
+	return true;
+
+}
+
+
+template<typename activeFunctor>
+bool pFlow::AdamsMoulton3::intRange(
+	real dt, 
+	realx3Vector_D& y,
+	realx3Vector_D& dy,
+	activeFunctor activeP )
+{
+
+	auto d_dy = dy.deviceVectorAll();
+	auto d_y  = y.deviceVectorAll();
+
+	auto d_dy0 = dy0_.deviceVectorAll();
+	auto d_y0  = y0_.deviceVectorAll();
+	auto d_dy1 = dy1_.deviceVectorAll();
+
+	auto activeRng = activeP.activeRange();
+
+	Kokkos::parallel_for(
+		"AdamsMoulton3::correct",
+		rpIntegration (activeRng.first, activeRng.second),
+		LAMBDA_HD(int32 i){
+			if( activeP(i))
+			{
+				auto corrct_y = d_y0[i] + dt*(
+					  static_cast<real>(5.0/12.0)*d_dy[i] 
+					+ static_cast<real>(8.0/12.0)*d_dy0[i] 
+					- static_cast<real>(1.0/12.0)*d_dy1[i]);
+				d_y[i]  = corrct_y;
+				d_y0[i] = corrct_y;
+				d_dy1[i]= d_dy0[i];
+			}
+		});
+	Kokkos::fence();
+
+
+	return true;
+}
+
+} // pFlow
+
+#endif //__integration_H__
diff --git a/src/Integration/CMakeLists.txt b/src/Integration/CMakeLists.txt
index c96de103..183b0477 100644
--- a/src/Integration/CMakeLists.txt
+++ b/src/Integration/CMakeLists.txt
@@ -5,6 +5,7 @@ AdamsBashforth5/AdamsBashforth5.C
 AdamsBashforth4/AdamsBashforth4.C
 AdamsBashforth3/AdamsBashforth3.C
 AdamsBashforth2/AdamsBashforth2.C
+AdamsMoulton3/AdamsMoulton3.C
 )
 
 set(link_libs Kokkos::kokkos phasicFlow)
diff --git a/src/Integration/integration/integration.H b/src/Integration/integration/integration.H
index 7050499d..d9b2f14d 100644
--- a/src/Integration/integration/integration.H
+++ b/src/Integration/integration/integration.H
@@ -78,6 +78,14 @@ public:
 
 		virtual bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) = 0;
 
+		virtual bool setInitialVals(
+				const int32IndexContainer& newIndices,
+				const realx3Vector& y) = 0;
+
+		virtual bool needSetInitialVals()const = 0;
+
+		virtual uniquePtr<integration> clone()const=0;
+
 		const word& baseName()const
 		{
 			return baseName_;
diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H
index 11a7f90c..276e1ae3 100644
--- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.H
+++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.H
@@ -162,7 +162,9 @@ public:
 		}
 		else
 			return false;
-		
+		//output<<" ContactSearch::broadSearch::PP fater.\n";
+
+		//output<<" ContactSearch::broadSearch::PW before.\n";
 		if(wallMapping_)
 		{
 			sphereWallTimer_.start();
@@ -171,6 +173,8 @@ public:
 		}
 		else
 			return false;
+
+		//output<<" ContactSearch::broadSearch::PW after.\n";
 		
 		return true;
 	}
diff --git a/src/Interaction/contactSearch/methods/NBS.H b/src/Interaction/contactSearch/methods/NBS.H
index 7729f6a8..cf5aed45 100644
--- a/src/Interaction/contactSearch/methods/NBS.H
+++ b/src/Interaction/contactSearch/methods/NBS.H
@@ -302,11 +302,13 @@ public:
 		performedSearch_ = false;
 
 		if( !performSearch() ) return true;
-		
+		//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before build"<<endInfo;
 		build(activeRange);
+		//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) after build"<<endInfo;
 			
+		//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) before findPairs"<<endInfo;
 		findPairs(pairs);
-		
+		//Info<<"NBS::broadSearch(PairsContainer& pairs, range activeRange, bool force=false) after findPairs"<<endInfo;
 		performedSearch_ = true;		
 		return true;
 	}
diff --git a/src/Interaction/contactSearch/methods/NBSLoop.H b/src/Interaction/contactSearch/methods/NBSLoop.H
index 6033ffa5..313832e7 100644
--- a/src/Interaction/contactSearch/methods/NBSLoop.H
+++ b/src/Interaction/contactSearch/methods/NBSLoop.H
@@ -42,6 +42,7 @@ while( m > -1 )
 		{
 			auto ln = n;
 			auto lm = m;
+			//output<<"lm "<< lm <<" ln "<< ln << " i j k "<< int32x3(i,j,k)<<endl;
 			if(lm>ln) Swap(lm,ln);
 			if( auto res = pairs.insert(lm,ln); res <0)
 			{
diff --git a/src/Interaction/sphereInteraction/sphereInteraction.H b/src/Interaction/sphereInteraction/sphereInteraction.H
index 89ca8cef..950276e2 100644
--- a/src/Interaction/sphereInteraction/sphereInteraction.H
+++ b/src/Interaction/sphereInteraction/sphereInteraction.H
@@ -134,19 +134,26 @@ public:
 
 	bool iterate() override
 	{
-		
+
+//Info<<"before contact search"<<endInfo;
 		if(this->contactSearch_)
 		{
+
 			if( this->contactSearch_().ppEnterBroadSearch())
 			{
+				//Info<<" before ppEnterBroadSearch"<<endInfo;
 				ppContactList_().beforeBroadSearch();
+				//Info<<" after ppEnterBroadSearch"<<endInfo;
 			}
 
 			if(this->contactSearch_().pwEnterBroadSearch())
 			{
+				//Info<<" before pwEnterBroadSearch"<<endInfo;
 				pwContactList_().beforeBroadSearch();
+				//Info<<" after pwEnterBroadSearch"<<endInfo;
 			}
 
+			//Info<<" before broadSearch"<<endInfo;
 			if( !contactSearch_().broadSearch(
 					ppContactList_(),
 					pwContactList_()) )
@@ -156,26 +163,37 @@ public:
 				fatalExit;
 			}
 
+			//Info<<" before broadSearch"<<endInfo;
+
+
 			if(this->contactSearch_().ppPerformedBroadSearch())
 			{
+				//Info<<" before afterBroadSearch"<<endInfo;
 				ppContactList_().afterBroadSearch();
+				//Info<<" after afterBroadSearch"<<endInfo;
 			}
 
 			if(this->contactSearch_().pwPerformedBroadSearch())
 			{
+				//Info<<" before pwContactList_().afterBroadSearch()"<<endInfo;
 				pwContactList_().afterBroadSearch();
+				//Info<<" after pwContactList_().afterBroadSearch()"<<endInfo;
 			}
 		}
+//Info<<"after contact search"<<endInfo;
 
 		if( sphParticles_.numActive()<=0)return true;
-
+//Info<<"before sphereSphereInteraction "<<endInfo;
 		ppInteractionTimer_.start();
 			sphereSphereInteraction();
 		ppInteractionTimer_.end();
-		
+//Info<<"after sphereSphereInteraction "<<endInfo;
+
+//Info<<"before sphereWallInteraction "<<endInfo;
 		pwInteractionTimer_.start();
 			sphereWallInteraction();
 		pwInteractionTimer_.end();
+//Info<<"after sphereWallInteraction "<<endInfo;
 		
 		return true;
 	}
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.C b/src/Particles/SphereParticles/sphereParticles/sphereParticles.C
index e3c0ba48..122a10d0 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.C
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.C
@@ -83,8 +83,15 @@ bool pFlow::sphereParticles::beforeIteration()
 	particles::beforeIteration();
 	
 	intPredictTimer_.start();
+	
+	//Info<<"before dyn predict"<<endInfo;
 		dynPointStruct_.predict(this->dt(), accelertion_);
+	//Info<<"after dyn predict"<<endInfo;
+
+	//Info<<"before revel predict"<<endInfo;
 		rVelIntegration_().predict(this->dt(),rVelocity_, rAcceleration_);
+	//Info<<"after rvel predict"<<endInfo;
+
 	intPredictTimer_.end();
 
 	return true;
@@ -93,7 +100,9 @@ bool pFlow::sphereParticles::beforeIteration()
 
 bool pFlow::sphereParticles::iterate() 
 {
+
 	accelerationTimer_.start();
+	//Info<<"before accelerationTimer_ "<<endInfo;
 		pFlow::sphereParticlesKernels::acceleration(
 			control().g(),
 			mass().deviceVectorAll(),
@@ -107,8 +116,11 @@ bool pFlow::sphereParticles::iterate()
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
+	//Info<<"before correct dyn "<<endInfo;
 		dynPointStruct_.correct(this->dt(), accelertion_);
+	//Info<<"after correct dyn "<<endInfo;
 		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
+	//Info<<"after correct rvel "<<endInfo;
 	intCorrectTimer_.end();
 	
 	return true;
@@ -285,6 +297,28 @@ pFlow::sphereParticles::sphereParticles(
 		fatalExit;
 	}
 
+	
+	if(rVelIntegration_->needSetInitialVals())
+	{
+		
+		auto [minInd, maxInd] = pStruct().activeRange();
+		int32IndexContainer indexHD(minInd, maxInd);
+	
+		auto n = indexHD.size();
+		auto index = indexHD.indicesHost();
+
+		realx3Vector rvel(n,RESERVE());
+		const auto hrVel = rVelocity_.hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			rvel.push_back( hrVel[index(i)]);
+		}
+		
+		Report(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endReport;
+		rVelIntegration_->setInitialVals(indexHD, rvel);
+
+	}
 
 	if(!initializeParticles())
 	{
@@ -293,6 +327,32 @@ pFlow::sphereParticles::sphereParticles(
 	
 } 
 
+bool pFlow::sphereParticles::update(const eventMessage& msg) 
+{
+	
+	if(rVelIntegration_->needSetInitialVals())
+	{
+		
+		
+		auto indexHD = pStruct().insertedPointIndex();
+	
+		auto n = indexHD.size();
+		auto index = indexHD.indicesHost();
+
+		realx3Vector rvel(n,RESERVE());
+		const auto hrVel = rVelocity_.hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			rvel.push_back( hrVel[index(i)]);
+		}
+		
+		rVelIntegration_->setInitialVals(indexHD, rvel);
+
+	}	
+
+	return true;
+}
 
 bool pFlow::sphereParticles::insertParticles
 (
@@ -312,7 +372,7 @@ bool pFlow::sphereParticles::insertParticles
 
 	auto exclusionListAllPtr = getFieldObjectList(); 
 
-	auto newInsertedPtr = pStruct().insertPoints( position, exclusionListAllPtr());
+	auto newInsertedPtr = pStruct().insertPoints( position, setField, time(), exclusionListAllPtr());
 	
 	if(!newInsertedPtr)
 	{
@@ -337,10 +397,6 @@ bool pFlow::sphereParticles::insertParticles
 		return false;
 	}
 	
-	for(auto sfEntry:setField)
-	{
-		if(!sfEntry.setPointFieldSelectedAll( time(), newInserted, false ))return false;
-	}
 
 	auto activeR = this->activeRange();
 
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.H b/src/Particles/SphereParticles/sphereParticles/sphereParticles.H
index 0b451ac6..0b208844 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.H
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.H
@@ -157,11 +157,8 @@ public:
 		shapes_.diameterMinMax(minDiam, maxDiam);
 	}
 
-	bool update(const eventMessage& msg) override
-	{
-		CONSUME_PARAM(msg);
-		return true;
-	}
+	bool update(const eventMessage& msg) override;
+	
 	
 	realx3PointField_D& rAcceleration() override
 	{
diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.C b/src/Particles/dynamicPointStructure/dynamicPointStructure.C
index 1828e87a..606ba9e1 100644
--- a/src/Particles/dynamicPointStructure/dynamicPointStructure.C
+++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.C
@@ -47,26 +47,28 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 				objectFile::READ_ALWAYS,
 				objectFile::WRITE_ALWAYS
 				),
-			pStruct_,
+			pStruct(),
 			zero3
 			)
 		)
 
 {
 
+	this->subscribe(pStruct());
+
 	Report(1)<< "Creating integration method "<<
 		greenText(integrationMethod_)<<" for dynamicPointStructure."<<endReport;
 
 	integrationPos_ = integration::create(
 		"pStructPosition",
 		time_.integration(),
-		pStruct_,
+		pStruct(),
 		integrationMethod_);
 	
 	integrationVel_ = integration::create(
 		"pStructVelocity",
 		time_.integration(),
-		pStruct_,
+		pStruct(),
 		integrationMethod_);
 
 	if( !integrationPos_ )
@@ -83,6 +85,36 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		fatalExit;
 	}
 
+	
+	if(!integrationPos_->needSetInitialVals()) return;
+
+	
+		
+	auto [minInd, maxInd] = pStruct().activeRange();
+	int32IndexContainer indexHD(minInd, maxInd);
+	
+	auto n = indexHD.size();
+	auto index = indexHD.indicesHost();
+
+	realx3Vector pos(n,RESERVE());
+	realx3Vector vel(n,RESERVE());
+	const auto hVel = velocity().hostVector();
+	const auto hPos = pStruct().pointPosition().hostVector();
+
+	for(auto i=0; i<n; i++)
+	{
+		pos.push_back( hPos[index(i)]);
+		vel.push_back( hVel[index(i)]);
+	}
+
+	//output<< "pos "<< pos<<endl;
+	//output<< "vel "<< vel<<endl;
+
+	Report(2)<< "Initializing the required vectors for position integratoin "<<endReport;
+	integrationPos_->setInitialVals(indexHD, pos);
+
+	Report(2)<< "Initializing the required vectors for velocity integratoin\n "<<endReport;
+	integrationVel_->setInitialVals(indexHD, vel);
 }
 
 bool pFlow::dynamicPointStructure::predict
@@ -91,7 +123,7 @@ bool pFlow::dynamicPointStructure::predict
 	realx3PointField_D& acceleration
 )
 {
-	auto& pos = pStruct_.pointPosition();
+	auto& pos = pStruct().pointPosition();
 
 	if(!integrationPos_().predict(dt, pos.VectorField(), velocity_.VectorField() ))return false;
 	if(!integrationVel_().predict(dt, velocity_.VectorField(), acceleration.VectorField()))return false;
@@ -105,11 +137,82 @@ bool pFlow::dynamicPointStructure::correct
 	realx3PointField_D& acceleration
 )
 {
-	auto& pos = pStruct_.pointPosition();
+	auto& pos = pStruct().pointPosition();
 	
 	if(!integrationPos_().correct(dt, pos.VectorField(), velocity_.VectorField() ))return false;
 	
 	if(!integrationVel_().correct(dt, velocity_.VectorField(), acceleration.VectorField()))return false;
 
 	return true;	
+}
+
+/*FUNCTION_H
+pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::dynamicPointStructure::insertPoints
+(
+	const realx3Vector& pos,
+	const List<eventObserver*>& exclusionList
+)
+{
+	auto newIndicesPtr = pointStructure::insertPoints(pos, exclusionList);
+
+	// no new point is inserted
+	if( !newIndicesPtr ) return newIndicesPtr;
+
+	if(!integrationPos_().needSetInitialVals()) return newIndicesPtr;
+
+	auto hVel = velocity_.hostVector();
+	auto n = newIndicesPtr().size();
+	auto index = newIndicesPtr().indicesHost();
+
+	realx3Vector velVec(n, RESERVE());
+	for(auto i=0; i<n; i++)
+	{
+		velVec.push_back( hVel[ index(i) ]);
+	}
+
+	integrationPos_->setInitialVals(newIndicesPtr(), pos );
+	integrationVel_->setInitialVals(newIndicesPtr(), velVec );
+	
+	return newIndicesPtr;
+
+}*/
+
+
+bool pFlow::dynamicPointStructure::update(
+	const eventMessage& msg) 
+{
+	if( msg.isInsert())
+	{
+		
+		if(!integrationPos_->needSetInitialVals())return true;
+
+		const auto indexHD = pStruct().insertedPointIndex();
+
+	
+		auto n = indexHD.size();
+
+		if(n==0) return true;
+
+		auto index = indexHD.indicesHost();
+
+		realx3Vector pos(n,RESERVE());
+		realx3Vector vel(n,RESERVE());
+		const auto hVel = velocity().hostVector();
+		const auto hPos = pStruct().pointPosition().hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			pos.push_back( hPos[index(i)]);
+			vel.push_back( hVel[index(i)]);
+		}
+
+
+		
+		integrationPos_->setInitialVals(indexHD, pos);
+
+		integrationVel_->setInitialVals(indexHD, vel);		
+
+	}
+
+	return true;
 }
\ No newline at end of file
diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.H b/src/Particles/dynamicPointStructure/dynamicPointStructure.H
index a62248d6..9a7c8649 100644
--- a/src/Particles/dynamicPointStructure/dynamicPointStructure.H
+++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.H
@@ -31,12 +31,15 @@ namespace pFlow
 {
 
 class dynamicPointStructure
+:
+	//public pointStructure
+public eventObserver
 {
 protected:
 
 	Time& time_;
 
-	const word 			integrationMethod_;
+	word 			integrationMethod_;
 
 	pointStructure& 	pStruct_;
 
@@ -52,8 +55,36 @@ public:
 
 	dynamicPointStructure(Time& time, const word& integrationMethod);
 
+	dynamicPointStructure(const dynamicPointStructure& ps) = default;
+
+	/*dynamicPointStructure(const dynamicPointStructure& ps):
+		pointStructure(ps),
+		time_(ps.time_),
+		integrationMethod_(ps.integrationMethod_),
+		velocity_(ps.velocity_),
+		integrationPos_(ps.integrationPos_->clone()),
+		integrationVel_(ps.integrationVel_->clone())
+	{
+
+	}*/
+
+		
+		// - no move construct 
+	dynamicPointStructure(dynamicPointStructure&&) = delete;
+
+	// - copy assignment 
+	//
+	//   should be changed, may causs undefined behavior 
+	//////////////////////////////////////////////////////////////
+	dynamicPointStructure& operator=(const dynamicPointStructure&) = default;
+
+	// - no move assignment 
+	dynamicPointStructure& operator=(dynamicPointStructure&&) = delete;
+
+	// - destructor 
 	virtual ~dynamicPointStructure() = default;
 
+
 	inline pointStructure& pStruct() 
 	{
 		return pStruct_;
@@ -79,6 +110,19 @@ public:
 	bool correct(real dt, realx3PointField_D& acceleration);
 
 
+	// - update data structure by inserting/setting new points 
+	//   Notifies all the fields in the registered list of data structure
+	//   and exclude the fields that re in the exclusionList
+	//   retrun nullptr if it fails 
+	/*FUNCTION_H
+	virtual uniquePtr<int32IndexContainer> insertPoints(
+		const realx3Vector& pos,
+		const List<eventObserver*>& exclusionList={nullptr}
+	)override;*/
+
+	bool update(const eventMessage& msg) override;
+	
+
 };
 
 }
diff --git a/src/Particles/particles/particles.C b/src/Particles/particles/particles.C
index 391c8797..ae8c3a68 100644
--- a/src/Particles/particles/particles.C
+++ b/src/Particles/particles/particles.C
@@ -31,10 +31,21 @@ pFlow::particles::particles
 	demParticles(control),
 	time_(control.time()),
 	integrationMethod_(integrationMethod),
+	/*dynPointStruct_(
+		time_.emplaceObject<dynamicPointStructure>(
+			objectFile(
+				pointStructureFile__,
+				"",
+				objectFile::READ_ALWAYS,
+				objectFile::WRITE_ALWAYS			
+				),
+			control.time(),
+			integrationMethod
+			)
+		),*/
 	dynPointStruct_(
 		control.time(),
-		integrationMethod
-		),
+		integrationMethod),
 	shapeName_(
 		control.time().emplaceObject<wordPointField>(
 			objectFile(
diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.H b/src/phasicFlow/containers/VectorHD/VectorSingle.H
index b2898e29..80bdbd81 100644
--- a/src/phasicFlow/containers/VectorHD/VectorSingle.H
+++ b/src/phasicFlow/containers/VectorHD/VectorSingle.H
@@ -559,6 +559,8 @@ public:
 		INLINE_FUNCTION_H
 		bool insertSetElement(const int32IndexContainer& indices, const Vector<T>& vals)
 		{
+
+			//Info<<"start of insertSetElement vecotsingle"<<endInfo;
 			if(indices.size() == 0)return true;
 			if(indices.size() != vals.size())return false;
 
@@ -569,38 +571,27 @@ public:
 				resize(maxInd+1);
 			} 
 
-			if constexpr (isHostAccessible_)
-			{
-				// TODO: remove const_cast
-				hostViewType1D<T> dVecVals( const_cast<T*>(vals.data()), vals.size());
+				
+			hostViewType1D<const T> hVecVals( vals.data(), vals.size());
+			deviceViewType1D<T> dVecVals("dVecVals", indices.size());
 
+			copy(dVecVals, hVecVals);
 
-				pFlow::algorithms::KOKKOS::fillSelected<T, int32, execution_space>(
-					deviceVectorAll().data(),
-					indices.hostView().data(),
-					dVecVals.data(),
-					indices.size());
+			using policy = Kokkos::RangePolicy<
+			execution_space,
+			Kokkos::IndexType<int32> >;
+			auto dVec = deviceVectorAll();
+			auto dIndex = indices.deviceView();
+		
+			Kokkos::parallel_for(
+				"insertSetElement",
+				policy(0,indices.size()), LAMBDA_HD(int32 i){	
+				dVec(dIndex(i))= dVecVals(i);
+				});
+			Kokkos::fence();
 
-				return true;
+			return true;
 			
-			}else
-			{
-				
-				// TODO: remove const_cast
-				hostViewType1D<T> hVecVals( const_cast<T*>(vals.data()), vals.size());
-				deviceViewType1D<T> dVecVals("dVecVals", indices.size());
-				copy(dVecVals, hVecVals);
-
-				
-				pFlow::algorithms::KOKKOS::fillSelected<T, int32, execution_space>(
-					deviceVectorAll().data(),
-					indices.deviceView().data(),
-					dVecVals.data(),
-					indices.size());
-				return true;
-			}
-
-			return false;
 		}
 
 		INLINE_FUNCTION_H
diff --git a/src/phasicFlow/setFieldList/setFieldEntry.C b/src/phasicFlow/setFieldList/setFieldEntry.C
index 7e1c4222..8ceb9b14 100644
--- a/src/phasicFlow/setFieldList/setFieldEntry.C
+++ b/src/phasicFlow/setFieldList/setFieldEntry.C
@@ -57,36 +57,50 @@ bool pFlow::setFieldEntry::checkForTypeAndValueAll()const
 	return true;
 }
 
-bool pFlow::setFieldEntry::setPointFieldDefaultValueNewAll
+void* pFlow::setFieldEntry::setPointFieldDefaultValueNewAll
 (
 	repository& owner,
 	pointStructure& pStruct,
 	bool verbose
 )
 {
-	if( 
-		!(
-			setPointFieldDefaultValueNew<int8> (owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<int16>(owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<int32> (owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<int64> (owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<uint32> (owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<label>(owner, pStruct, verbose) ||
-			setPointFieldDefaultValueNew<real>(owner, pStruct, verbose)  ||
-			setPointFieldDefaultValueNew<realx3>(owner, pStruct, verbose) ||
-			setPointFieldDefaultValueStdNew<word>(owner, pStruct, verbose)
-		 )
-	   )
+	if( void* res =  setPointFieldDefaultValueNew<int8> (owner, pStruct, verbose) ; res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<int16>(owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<int32> (owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<int64> (owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<uint32> (owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<label>(owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<real>(owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueNew<realx3>(owner, pStruct, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldDefaultValueStdNew<word>(owner, pStruct, verbose); res)
+	{
+		return res;
+	}else
 	{
 		fatalErrorInFunction<< 
 		"  un-supported data type "<<entry_.firstPart() << " in setField for field " << fieldName() <<endl;
-		return false;
+		return nullptr;
 	}
 
-	return true;
 }
 
-bool pFlow::setFieldEntry::setPointFieldSelectedAll
+void* pFlow::setFieldEntry::setPointFieldSelectedAll
 (
 	repository& owner,
 	int32IndexContainer& selected,
@@ -94,25 +108,38 @@ bool pFlow::setFieldEntry::setPointFieldSelectedAll
 )
 {
 
-	if( 
-		!(
-			setPointFieldSelected<int8> (owner, selected, verbose) ||
-			setPointFieldSelected<int16>(owner, selected, verbose) ||
-			setPointFieldSelected<int32> (owner, selected, verbose) ||
-			setPointFieldSelected<int64> (owner, selected, verbose) ||
-			setPointFieldSelected<uint32> (owner, selected, verbose) ||
-			setPointFieldSelected<label>(owner, selected, verbose) ||
-			setPointFieldSelected<real>(owner, selected, verbose)  ||
-			setPointFieldSelected<realx3>(owner, selected, verbose) ||
-			setPointFieldSelectedStd<word>(owner, selected, verbose)
-		 )
-	   )
+	if( void* res =  setPointFieldSelected<int8> (owner, selected, verbose) ; res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<int16>(owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<int32> (owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<int64> (owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<uint32> (owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<label>(owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<real>(owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelected<realx3>(owner, selected, verbose); res)
+	{
+		return res;
+	}else if(void* res = setPointFieldSelectedStd<word>(owner, selected, verbose); res)
+	{
+		return res;
+	}else
 	{
 		fatalErrorInFunction<< 
 		"  un-supported data type "<<entry_.firstPart() << " in setField for field " << fieldName() <<endl;
-		return false;
+		return nullptr;
 	}
 
-	return true;
-
 }
diff --git a/src/phasicFlow/setFieldList/setFieldEntry.H b/src/phasicFlow/setFieldList/setFieldEntry.H
index 62f16112..07fdae6a 100644
--- a/src/phasicFlow/setFieldList/setFieldEntry.H
+++ b/src/phasicFlow/setFieldList/setFieldEntry.H
@@ -65,7 +65,7 @@ public:
 		bool checkForTypeAndValueAll()const;
 		
 		template <typename Type>
-		bool setPointFieldDefaultValueNew
+		void* setPointFieldDefaultValueNew
 		(
 			repository& owner,
 			pointStructure& pStruct,
@@ -73,14 +73,14 @@ public:
 		);
 
 		template<typename Type>
-		bool setPointFieldDefaultValueStdNew
+		void* setPointFieldDefaultValueStdNew
 		(
 			repository& owner,
 			pointStructure& pStruct,
 			bool verbose = false
 		);
 
-		bool setPointFieldDefaultValueNewAll
+		void* setPointFieldDefaultValueNewAll
 		(
 			repository& owner,
 			pointStructure& pStruct,
@@ -88,7 +88,7 @@ public:
 		);
 
 		template <typename Type>
-		bool setPointFieldSelected
+		void* setPointFieldSelected
 		(
 			repository& owner,
 			int32IndexContainer& selected,
@@ -96,14 +96,14 @@ public:
 		);
 
 		template <typename Type>
-		bool setPointFieldSelectedStd
+		void* setPointFieldSelectedStd
 		(
 			repository& owner,
 			int32IndexContainer& selected,
 			bool verbose = false
 		);
 
-		bool setPointFieldSelectedAll
+		void* setPointFieldSelectedAll
 		(
 			repository& owner,
 			int32IndexContainer& selected,
diff --git a/src/phasicFlow/setFieldList/setFieldEntryTemplates.C b/src/phasicFlow/setFieldList/setFieldEntryTemplates.C
index 99ed6745..d3bda996 100644
--- a/src/phasicFlow/setFieldList/setFieldEntryTemplates.C
+++ b/src/phasicFlow/setFieldList/setFieldEntryTemplates.C
@@ -37,7 +37,7 @@ bool pFlow::setFieldEntry::checkForTypeAndValue()const
 }
 
 template <typename Type>
-bool pFlow::setFieldEntry::setPointFieldDefaultValueNew
+void* pFlow::setFieldEntry::setPointFieldDefaultValueNew
 (
 	repository& owner,
 	pointStructure& pStruct,
@@ -45,7 +45,7 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueNew
 )
 {
 
-	if( !checkForType<Type>() ) return false;
+	if( !checkForType<Type>() ) return nullptr;
 	
 	Type defValue = entry_.secondPartVal<Type>();
 	
@@ -54,6 +54,7 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueNew
 	             " in repository "<< owner.name() <<endReport;
 
 	
+	auto& field = 
 	owner.emplaceObject<pointField<VectorSingle,Type>>
 	(
 		objectFile
@@ -67,11 +68,11 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueNew
 		defValue
 	);
 
-	return true;
+	return &field;
 }
 
 template <typename Type>
-bool pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
+void* pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
 (
 	repository& owner,
 	pointStructure& pStruct,
@@ -79,7 +80,7 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
 )
 {
 
-	if( !checkForType<Type>() ) return false;
+	if( !checkForType<Type>() ) return nullptr;
 	
 	Type defValue = entry_.secondPartVal<Type>();
 	
@@ -88,6 +89,7 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
 	             " in repository "<< owner.name() <<endReport;
 
 	// by default we perform operations on host
+	auto& field = 
 	owner.emplaceObject<pointField<Vector,Type, vecAllocator<Type>>>
 	(
 		objectFile
@@ -101,18 +103,18 @@ bool pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
 		defValue
 	);
 
-	return true;
+	return &field;
 }
 
 template <typename Type>
-bool pFlow::setFieldEntry::setPointFieldSelected
+void* pFlow::setFieldEntry::setPointFieldSelected
 (
 	repository& owner,
 	int32IndexContainer& selected,
 	bool verbose
 )
 {
-	if( !checkForType<Type>() ) return false;
+	if( !checkForType<Type>() ) return nullptr;
 
 	
 	auto fName = fieldName();
@@ -121,7 +123,7 @@ bool pFlow::setFieldEntry::setPointFieldSelected
 	{
 		fatalErrorInFunction<<
 		"Cannot find "<< fName << " in repository " << owner.name() << ". \n";
-		return false;
+		return nullptr;
 	}
 
 	Type value = entry_.secondPartVal<Type>();
@@ -137,31 +139,40 @@ bool pFlow::setFieldEntry::setPointFieldSelected
 	{
 		
 		auto& field = owner.lookupObject<pointField<VectorSingle,Type>>(fName);
-		return field.insertSetElement(selected, value);	
+		if(field.insertSetElement(selected, value))
+			return &field;
+		else
+			return nullptr;
 	}
 
 	if( pointField<VectorSingle,Type, HostSpace>::TYPENAME() == fieldTypeName )
 	{
 		
 		auto& field = owner.lookupObject<pointField<VectorSingle,Type,HostSpace>>(fName);
-		return field.insertSetElement(selected, value);	
+		if(field.insertSetElement(selected, value))
+			return &field;
+		else
+			return nullptr;
 	}
 
 	if( pointField<VectorDual,Type>::TYPENAME() == fieldTypeName )
 	{
 		
 		auto& field = owner.lookupObject<pointField<VectorDual,Type>>(fName);
-		return field.insertSetElement(selected, value);	
+		if(field.insertSetElement(selected, value))
+			return &field;
+		else
+			return nullptr;
 	}
 
 	fatalErrorInFunction<<
 	fieldTypeName<< " is not a supported field type for setFieldEntry.\n";
-	return false;
+	return nullptr;
 	
 }
 
 template <typename Type>
-bool pFlow::setFieldEntry::setPointFieldSelectedStd
+void* pFlow::setFieldEntry::setPointFieldSelectedStd
 (
 	repository& owner,
 	int32IndexContainer& selected,
@@ -169,7 +180,7 @@ bool pFlow::setFieldEntry::setPointFieldSelectedStd
 )
 {
 	
-	if( !checkForType<Type>() ) return false;
+	if( !checkForType<Type>() ) return nullptr;
 
 	
 	auto fName = fieldName();
@@ -178,7 +189,7 @@ bool pFlow::setFieldEntry::setPointFieldSelectedStd
 	{
 		fatalErrorInFunction<<
 		"  Cannot find "<< fName << " in repository " << owner.name() << ". \n";
-		return false;
+		return nullptr;
 	}
 
 	
@@ -194,9 +205,11 @@ bool pFlow::setFieldEntry::setPointFieldSelectedStd
 	if( pointField<Vector, Type, vecAllocator<Type>>::TYPENAME() == fieldTypeName )
 	{
 		auto& field = owner.lookupObject<pointField<Vector,Type, vecAllocator<Type>>>(fName);
-		return field.insertSetElement(selected, value);	
+		if(field.insertSetElement(selected, value))
+			return &field;
+		else
+			return nullptr;
 	}
 
-
-	return false;
+	return nullptr;
 }
\ No newline at end of file
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.C b/src/phasicFlow/structuredData/pointStructure/pointStructure.C
index 9126c384..37c19f7b 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.C
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.C
@@ -21,6 +21,7 @@ Licence:
 
 #include "pointStructure.H"
 #include "pointStructureKernels.H"
+#include "setFieldList.H"
 #include "error.H"
 #include "iOstream.H"
 #include "Time.H"
@@ -298,10 +299,13 @@ FUNCTION_H
 pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
 (
 	const realx3Vector& pos,
+	const setFieldList& setField,
+	repository& owner,
 	const List<eventObserver*>& exclusionList
 )
 {
 
+
 	auto numNew = pos.size();
 	if( numNew==0)
 	{
@@ -312,27 +316,53 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
 
 	if(!newPointsPtr)return nullptr;
 
-	auto oldSize = size();
-	auto oldCapacity = capacity();
-	auto oldRange = activeRange();
+	auto oldSize 			= size();
+	auto oldCapacity 	= capacity();
+	auto oldRange 		= activeRange();
 
 	tobeInsertedIndex_ = newPointsPtr();
-	
+
 	// set the position of new points
-	if(!pointPosition_.insertSetElement(newPointsPtr(), pos))return nullptr;
-	
-	if(!pointFlag_.insertSetElement( newPointsPtr(), static_cast<int8>(PointFlag::ACTIVE)))return nullptr;
-	
+
+	if(!pointPosition_.insertSetElement(
+		newPointsPtr(),
+		pos)
+		)return nullptr;
+
+	if(!pointFlag_.insertSetElement(
+		newPointsPtr(),
+		static_cast<int8>(PointFlag::ACTIVE))
+		)return nullptr;
+
+
 
 	setNumMaxPoints();
 	auto minInd = newPointsPtr().min();
 	auto maxInd = newPointsPtr().max();
 
+
+	List<eventObserver*> localExlusion(exclusionList);
+
+	for(auto sfEntry:setField)
+	{
+		if(void* fieldPtr = 
+			sfEntry.setPointFieldSelectedAll( 
+				owner, 
+				newPointsPtr(), 
+				false );
+			fieldPtr)
+
+			localExlusion.push_back( 
+				static_cast<eventObserver*>(fieldPtr)
+				);
+		else
+			return nullptr;
+	}
+
 	// changes the active rage based on the new inserted points 
 	activeRange_ = { min(activeRange_.first, minInd ),
 					max(activeRange_.second, maxInd+1)};
 
-
 	numActivePoints_ += numNew;
 
 	eventMessage msg(eventMessage::INSERT);
@@ -347,7 +377,7 @@ pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::pointStructure::insertPoints
 		msg.add(eventMessage::CAP_CHANGED);
 	
 	// notify all the registered objects except the exclusionList 
-	if( !this->notify(msg, exclusionList) ) return nullptr;
+	if( !this->notify(msg, localExlusion) ) return nullptr;
 
 	return newPointsPtr;
 }
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.H b/src/phasicFlow/structuredData/pointStructure/pointStructure.H
index 9d53890c..1cf74045 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.H
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.H
@@ -38,6 +38,8 @@ namespace pFlow
 
 //forward
 class box;
+class setFieldList;
+class repository;
 
 class pointStructure
 :
@@ -319,6 +321,8 @@ public:
 		FUNCTION_H
 		virtual uniquePtr<int32IndexContainer> insertPoints(
 			const realx3Vector& pos,
+			const setFieldList& setField,
+			repository& owner,
 			const List<eventObserver*>& exclusionList={nullptr}
 		);