From 853d50e96f446ac5badb468221472121f621086d Mon Sep 17 00:00:00 2001
From: hamidrezanorouzi <hamid.r.norouzi@gmail.com>
Date: Sat, 24 Dec 2022 15:00:00 +0330
Subject: [PATCH] correction in the macros to be compatible with OpenFOAM

---
 DEMSystems/CMakeLists.txt                     |   2 +
 DEMSystems/DEMSystem/DEMSystem.cpp            |  75 ++++++++++
 DEMSystems/DEMSystem/DEMSystem.hpp            | 128 ++++++++++++++++++
 .../domainDistribute/domainDistribute.cpp     | 101 ++++++++++++++
 .../domainDistribute/domainDistribute.hpp     |  89 ++++++++++++
 .../sphereDEMSystem/sphereDEMSystem.cpp       |  78 ++++++++---
 .../sphereDEMSystem/sphereDEMSystem.hpp       |  76 +++++------
 .../sphereGranFlow/createDEMComponents.hpp    |   6 +-
 solvers/sphereGranFlow/sphereGranFlow.cpp     |   4 +-
 src/Geometry/geometry/geometry.cpp            |  10 +-
 .../geometryMotion/geometryMotion.cpp         |   4 +-
 .../Models/contactForce/linearCF.hpp          |   4 +-
 .../Models/contactForce/nonLinearCF.hpp       |   4 +-
 .../Models/contactForce/nonLinearMod.hpp      |   2 +-
 .../ContactSearch/ContactSearch.hpp           |   8 +-
 .../contactSearch/contactSearch.cpp           |   4 +-
 .../contactSearch/methods/NBSLevel0.hpp       |   4 +-
 .../contactSearch/methods/NBSLevels.hpp       |  10 +-
 .../wallMappings/cellsWallLevel0.hpp          |   4 +-
 .../wallMappings/cellsWallLevels.hpp          |   4 +-
 .../wallMappings/multiGridMapping.hpp         |   4 +-
 src/Interaction/interaction/interaction.cpp   |   4 +-
 .../sphereInteraction/sphereInteraction.cpp   |   2 +-
 .../rotatingAxisMotion/rotatingAxisMotion.cpp |   2 +-
 .../Insertion/Insertion/Insertion.cpp         |  20 +--
 .../InsertionRegion/InsertionRegion.cpp       |   2 +-
 .../Insertion/insertion/insertion.cpp         |   8 +-
 .../Insertion/shapeMixture/shapeMixture.cpp   |   6 +-
 .../sphereParticles/sphereParticles.cpp       |  28 ++--
 .../dynamicPointStructure.cpp                 |   8 +-
 src/Particles/particles/particleIdHandler.hpp |   2 +-
 src/phasicFlow/Kokkos/KokkosTypes.hpp         |   2 +-
 src/phasicFlow/containers/Vector/Vector.cpp   |   8 +-
 .../containers/Vector/VectorAlgorithm.hpp     |   2 +-
 .../pointField/pointFieldAlgorithms.hpp       |  14 +-
 src/phasicFlow/globals/error.hpp              |   4 +-
 src/phasicFlow/globals/pFlowMacros.hpp        |   2 +-
 src/phasicFlow/repository/Time/Time.cpp       |   2 +-
 .../repository/Time/timeControl.hpp           |   2 +-
 .../repository/repository/repository.cpp      |   2 +-
 .../setFieldList/setFieldEntryTemplates.cpp   |  16 +--
 src/phasicFlow/streams/streams.hpp            |  17 +--
 src/phasicFlow/structuredData/box/box.hpp     |   5 +
 .../pointStructure/pointStructure.cpp         |   4 +-
 .../selectors/selectBox/selectBox.cpp         |   2 +-
 .../trisurfaceStructure/multiTriSurface.cpp   |   4 +-
 src/setHelpers/initialize.hpp                 |   6 +-
 src/setHelpers/initialize_Control.hpp         |   2 +-
 src/setHelpers/setProperty.hpp                |   2 +-
 src/setHelpers/setSurfaceGeometry.hpp         |   2 +-
 utilities/checkPhasicFlow/checkPhasicFlow.cpp |   4 +-
 .../geometryPhasicFlow/geometryPhasicFlow.cpp |  16 +--
 utilities/pFlowToVTK/geometric.hpp            |   2 +-
 utilities/pFlowToVTK/pointFieldToVTK.hpp      |  16 +--
 utilities/pFlowToVTK/triSurfaceFieldToVTK.hpp |   4 +-
 .../particlesPhasicFlow.cpp                   |  22 +--
 .../positionParticles/positionParticles.cpp   |   2 +-
 .../positionRandom/positionRandom.cpp         |  10 +-
 .../postprocessPhasicFlow/includeMask.cpp     |   2 +-
 .../postprocessPhasicFlow/postprocess.cpp     |  16 +--
 .../postprocessPhasicFlow/processField.cpp    |   2 +-
 61 files changed, 670 insertions(+), 227 deletions(-)
 create mode 100644 DEMSystems/DEMSystem/DEMSystem.cpp
 create mode 100644 DEMSystems/DEMSystem/DEMSystem.hpp
 create mode 100644 DEMSystems/domainDistribute/domainDistribute.cpp
 create mode 100644 DEMSystems/domainDistribute/domainDistribute.hpp

diff --git a/DEMSystems/CMakeLists.txt b/DEMSystems/CMakeLists.txt
index 156e23bd..243df700 100644
--- a/DEMSystems/CMakeLists.txt
+++ b/DEMSystems/CMakeLists.txt
@@ -1,6 +1,8 @@
 
 set(SourceFiles 
+DEMSystem/DEMSystem.cpp
 sphereDEMSystem/sphereDEMSystem.cpp
+domainDistribute/domainDistribute.cpp
 )
 
 set(link_libs Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities)
diff --git a/DEMSystems/DEMSystem/DEMSystem.cpp b/DEMSystems/DEMSystem/DEMSystem.cpp
new file mode 100644
index 00000000..9c9b0c02
--- /dev/null
+++ b/DEMSystems/DEMSystem/DEMSystem.cpp
@@ -0,0 +1,75 @@
+/*------------------------------- 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 "DEMSystem.hpp"
+
+
+pFlow::coupling::DEMSystem::DEMSystem(
+	word  demSystemName,
+	int32 numDomains, 
+	const std::vector<box>& domains, 
+	int argc, 
+	char* argv[])
+:
+	ControlDict_(),
+	domains_(domains)
+{
+  
+	REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
+	Control_ = makeUnique<systemControl>(
+		ControlDict_.startTime(),
+		ControlDict_.endTime(),
+		ControlDict_.saveInterval(),
+		ControlDict_.startTimeName());
+
+}
+
+pFlow::coupling::DEMSystem::~DEMSystem()
+{}
+
+
+pFlow::uniquePtr<pFlow::coupling::DEMSystem>
+	pFlow::coupling::DEMSystem::create(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains, 
+		int argc, 
+		char* argv[]
+		)
+{
+	if( wordvCtorSelector_.search(demSystemName) )
+	{
+		return wordvCtorSelector_[demSystemName] (demSystemName, numDomains, domains, argc, argv);
+	}
+	else
+	{
+		printKeys
+		( 
+			fatalError << "Ctor Selector "<< demSystemName << " dose not exist. \n"
+			<<"Avaiable ones are: \n\n"
+			,
+			wordvCtorSelector_
+		);
+		return nullptr;
+	}
+
+	return nullptr;
+}
+
diff --git a/DEMSystems/DEMSystem/DEMSystem.hpp b/DEMSystems/DEMSystem/DEMSystem.hpp
new file mode 100644
index 00000000..a2731261
--- /dev/null
+++ b/DEMSystems/DEMSystem/DEMSystem.hpp
@@ -0,0 +1,128 @@
+/*------------------------------- 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 __DEMSystem_hpp__
+#define __DEMSystem_hpp__
+
+#include <vector>
+
+#include "types.hpp"
+#include "virtualConstructor.hpp"
+#include "uniquePtr.hpp"
+#include "systemControl.hpp"
+#include "readControlDict.hpp"
+
+
+namespace pFlow::coupling
+{
+
+
+class DEMSystem
+{
+protected:
+	
+	readControlDict 				ControlDict_;
+
+	uniquePtr<systemControl>		Control_ = nullptr;
+
+	std::vector<box> 				domains_;
+
+
+	// methods 
+	auto& Control()
+	{
+		return Control_();
+	}
+
+
+public:
+
+	TypeInfo("DEMSystem");
+
+	DEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains, 
+		int argc, 
+		char* argv[]);
+
+	virtual ~DEMSystem();
+
+	DEMSystem(const DEMSystem&)=delete;
+
+	DEMSystem& operator = (const DEMSystem&)=delete;
+
+	create_vCtor(
+		DEMSystem,
+		word,
+		(
+			word  demSystemName,
+			int32 numDomains, 
+			const std::vector<box>& domains, 
+			int argc, 
+			char* argv[]
+		),
+		(
+			demSystemName,
+			numDomains, 
+			domains, 
+			argc, 
+			argv
+		));
+
+	realx3 g()const
+	{
+		return Control_->g();
+	}
+
+	auto inline constexpr usingDoulle()const
+	{
+		return pFlow::usingDouble__;
+	}
+
+	virtual 
+	int32 numParInDomain(int32 di)const = 0;
+	
+	virtual
+	std::vector<int32> numParInDomain()const = 0;
+
+	virtual 
+	bool iterate(int32 n, real timeToWrite, word timeName) = 0;
+	
+	virtual 
+	real maxBounndingSphereSize()const = 0;
+
+	static 
+	uniquePtr<DEMSystem>
+		create(
+			word  demSystemName,
+			int32 numDomains, 
+			const std::vector<box>& domains, 
+			int argc, 
+			char* argv[]);
+
+
+
+};
+
+
+} // pFlow
+
+#endif // __DEMSystem_hpp__
\ No newline at end of file
diff --git a/DEMSystems/domainDistribute/domainDistribute.cpp b/DEMSystems/domainDistribute/domainDistribute.cpp
new file mode 100644
index 00000000..18e43832
--- /dev/null
+++ b/DEMSystems/domainDistribute/domainDistribute.cpp
@@ -0,0 +1,101 @@
+/*------------------------------- 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 "domainDistribute.hpp"
+
+
+
+pFlow::coupling::domainDistribute::domainDistribute(
+	int32 numDomains, 
+	const Vector<box>& domains_,
+	real maxBoundingBox)
+:
+numDomains_(numDomains),
+extDomains_("extDomains", numDomains),
+particlesInDomains_("particlesInDomains", numDomains),
+numParInDomain_("numParInDomain", numDomains, 0),
+maxBoundingBoxSize_(maxBoundingBox)
+{
+
+	realx3 dl = domainExtension_ * maxBoundingBoxSize_;
+
+	for(int32 i=0; i<numDomains_; i++)
+	{
+		extDomains_[i] = extendBox(domains_[i], dl);
+	}
+}
+
+bool pFlow::coupling::domainDistribute::locateParticles(
+		ViewType1D<realx3,HostSpace> points, includeMask mask)
+	{
+		range active = mask.activeRange();
+		auto numInDomain = numParInDomain_.deviceVectorAll();
+		auto numDomains = numDomains_;
+
+		using policy = Kokkos::RangePolicy<
+			DefaultHostExecutionSpace,
+			Kokkos::IndexType<int32> >;
+
+		Kokkos::parallel_for("locateParticles",
+			policy(active.first, active.second),
+			LAMBDA_HD(int32 i){
+				if(mask(i))
+				{
+					for(int32 di=0; di<numDomains; di++)
+					{
+						if(extDomains_[i].isInside(points[i]))
+							Kokkos::atomic_add(&numInDomain[di],1);
+					}
+				}
+			});
+		Kokkos::fence();
+
+		
+		for(int32 i=0; i<numDomains_; i++)
+		{
+			particlesInDomains_[i].resize(numParInDomain_[i]);
+		}
+
+		numParInDomain_.fill(0);
+
+		auto particlesInDomainsPtr = particlesInDomains_.data();
+
+		Kokkos::parallel_for("locateParticles",
+			policy(active.first, active.second),
+			LAMBDA_HD(int32 i){
+				if(mask(i))
+				{
+					for(int32 di=0; di<numDomains; di++)
+					{
+						if(extDomains_[i].isInside(points[i]))
+						{
+
+							particlesInDomainsPtr[di][numInDomain[di]] = i;
+							Kokkos::atomic_add(&numInDomain[di],1);
+
+						}
+					}
+				}
+			});
+		Kokkos::fence();		
+
+
+		return true;
+	}
\ No newline at end of file
diff --git a/DEMSystems/domainDistribute/domainDistribute.hpp b/DEMSystems/domainDistribute/domainDistribute.hpp
new file mode 100644
index 00000000..0095feef
--- /dev/null
+++ b/DEMSystems/domainDistribute/domainDistribute.hpp
@@ -0,0 +1,89 @@
+/*------------------------------- 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 __domainDistribute_hpp__
+#define __domainDistribute_hpp__
+
+#include "box.hpp"
+#include "Vectors.hpp"
+#include "VectorSingles.hpp"
+#include "pointStructure.hpp"
+
+namespace pFlow::coupling
+{
+
+class domainDistribute
+{
+protected:
+	
+	// protected members 
+
+	int32 									numDomains_;
+
+	VectorSingle<box,HostSpace> 			extDomains_;
+
+	Vector<VectorSingle<int32,HostSpace>>	particlesInDomains_;
+	
+	int32Vector_H 							numParInDomain_;
+
+
+	real maxBoundingBoxSize_;
+
+	int32 lastUpdateIter_ = 0;
+
+	int32 updateInterval_ = 1;
+
+	real domainExtension_ = 1.0;
+
+	using includeMask = typename pointStructure::activePointsHost;
+
+public:
+
+	domainDistribute(
+		int32 numDomains, 
+		const Vector<box>& domains_,
+		real maxBoundinBox);
+
+	~domainDistribute()=default;
+
+	domainDistribute(const domainDistribute&)=delete;
+
+	domainDistribute& operator=(const domainDistribute&)=delete;
+
+	int32 numDomains()const
+	{
+		return numDomains_;
+	}
+	
+	int32 numParInDomain(int32 di)const
+	{
+		return numParInDomain_[di];
+	}
+
+	//template<typename includeMask>
+	bool locateParticles(
+		ViewType1D<realx3,HostSpace> points, includeMask mask);
+	
+};
+
+
+} // pFlow::coupling
+
+#endif //__domainDistribute_hpp__
\ No newline at end of file
diff --git a/DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp b/DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp
index 9d145494..14139f18 100644
--- a/DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp
+++ b/DEMSystems/sphereDEMSystem/sphereDEMSystem.cpp
@@ -21,37 +21,37 @@ Licence:
 #include "sphereDEMSystem.hpp"
 
 
-pFlow::sphereDEMSystem::sphereDEMSystem(int argc, char* argv[])
+pFlow::coupling::sphereDEMSystem::sphereDEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[])
 :
-	ControlDict_()
+	DEMSystem(demSystemName, numDomains, domains, argc, argv)
 {
   
-	Report(0)<<"Initializing host/device execution spaces . . . \n";
-	Report(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endReport;
-	Report(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endReport;
+	REPORT(0)<<"Initializing host/device execution spaces . . . \n";
+	REPORT(1)<<"Host execution space is "<< greenText(DefaultHostExecutionSpace::name())<<endREPORT;
+	REPORT(1)<<"Device execution space is "<<greenText(DefaultExecutionSpace::name())<<endREPORT;
 
  	// initialize Kokkos
 	Kokkos::initialize( argc, argv );  
   
-	Report(0)<<"\nCreating Control repository . . ."<<endReport;
-	Control_ = makeUnique<systemControl>(
-		ControlDict_.startTime(),
-		ControlDict_.endTime(),
-		ControlDict_.saveInterval(),
-		ControlDict_.startTimeName());
+	
 
-	Report(0)<<"\nReading proprties . . . "<<endReport;
+	REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
 	property_ = makeUnique<property>(
     	Control().caseSetup().path()+propertyFile__);
 
-	Report(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endReport;
+	REPORT(0)<< "\nCreating surface geometry for sphereDEMSystem . . . "<<endREPORT;
 	geometry_ = geometry::create(Control(), Property());
 
-	Report(0)<<"\nReading sphere particles . . ."<<endReport;
+	REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
 	particles_ = makeUnique<sphereParticles>(Control(), Property());
 
 
-	//Report(0)<<"\nCreating particle insertion for spheres. . ."<<endReport;
+	//REPORT(0)<<"\nCreating particle insertion for spheres. . ."<<endREPORT;
 	/*insertion_ = 
 	Control().caseSetup().emplaceObject<sphereInsertion>(
 		objectFile(
@@ -64,15 +64,20 @@ pFlow::sphereDEMSystem::sphereDEMSystem(int argc, char* argv[])
 		sphParticles.shapes()
 	);*/
 
-	Report(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endReport;
+	REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
 	interaction_ = interaction::create(
 		Control(),
 		Particles(),
 		Geometry());
 
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	particleDistribution_ = makeUnique<domainDistribute>(numDomains, domains, maxD);
+
 }
 
-pFlow::sphereDEMSystem::~sphereDEMSystem()
+pFlow::coupling::sphereDEMSystem::~sphereDEMSystem()
 {
 	interaction_.reset();
 	insertion_.reset();
@@ -85,3 +90,42 @@ pFlow::sphereDEMSystem::~sphereDEMSystem()
 	Kokkos::finalize();
 }
 
+
+
+pFlow::int32 
+	pFlow::coupling::sphereDEMSystem::numParInDomain(int32 di)const
+{
+	return particleDistribution_().numParInDomain(di);
+}
+
+std::vector<pFlow::int32> 
+	pFlow::coupling::sphereDEMSystem::numParInDomain()const
+{
+	const auto& distribute = particleDistribution_();
+	int32 numDomains = distribute.numDomains();
+	std::vector<int32> nums(numDomains);
+	for(int32 i=0; i<numDomains; i++)
+	{
+		nums[i] = distribute.numParInDomain(i);
+	}
+
+	return nums;
+}
+
+
+bool pFlow::coupling::sphereDEMSystem::iterate(
+	int32 n, 
+	real timeToWrite, 
+	word timeName) 
+{
+	return true;
+}
+
+pFlow::real 
+	pFlow::coupling::sphereDEMSystem::maxBounndingSphereSize()const
+{
+	real minD, maxD;
+	particles_->boundingSphereMinMax(minD, maxD);
+
+	return maxD;
+}
\ No newline at end of file
diff --git a/DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp b/DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp
index e6233bdd..5ddc04e5 100644
--- a/DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp
+++ b/DEMSystems/sphereDEMSystem/sphereDEMSystem.hpp
@@ -21,54 +21,40 @@ Licence:
 #ifndef __sphereDEMSystem_hpp__
 #define __sphereDEMSystem_hpp__
 
-#include <array>
-
-#include "systemControl.hpp"
+#include "DEMSystem.hpp"
 #include "property.hpp"
 #include "uniquePtr.hpp"
 #include "geometry.hpp"
 #include "sphereParticles.hpp"
 #include "interaction.hpp"
 #include "Insertions.hpp"
-#include "readControlDict.hpp"
+#include "domainDistribute.hpp"
 
 
-
-
-namespace pFlow
+namespace pFlow::coupling
 {
 
 class sphereDEMSystem
+:
+	public DEMSystem
 {
 protected:
 	
-	readControlDict 						ControlDict_;
+	// protected members 
 
-	uniquePtr<systemControl>		Control_ = nullptr;
+	uniquePtr<property> 		property_ = nullptr;
 
-	uniquePtr<property> 				property_ = nullptr;
-
-	uniquePtr<geometry> 			  geometry_ = nullptr;
+	uniquePtr<geometry> 		geometry_ = nullptr;
 
 	uniquePtr<sphereParticles> 	particles_ = nullptr;
 
 	uniquePtr<sphereInsertion> 	insertion_ = nullptr;
 
-	uniquePtr<interaction> 			interaction_ = nullptr;
+	uniquePtr<interaction> 		interaction_ = nullptr;
 
+	uniquePtr<domainDistribute> particleDistribution_=nullptr;
 
-	int32 									numDomains_;
-
-	VectorDual<box> 				domains_;
-
-	Vector<int32Vector_H> 	parIndexInDomain_;
-
-	
-	auto& Control()
-	{
-		return Control_();
-	}
-
+//  protected member functions
 	auto& Property()
 	{
 		return property_();
@@ -91,27 +77,39 @@ protected:
 
 public:
 
-	sphereDEMSystem(int argc, char* argv[]);
+	TypeInfo("sphereDEMSystem");
 
-	~sphereDEMSystem();
+	sphereDEMSystem(
+		word  demSystemName,
+		int32 numDomains, 
+		const std::vector<box>& domains,
+		int argc, 
+		char* argv[]);
+
+	virtual ~sphereDEMSystem();
 
 	sphereDEMSystem(const sphereDEMSystem&)=delete;
 
 	sphereDEMSystem& operator = (const sphereDEMSystem&)=delete;
 
-	
-	std::array<double,3> g()const
-	{
-		return {
-			Control_->g().x(), 
-			Control_->g().y(), 
-			Control_->g().z()};
-	}
+	add_vCtor(
+		DEMSystem,
+		sphereDEMSystem,
+		word);
 
-	bool inline usingDoulle()const
-	{
-		return pFlow::usingDouble__;
-	}
+
+	 
+	int32 numParInDomain(int32 di)const override;
+
+	
+	std::vector<int32> numParInDomain()const override;
+	
+	
+	virtual 
+	bool iterate(int32 n, real timeToWrite, word timeName) override;
+	
+	virtual 
+	real maxBounndingSphereSize()const override;
 
 };
 
diff --git a/solvers/sphereGranFlow/createDEMComponents.hpp b/solvers/sphereGranFlow/createDEMComponents.hpp
index a407a9e1..b4ade9c8 100755
--- a/solvers/sphereGranFlow/createDEMComponents.hpp
+++ b/solvers/sphereGranFlow/createDEMComponents.hpp
@@ -19,11 +19,11 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 //
-Report(0)<<"\nReading sphere particles . . ."<<endReport;
+REPORT(0)<<"\nReading sphere particles . . ."<<endREPORT;
 sphereParticles sphParticles(Control, proprties);
 
 //
-Report(0)<<"\nCreating particle insertion object . . ."<<endReport;
+REPORT(0)<<"\nCreating particle insertion object . . ."<<endREPORT;
 auto& sphInsertion = 
 	Control.caseSetup().emplaceObject<sphereInsertion>(
 		objectFile(
@@ -38,7 +38,7 @@ auto& sphInsertion =
 
 
 
-Report(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endReport;
+REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."<<endREPORT;
 auto interactionPtr = interaction::create(
 	Control,
 	sphParticles,
diff --git a/solvers/sphereGranFlow/sphereGranFlow.cpp b/solvers/sphereGranFlow/sphereGranFlow.cpp
index e00020fc..a4920b86 100755
--- a/solvers/sphereGranFlow/sphereGranFlow.cpp
+++ b/solvers/sphereGranFlow/sphereGranFlow.cpp
@@ -61,7 +61,7 @@ if(!cmds.parse(argc, argv)) return 0;
 
 	#include "createDEMComponents.hpp"
 	
-	Report(0)<<"\nStart of time loop . . .\n"<<endReport;
+	REPORT(0)<<"\nStart of time loop . . .\n"<<endREPORT;
 
 	do 
 	{
@@ -98,7 +98,7 @@ if(!cmds.parse(argc, argv)) return 0;
 
 	}while(Control++);
 
-	Report(0)<<"\nEnd of time loop.\n"<<endReport;
+	REPORT(0)<<"\nEnd of time loop.\n"<<endREPORT;
 
 // this should be palced in each main 
 #include "finalize.hpp"
diff --git a/src/Geometry/geometry/geometry.cpp b/src/Geometry/geometry/geometry.cpp
index 03d7fd4d..d663ce78 100644
--- a/src/Geometry/geometry/geometry.cpp
+++ b/src/Geometry/geometry/geometry.cpp
@@ -29,7 +29,7 @@ bool pFlow::geometry::findPropertyId()
 	propId.clear();
 
 	uint32 pId;
-	forAll(matI, materialName_)
+	ForAll(matI, materialName_)
 	{
 		
 		if( !wallProperty_.nameToIndex( materialName_[matI], pId ) )
@@ -249,11 +249,11 @@ pFlow::uniquePtr<pFlow::geometry>
 	
 	auto geomModel = angleBracketsNames("geometry", model);
 
-	Report(1)<< "Selecting geometry model . . ."<<endReport;
+	REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
 	if( systemControlvCtorSelector_.search(geomModel) )
 	{
 		auto objPtr = systemControlvCtorSelector_[geomModel] (control, prop);
-		Report(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endReport;
+		REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
 		return objPtr;
 	}
 	else
@@ -285,7 +285,7 @@ pFlow::uniquePtr<pFlow::geometry>
 
 	auto geomModel = angleBracketsNames("geometry", model);
 
-	Report(1)<< "Selecting geometry model . . ."<<endReport;
+	REPORT(1)<< "Selecting geometry model . . ."<<endREPORT;
 
 	if( dictionaryvCtorSelector_.search(geomModel) )
 	{
@@ -298,7 +298,7 @@ pFlow::uniquePtr<pFlow::geometry>
 				motionCompName,
 				propName
 			);
-		Report(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endReport;
+		REPORT(2)<<"Model "<< greenText(geomModel)<<" is created.\n"<<endREPORT;
 		return objPtr;
 	}
 	else
diff --git a/src/Geometry/geometryMotion/geometryMotion.cpp b/src/Geometry/geometryMotion/geometryMotion.cpp
index 5e6136f5..20fd0cd7 100644
--- a/src/Geometry/geometryMotion/geometryMotion.cpp
+++ b/src/Geometry/geometryMotion/geometryMotion.cpp
@@ -53,7 +53,7 @@ bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 	triMotionIndex_.reserve( this->surface().capacity() );
 	triMotionIndex_.clear();
 
-	forAll( surfI, motionComponentName_)
+	ForAll( surfI, motionComponentName_)
 	{
 		auto mName = motionComponentName_[surfI];
 		auto mInd = motionModel_.nameToIndex(mName);
@@ -71,7 +71,7 @@ bool pFlow::geometryMotion<MotionModel>::findMotionIndex()
 	pointMotionIndex_.reserve(triSurface_.numPoints());
 	pointMotionIndex_.clear();
 
-	forAll(surfI, motionIndex_)
+	ForAll(surfI, motionIndex_)
 	{
 		auto nP = triSurface_.surfNumPoints(surfI);
 		for(int32 i=0; i<nP; i++)
diff --git a/src/Interaction/Models/contactForce/linearCF.hpp b/src/Interaction/Models/contactForce/linearCF.hpp
index 6f70489f..5ef6da3e 100644
--- a/src/Interaction/Models/contactForce/linearCF.hpp
+++ b/src/Interaction/Models/contactForce/linearCF.hpp
@@ -135,7 +135,7 @@ protected:
 
 		realVector etha_n(nElem);
 		realVector etha_t(nElem);
-		forAll(i , kn)
+		ForAll(i , kn)
 		{
 			etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/
 					sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0));
@@ -145,7 +145,7 @@ protected:
 		}
 
 		Vector<linearProperties> prop(nElem);
-		forAll(i,kn)
+		ForAll(i,kn)
 		{
 			prop[i] = {kn[i], kt[i], etha_n[i], etha_t[i], mu[i]};
 		}
diff --git a/src/Interaction/Models/contactForce/nonLinearCF.hpp b/src/Interaction/Models/contactForce/nonLinearCF.hpp
index 5e7f27f3..b9f42d69 100644
--- a/src/Interaction/Models/contactForce/nonLinearCF.hpp
+++ b/src/Interaction/Models/contactForce/nonLinearCF.hpp
@@ -123,7 +123,7 @@ protected:
 
 		realVector etha_n(nElem);
 		
-		forAll(i , en)
+		ForAll(i , en)
 		{
 			//K_hertz = 4.0/3.0*Yeff*sqrt(Reff); 
 			//-2.2664*log(en)*sqrt(meff*K_hertz)/sqrt( log(en)**2 + 10.1354);
@@ -138,7 +138,7 @@ protected:
 		}
 
 		Vector<nonLinearProperties> prop(nElem);
-		forAll(i,Yeff)
+		ForAll(i,Yeff)
 		{
 			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
 		}
diff --git a/src/Interaction/Models/contactForce/nonLinearMod.hpp b/src/Interaction/Models/contactForce/nonLinearMod.hpp
index c9dbdfa9..3bfbd29b 100644
--- a/src/Interaction/Models/contactForce/nonLinearMod.hpp
+++ b/src/Interaction/Models/contactForce/nonLinearMod.hpp
@@ -122,7 +122,7 @@ protected:
 
 
 		Vector<nonLinearProperties> prop(nElem);
-		forAll(i,Yeff)
+		ForAll(i,Yeff)
 		{
 			prop[i] = {Yeff[i], Geff[i], etha_n[i], mu[i]};
 		}
diff --git a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
index d8065d94..d6428064 100644
--- a/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
+++ b/src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
@@ -98,8 +98,8 @@ public:
 				position,
 				diam
 			);
-		Report(2)<<"Contact search algorithm for particle-particle is "<<
-				 greenText(particleContactSearch_().typeName())<<endReport;
+		REPORT(2)<<"Contact search algorithm for particle-particle is "<<
+				 greenText(particleContactSearch_().typeName())<<endREPORT;
 
 
 		auto wmDict = dict().subDict(wmMethod+"Info");
@@ -120,8 +120,8 @@ public:
 				wPoints,
 				wVertices
 				);
-		Report(2)<<"Wall mapping algorithm for particle-wall is "<<
-				 greenText(wallMapping_().typeName())<< endReport;
+		REPORT(2)<<"Wall mapping algorithm for particle-wall is "<<
+				 greenText(wallMapping_().typeName())<< endREPORT;
 
 	}
 
diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp
index b26215b4..17c1dbd3 100644
--- a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp
+++ b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp
@@ -53,13 +53,13 @@ pFlow::uniquePtr<pFlow::contactSearch> pFlow::contactSearch::create(
 	auto model = angleBracketsNames2("ContactSearch", baseMethName, wallMethod);
 
 
-	Report(1)<<"Selecting contact search model . . ."<<endReport;
+	REPORT(1)<<"Selecting contact search model . . ."<<endREPORT;
 	
 
 	if( dictionaryvCtorSelector_.search(model))
 	{
 		auto objPtr = dictionaryvCtorSelector_[model] (dict, domain, prtcl, geom, timers);
-		Report(2)<<"Model "<< greenText(model)<<" is created."<<endReport;
+		REPORT(2)<<"Model "<< greenText(model)<<" is created."<<endREPORT;
 		return objPtr;
 	}
 	else
diff --git a/src/Interaction/contactSearch/methods/NBSLevel0.hpp b/src/Interaction/contactSearch/methods/NBSLevel0.hpp
index 2e23f961..d9ae2f83 100644
--- a/src/Interaction/contactSearch/methods/NBSLevel0.hpp
+++ b/src/Interaction/contactSearch/methods/NBSLevel0.hpp
@@ -200,8 +200,8 @@ public:
 				
 				pairs.increaseCapacityBy(len);
 
-				Info<< "The contact pair container capacity increased from "<<
-				oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<endInfo;
+				INFORMATION<< "The contact pair container capacity increased from "<<
+				oldCap << " to "<<pairs.capacity()<<" in NBSLevel0."<<endINFO;
 				
 			}
 
diff --git a/src/Interaction/contactSearch/methods/NBSLevels.hpp b/src/Interaction/contactSearch/methods/NBSLevels.hpp
index ca7bd730..30603bd3 100644
--- a/src/Interaction/contactSearch/methods/NBSLevels.hpp
+++ b/src/Interaction/contactSearch/methods/NBSLevels.hpp
@@ -108,11 +108,11 @@ protected:
 		copy(sizeRangeLevels_, sizeRangeLevelsHost_);
 		copy(maxSizeLevels_, maxSizeLevelsHost_);
 
-		Report(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<<endReport;
+		REPORT(2)<<"Grids with "<< yellowText(numLevels_)<< " levels have been created."<<endREPORT;
 		for(int32 lvl=0; lvl<numLevels_; lvl++)
 		{
-			Report(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
-				<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endReport;
+			REPORT(3)<<"Cell gird No "<< yellowText(lvl)<<" with size range ("
+				<<sizeRangeLevelsHost_[lvl].first<<","<<sizeRangeLevelsHost_[lvl].second<<"]."<<endREPORT;
 		}
 		
 		return true;
@@ -235,8 +235,8 @@ public:
 				
 				pairs.increaseCapacityBy(len);
 
-				Info<< "The contact pair container capacity increased from "<<
-				oldCap << " to "<<pairs.capacity()<<" in NBSLevels."<<endInfo;
+				INFORMATION<< "The contact pair container capacity increased from "<<
+				oldCap << " to "<<pairs.capacity()<<" in NBSLevels."<<endINFO;
 				
 			}
 
diff --git a/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp b/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp
index 9ea9ce67..65034507 100644
--- a/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp
+++ b/src/Interaction/contactSearch/wallMappings/cellsWallLevel0.hpp
@@ -196,9 +196,9 @@ public:
 				auto oldCap = pairs.capacity();
 				pairs.increaseCapacityBy(len);
 
-				Info<<"Contact pair container capacity increased from "<<
+				INFORMATION<<"Contact pair container capacity increased from "<<
 						oldCap << " to " 
-						<< pairs.capacity() <<" in cellsWallLevel0."<<endInfo;
+						<< pairs.capacity() <<" in cellsWallLevel0."<<endINFO;
 				
 				Kokkos::fence();
 			}
diff --git a/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp b/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp
index 56f39461..9f2c80ec 100644
--- a/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp
+++ b/src/Interaction/contactSearch/wallMappings/cellsWallLevels.hpp
@@ -132,9 +132,9 @@ public:
 				auto oldCap = pairs.capacity();
 				pairs.increaseCapacityBy(len);
 
-				Info<<"Contact pair container capacity increased from "<<
+				INFORMATION<<"Contact pair container capacity increased from "<<
 						oldCap << " to " 
-						<< pairs.capacity() <<" in cellsWallLevels."<<endInfo;
+						<< pairs.capacity() <<" in cellsWallLevels."<<endINFO;
 				
 				Kokkos::fence();
 			}
diff --git a/src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp b/src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp
index 10a55001..3a266f1f 100644
--- a/src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp
+++ b/src/Interaction/contactSearch/wallMappings/multiGridMapping.hpp
@@ -115,8 +115,8 @@ public:
 			)
 	{
 
-		Report(3)<<"Multi-grid wall mapping with "<< 
-			yellowText(numLevels)<<" levels has been created."<<endReport;
+		REPORT(3)<<"Multi-grid wall mapping with "<< 
+			yellowText(numLevels)<<" levels has been created."<<endREPORT;
 	}
 
 
diff --git a/src/Interaction/interaction/interaction.cpp b/src/Interaction/interaction/interaction.cpp
index 63e797e7..8d368c29 100644
--- a/src/Interaction/interaction/interaction.cpp
+++ b/src/Interaction/interaction/interaction.cpp
@@ -79,14 +79,14 @@ pFlow::uniquePtr<pFlow::interaction> pFlow::interaction::create
 			clType);
 
 
-	Report(1)<< "Selecting interaction model..."<<endReport;
+	REPORT(1)<< "Selecting interaction model..."<<endREPORT;
 	if( systemControlvCtorSelector_.search(interactionModel) )
 	{
 		auto objPtr = 
 			systemControlvCtorSelector_[interactionModel]
 			(control, prtcl, geom);
 
-		Report(2)<<"Model "<<greenText(interactionModel)<<" is created."<<endReport;
+		REPORT(2)<<"Model "<<greenText(interactionModel)<<" is created."<<endREPORT;
 		return objPtr;
 	}
 	else
diff --git a/src/Interaction/sphereInteraction/sphereInteraction.cpp b/src/Interaction/sphereInteraction/sphereInteraction.cpp
index eee24bda..9426041b 100644
--- a/src/Interaction/sphereInteraction/sphereInteraction.cpp
+++ b/src/Interaction/sphereInteraction/sphereInteraction.cpp
@@ -30,7 +30,7 @@ bool pFlow::sphereInteraction<contactForceModel,geometryMotionModel, contactList
 
 	auto modelDict = this->fileDict().subDict("model");
 
-	Report(1)<<"Createing contact force model . . ."<<endReport;
+	REPORT(1)<<"Createing contact force model . . ."<<endREPORT;
 	forceModel_ = makeUnique<ContactForceModel>(
 		this->numMaterials(),
 		rhoD.deviceVector(),
diff --git a/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.cpp b/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.cpp
index 8d8fa3c7..9ef21e44 100644
--- a/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.cpp
+++ b/src/MotionModel/rotatingAxisMotion/rotatingAxisMotion.cpp
@@ -96,7 +96,7 @@ bool pFlow::rotatingAxisMotion::writeDictionary
 
 	auto& motionInfo = dict.subDictOrCreate("rotatingAxisMotionInfo");
 	
-	forAll(i, axis_)
+	ForAll(i, axis_)
 	{
 		
 		auto& axDict = motionInfo.subDictOrCreate(axisName_[i]);
diff --git a/src/Particles/Insertion/Insertion/Insertion.cpp b/src/Particles/Insertion/Insertion/Insertion.cpp
index 08a7f4b9..a8a406e3 100644
--- a/src/Particles/Insertion/Insertion/Insertion.cpp
+++ b/src/Particles/Insertion/Insertion/Insertion.cpp
@@ -37,7 +37,7 @@ bool pFlow::Insertion<ShapeType>::readInsertionDict
 
 	for(auto& name:regionDicNames)
 	{
-		Report(2)<<"reading insertion region "<< greenText(name)<<endReport;
+		REPORT(2)<<"reading insertion region "<< greenText(name)<<endREPORT;
 		regions_.push_backSafe(dict.subDict(name), shapes_);
 	}
 
@@ -54,7 +54,7 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
 
 	if( !this->isActive() ) return true;
 
-	forAll(i,regions_)
+	ForAll(i,regions_)
 	{
 		auto& rgnDict = dict.subDictOrCreate(regions_[i].name());
 
@@ -92,7 +92,7 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 	if(!isActive()) return true;
 	
 
-	forAll(i,regions_)
+	ForAll(i,regions_)
 	{
 		bool insertionOccured = false;
 		auto& rgn = regions_[i];
@@ -106,9 +106,9 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 
 				if(insertionOccured)
 				{
-					Report(0)<<"\nParticle insertion from "<< greenText(rgn.name())<<endReport;
-					Report(1)<< cyanText(pos.size()) << " new particles is being inserted at Time: "<<
-						cyanText(currentTime) <<" s."<<endReport;
+					REPORT(0)<<"\nParticle insertion from "<< greenText(rgn.name())<<endREPORT;
+					REPORT(1)<< cyanText(pos.size()) << " new particles is being inserted at Time: "<<
+						cyanText(currentTime) <<" s."<<endREPORT;
 						
 					if(!particles_.insertParticles(pos, shapes, rgn.setFields()))
 					{
@@ -117,8 +117,8 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 						" to particles. \n";
 						return false;
 					}
-					Report(1)<<"Total number of particles inserted from this region is "<<
-						cyanText(rgn.totalInserted())<<'\n'<<endReport;
+					REPORT(1)<<"Total number of particles inserted from this region is "<<
+						cyanText(rgn.totalInserted())<<'\n'<<endREPORT;
 				}
 				else
 				{
@@ -130,8 +130,8 @@ bool pFlow::Insertion<ShapeType>::insertParticles
 			{
 				if(insertionOccured)
 				{
-					Warning<< "\n fewer number of particles are inserted from region "<< rgn.name() <<
-					" than expected. You may stop the simulation to change settings."<<endWarning;
+					yWARNING<< "\n fewer number of particles are inserted from region "<< rgn.name() <<
+					" than expected. You may stop the simulation to change settings."<<endyWARNING;
 					continue;
 				}
 				else
diff --git a/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp b/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp
index 5af830c9..57404bf7 100644
--- a/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp
+++ b/src/Particles/Insertion/InsertionRegion/InsertionRegion.cpp
@@ -28,7 +28,7 @@ bool pFlow::InsertionRegion<ShapeType>::checkForContact
 )
 {
 
-	forAll(i, pos)
+	ForAll(i, pos)
 	{
 		if( length(pos[i]-p) < 0.5*(diams[i]+d) ) return true;
 	}
diff --git a/src/Particles/Insertion/insertion/insertion.cpp b/src/Particles/Insertion/insertion/insertion.cpp
index 441eb9ab..5737d521 100644
--- a/src/Particles/Insertion/insertion/insertion.cpp
+++ b/src/Particles/Insertion/insertion/insertion.cpp
@@ -33,11 +33,11 @@ bool pFlow::insertion::readInsertionDict
 	active_ = dict.getVal<Logical>("active");
 
 	if(active_)
-		Report(1)<< "Particle insertion mechanism is "<<
-		yellowText("active")<<" in the simulation."<<endReport;
+		REPORT(1)<< "Particle insertion mechanism is "<<
+		yellowText("active")<<" in the simulation."<<endREPORT;
 	else
-		Report(1)<< "Particle insertion mechanism is "<<
-		yellowText("not active")<<" in the simulation."<<endReport;
+		REPORT(1)<< "Particle insertion mechanism is "<<
+		yellowText("not active")<<" in the simulation."<<endREPORT;
 
 	
 	return true;
diff --git a/src/Particles/Insertion/shapeMixture/shapeMixture.cpp b/src/Particles/Insertion/shapeMixture/shapeMixture.cpp
index 3f04b069..041f71dc 100755
--- a/src/Particles/Insertion/shapeMixture/shapeMixture.cpp
+++ b/src/Particles/Insertion/shapeMixture/shapeMixture.cpp
@@ -37,7 +37,7 @@ pFlow::shapeMixture::shapeMixture
 pFlow::word pFlow::shapeMixture::getNextShapeName()
 {
 	
-	forAll(i, names_)
+	ForAll(i, names_)
 	{
 		if(current_[i]< number_[i])
 		{
@@ -116,7 +116,7 @@ bool pFlow::shapeMixture::read(const dictionary& dict)
     }
 
     current_.clear();
-    forAll(i, numberInserted_)
+    ForAll(i, numberInserted_)
     {
     	current_.push_back(numberInserted_[i]%number_[i]);
     }
@@ -130,7 +130,7 @@ bool pFlow::shapeMixture::write
 ) const
 {
 
-	forAll(i, names_)
+	ForAll(i, names_)
 	{
 		if(!dict.add(names_[i],number_[i]))
 		{
diff --git a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
index 13cbe20a..6191a22e 100644
--- a/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
+++ b/src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp
@@ -84,13 +84,13 @@ bool pFlow::sphereParticles::beforeIteration()
 	
 	intPredictTimer_.start();
 	
-	//Info<<"before dyn predict"<<endInfo;
+	//INFO<<"before dyn predict"<<endINFO;
 		dynPointStruct_.predict(this->dt(), accelertion_);
-	//Info<<"after dyn predict"<<endInfo;
+	//INFO<<"after dyn predict"<<endINFO;
 
-	//Info<<"before revel predict"<<endInfo;
+	//INFO<<"before revel predict"<<endINFO;
 		rVelIntegration_().predict(this->dt(),rVelocity_, rAcceleration_);
-	//Info<<"after rvel predict"<<endInfo;
+	//INFO<<"after rvel predict"<<endINFO;
 
 	intPredictTimer_.end();
 
@@ -102,7 +102,7 @@ bool pFlow::sphereParticles::iterate()
 {
 
 	accelerationTimer_.start();
-	//Info<<"before accelerationTimer_ "<<endInfo;
+	//INFO<<"before accelerationTimer_ "<<endINFO;
 		pFlow::sphereParticlesKernels::acceleration(
 			control().g(),
 			mass().deviceVectorAll(),
@@ -116,11 +116,11 @@ bool pFlow::sphereParticles::iterate()
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
-	//Info<<"before correct dyn "<<endInfo;
+	//INFO<<"before correct dyn "<<endINFO;
 		dynPointStruct_.correct(this->dt(), accelertion_);
-	//Info<<"after correct dyn "<<endInfo;
+	//INFO<<"after correct dyn "<<endINFO;
 		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
-	//Info<<"after correct rvel "<<endInfo;
+	//INFO<<"after correct rvel "<<endINFO;
 	intCorrectTimer_.end();
 	
 	return true;
@@ -158,7 +158,7 @@ bool pFlow::sphereParticles::insertSphereParticles(
 	real d, m, I;
 	int8 pId;
 
-	forAll(i, names )
+	ForAll(i, names )
 	{
 		if (diameterMassInertiaPropId(names[i], d, m, I, pId))
 		{
@@ -280,8 +280,8 @@ pFlow::sphereParticles::sphereParticles(
 
 {
 
-	Report(1)<<"Creating integration method "<<greenText(this->integrationMethod())
-		  << " for rotational velocity."<<endReport;
+	REPORT(1)<<"Creating integration method "<<greenText(this->integrationMethod())
+		  << " for rotational velocity."<<endREPORT;
 		  
 	rVelIntegration_ = 
 		integration::create(
@@ -315,7 +315,7 @@ pFlow::sphereParticles::sphereParticles(
 			rvel.push_back( hrVel[index(i)]);
 		}
 		
-		Report(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endReport;
+		REPORT(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endREPORT;
 		rVelIntegration_->setInitialVals(indexHD, rvel);
 
 	}
@@ -400,9 +400,9 @@ bool pFlow::sphereParticles::insertParticles
 
 	auto activeR = this->activeRange();
 
-	Report(1)<< "Active range is "<<yellowText("["<<activeR.first<<", "<<activeR.second<<")")<<
+	REPORT(1)<< "Active range is "<<yellowText("["<<activeR.first<<", "<<activeR.second<<")")<<
 	" and number of active points is "<< cyanText(this->numActive())<<
-	" and pointStructure capacity is "<<cyanText(this->capacity())<<endReport;
+	" and pointStructure capacity is "<<cyanText(this->capacity())<<endREPORT;
 
 	return true;
 		
diff --git a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
index 7f2bba65..16e5b849 100644
--- a/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
+++ b/src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
@@ -56,8 +56,8 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 
 	this->subscribe(pStruct());
 
-	Report(1)<< "Creating integration method "<<
-		greenText(integrationMethod_)<<" for dynamicPointStructure."<<endReport;
+	REPORT(1)<< "Creating integration method "<<
+		greenText(integrationMethod_)<<" for dynamicPointStructure."<<endREPORT;
 
 	integrationPos_ = integration::create(
 		"pStructPosition",
@@ -110,10 +110,10 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 	//output<< "pos "<< pos<<endl;
 	//output<< "vel "<< vel<<endl;
 
-	Report(2)<< "Initializing the required vectors for position integratoin "<<endReport;
+	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;
+	REPORT(2)<< "Initializing the required vectors for velocity integratoin\n "<<endREPORT;
 	integrationVel_->setInitialVals(indexHD, vel);
 }
 
diff --git a/src/Particles/particles/particleIdHandler.hpp b/src/Particles/particles/particleIdHandler.hpp
index 1e2e6d33..abe708b4 100644
--- a/src/Particles/particles/particleIdHandler.hpp
+++ b/src/Particles/particles/particleIdHandler.hpp
@@ -45,7 +45,7 @@ public:
 			nextId_ = 0; 
 			id.modifyOnHost();
 
-			forAll(i,id)
+			ForAll(i,id)
 			{
 				if(id.isActive(i))
 				{
diff --git a/src/phasicFlow/Kokkos/KokkosTypes.hpp b/src/phasicFlow/Kokkos/KokkosTypes.hpp
index 74bca64f..b378124d 100644
--- a/src/phasicFlow/Kokkos/KokkosTypes.hpp
+++ b/src/phasicFlow/Kokkos/KokkosTypes.hpp
@@ -62,7 +62,7 @@ template<typename T, typename... properties>
 	using ViewType1D = Kokkos::View<T*,properties...>;
 
 template<typename T, typename... properties>
-	using DualViewType1D = Kokkos::View<T*,properties...>;
+	using DualViewType1D = Kokkos::DualView<T*,properties...>;
 
 template<typename T, typename... properties>
 	using ViewType3D = Kokkos::View<T***,properties...>;
diff --git a/src/phasicFlow/containers/Vector/Vector.cpp b/src/phasicFlow/containers/Vector/Vector.cpp
index 80c942c6..f644fc25 100644
--- a/src/phasicFlow/containers/Vector/Vector.cpp
+++ b/src/phasicFlow/containers/Vector/Vector.cpp
@@ -208,7 +208,7 @@ bool pFlow::Vector<T, Allocator>::insertSetElement(
 	
 	auto hIndices = indices.hostView();
 
-	forAll(i, indices)
+	ForAll(i, indices)
 	{
 		auto s = size();
 		auto idx = hIndices[i];
@@ -240,7 +240,7 @@ bool pFlow::Vector<T, Allocator>::insertSetElement(
 
 	auto hIndices = indices.hostView();
 	
-	forAll(i, indices)
+	ForAll(i, indices)
 	{
 		auto s = size();
 		auto idx = hIndices[i];
@@ -271,7 +271,7 @@ bool pFlow::Vector<T, Allocator>::insertSetElement
 {
 	if(indices.size() == 0)return true;
 	
-	forAll(i, indices)
+	ForAll(i, indices)
 	{
 		auto s = size();
 		auto idx = indices[i];
@@ -303,7 +303,7 @@ bool pFlow::Vector<T, Allocator>::insertSetElement
 	if(indices.size() == 0)return true;
 	if(indices.size() != vals.size())return false;
 
-	forAll(i, indices)
+	ForAll(i, indices)
 	{
 		auto s = size();
 		auto idx = indices[i];
diff --git a/src/phasicFlow/containers/Vector/VectorAlgorithm.hpp b/src/phasicFlow/containers/Vector/VectorAlgorithm.hpp
index 80c5b9be..c6dd9bd1 100644
--- a/src/phasicFlow/containers/Vector/VectorAlgorithm.hpp
+++ b/src/phasicFlow/containers/Vector/VectorAlgorithm.hpp
@@ -68,7 +68,7 @@ inline void sort(Vector<T, Allocator>& vec)
 template<typename T, typename Allocator>
 inline int64 find(Vector<T, Allocator>& vec, const T& val)
 {
-	forAll( i,vec)
+	ForAll( i,vec)
 	{
 		if ( vec[i] == val) return static_cast<int64>(i);
 	}
diff --git a/src/phasicFlow/containers/pointField/pointFieldAlgorithms.hpp b/src/phasicFlow/containers/pointField/pointFieldAlgorithms.hpp
index 2d38e1bf..022e35f1 100644
--- a/src/phasicFlow/containers/pointField/pointFieldAlgorithms.hpp
+++ b/src/phasicFlow/containers/pointField/pointFieldAlgorithms.hpp
@@ -181,7 +181,7 @@ T maxActive(const pointField<VectorDual, T, MemorySpace>& pField)
 		fill_n(field, field.size(), val);
 		return;
 	}
-	forAll(i, field)
+	ForAll(i, field)
 	{
 		if(field.pStruct().isActive(i)) field[i] = val;
 	}
@@ -191,7 +191,7 @@ template<typename T>
 	void inline fillMarkedDelete( pointField<T>& field, const T& val)
 {
 	if(field.pStruct().allActive())return;
-	forAll(i,field)
+	ForAll(i,field)
 	{
 		if(!field.pStruct().isActive(i)) field[i] = val;
 	}
@@ -215,14 +215,14 @@ inline auto for_eachActive(pointField<T>& field, UnaryPredicate func)
 {
 	if(field.pStruct().allActive())
 	{
-		forAll(i,field)
+		ForAll(i,field)
 		{
 			func(i);
 		}
 	}
 	else
 	{
-		forAll(i, field)
+		ForAll(i, field)
 		{
 			if(field.pStruct().isActive(i)) func(i);
 		}
@@ -235,14 +235,14 @@ inline bool for_eachActiveBreak(pointField<T>& field, UnaryPredicate func)
 {
 	if(field.pStruct().allActive())
 	{
-		forAll(i,field)
+		ForAll(i,field)
 		{
 			if(!func(i))return false;
 		}
 	}
 	else
 	{
-		forAll(i, field)
+		ForAll(i, field)
 		{
 			if(field.pStruct().isActive(i))
 			{
@@ -258,7 +258,7 @@ inline auto for_eachMarkedDelete(pointField<T>& field, UnaryPredicate func)
 {
 	if(field.pStruct().allActive()) return func;
 		
-	forAll(i, field)
+	ForAll(i, field)
 	{
 		if(!field.pStruct().isActive(i)) func(i);
 	}
diff --git a/src/phasicFlow/globals/error.hpp b/src/phasicFlow/globals/error.hpp
index 0a375725..52f56157 100644
--- a/src/phasicFlow/globals/error.hpp
+++ b/src/phasicFlow/globals/error.hpp
@@ -41,10 +41,10 @@ namespace pFlow
 
 #define fatalErrorInFunction fatalErrorIn(FUNCTION_NAME)
 
-#define notImplemented(functionName) \
+#define Not_Implemented(functionName) \
     notImplementedErrorMessage ((functionName), __FILE__, __LINE__ )
         
-#define notImplementedFunction notImplemented(FUNCTION_NAME);
+#define notImplementedFunction Not_Implemented(FUNCTION_NAME);
 
 #define ioErrorInFile( fileName, lineNumber) \
 	ioErrorMessage( fileName, lineNumber, FUNCTION_NAME, __FILE__, __LINE__ )
diff --git a/src/phasicFlow/globals/pFlowMacros.hpp b/src/phasicFlow/globals/pFlowMacros.hpp
index f0f78229..fd3a5bf6 100755
--- a/src/phasicFlow/globals/pFlowMacros.hpp
+++ b/src/phasicFlow/globals/pFlowMacros.hpp
@@ -68,7 +68,7 @@ Licence:
 #endif
 
 
-#define forAll(i, container) for(auto i=0; i < container.size(); ++i)
+#define ForAll(i, container) for(auto i=0; i < container.size(); ++i)
 
 #ifdef USE_STD_PARALLEL_ALG
 static inline const bool useStdParallel__ = true;
diff --git a/src/phasicFlow/repository/Time/Time.cpp b/src/phasicFlow/repository/Time/Time.cpp
index ebddd80c..4f80061b 100644
--- a/src/phasicFlow/repository/Time/Time.cpp
+++ b/src/phasicFlow/repository/Time/Time.cpp
@@ -85,7 +85,7 @@ bool pFlow::Time::write
 {
 	if(outputToFile())
 	{
-		Report(0)<<"\nWriting to file at time: "<< cyanText(currentTimeWord())<<endReport;
+		REPORT(0)<<"\nWriting to file at time: "<< cyanText(currentTimeWord())<<endREPORT;
 		return repository::write(verbose);
 	}
 	return true;
diff --git a/src/phasicFlow/repository/Time/timeControl.hpp b/src/phasicFlow/repository/Time/timeControl.hpp
index 4a7eb94f..9ad2916e 100644
--- a/src/phasicFlow/repository/Time/timeControl.hpp
+++ b/src/phasicFlow/repository/Time/timeControl.hpp
@@ -158,7 +158,7 @@ public:
 		currentTime_ += dt_;
 		if(screenReport())
 		{
-			Report(0)<<"Time (s): "<<cyanText( currentTimeWord() )<<endReport;
+			REPORT(0)<<"Time (s): "<<cyanText( currentTimeWord() )<<endREPORT;
 		}
 		// switch outputToFile_ on/off
 		checkForOutputToFile();
diff --git a/src/phasicFlow/repository/repository/repository.cpp b/src/phasicFlow/repository/repository/repository.cpp
index 79944974..520afa38 100644
--- a/src/phasicFlow/repository/repository/repository.cpp
+++ b/src/phasicFlow/repository/repository/repository.cpp
@@ -245,7 +245,7 @@ bool pFlow::repository::write
 	{
 		if(verbose)
 		{
-			Report(1)<< "Writing to " << obj.second.path()<<endReport;
+			REPORT(1)<< "Writing to " << obj.second.path()<<endREPORT;
 		}
 		
 		if(!obj.second.write())
diff --git a/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp b/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp
index d3bda996..32e63e96 100644
--- a/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp
+++ b/src/phasicFlow/setFieldList/setFieldEntryTemplates.cpp
@@ -50,8 +50,8 @@ void* pFlow::setFieldEntry::setPointFieldDefaultValueNew
 	Type defValue = entry_.secondPartVal<Type>();
 	
 	if(verbose)
-		Report(2)<<"Creating pointField " << greenText(fieldName())<< " with default value " << cyanText(defValue)<<
-	             " in repository "<< owner.name() <<endReport;
+		REPORT(2)<<"Creating pointField " << greenText(fieldName())<< " with default value " << cyanText(defValue)<<
+	             " in repository "<< owner.name() <<endREPORT;
 
 	
 	auto& field = 
@@ -85,8 +85,8 @@ void* pFlow::setFieldEntry::setPointFieldDefaultValueStdNew
 	Type defValue = entry_.secondPartVal<Type>();
 	
 	if(verbose)
-		Report(2)<<"Creating pointField " << greenText(fieldName())<< " with default value " << cyanText(defValue)<<
-	             " in repository "<< owner.name() <<endReport;
+		REPORT(2)<<"Creating pointField " << greenText(fieldName())<< " with default value " << cyanText(defValue)<<
+	             " in repository "<< owner.name() <<endREPORT;
 
 	// by default we perform operations on host
 	auto& field = 
@@ -129,8 +129,8 @@ void* pFlow::setFieldEntry::setPointFieldSelected
 	Type value = entry_.secondPartVal<Type>();
 
 	if(verbose)
-		Report(2)<< "Setting selected points of " << greenText(fName)
-		 		 << " to value " << cyanText(value) <<endReport;
+		REPORT(2)<< "Setting selected points of " << greenText(fName)
+		 		 << " to value " << cyanText(value) <<endREPORT;
 	
 	
 	auto fieldTypeName = owner.lookupObjectTypeName(fName);
@@ -196,8 +196,8 @@ void* pFlow::setFieldEntry::setPointFieldSelectedStd
 	Type value = entry_.secondPartVal<Type>();
 
 	if(verbose)
-		Report(2)<< "Setting selected points of " << greenText(fName)
-		 		 << " to value " << cyanText(value) <<endReport;
+		REPORT(2)<< "Setting selected points of " << greenText(fName)
+		 		 << " to value " << cyanText(value) <<endREPORT;
 	
 	
 	auto fieldTypeName = owner.lookupObjectTypeName(fName);
diff --git a/src/phasicFlow/streams/streams.hpp b/src/phasicFlow/streams/streams.hpp
index 29307da4..de60dc6f 100755
--- a/src/phasicFlow/streams/streams.hpp
+++ b/src/phasicFlow/streams/streams.hpp
@@ -34,16 +34,17 @@ namespace pFlow
 #define cyanText(text) 		cyanColor<<text<<defaultColor
 #define boldText(text) 		boldChar<<text<<defaultColor
 
-#define Info pFlow::output<<boldChar<<magentaColor<<"> Info: "<<defaultColor<<magentaColor
-#define endInfo defaultColor<<pFlow::nl
+#define INFORMATION pFlow::output<<boldChar<<magentaColor<<"> INFO: "<<defaultColor<<magentaColor
+#define endINFO defaultColor<<pFlow::nl
 
-#define Report(n) pFlow::output.space(2*n)
-#define endReport pFlow::nl
+#define REPORT(n) pFlow::output.space(2*n)
+#define endREPORT pFlow::nl
 
-#define Warning pFlow::output<<boldChar<<yellowColor<<"> Warning\n"<<defaultColor<<yellowColor<<"  " 
-#define endWarning defaultColor<<pFlow::nl
 
-#define Err pFlow::output<<boldChar<<redColor<<"> Error\n"<<defaultColor<<redColor<<"  " 
-#define endErr defaultColor<<pFlow::nl
+#define yWARNING pFlow::output<<boldChar<<yellowColor<<"> WARNING\n"<<defaultColor<<yellowColor<<"  " 
+#define endyWARNING defaultColor<<pFlow::nl
+
+#define ERR pFlow::output<<boldChar<<redColor<<"> ERROR\n"<<defaultColor<<redColor<<"  " 
+#define endERR defaultColor<<pFlow::nl
 
 #endif
diff --git a/src/phasicFlow/structuredData/box/box.hpp b/src/phasicFlow/structuredData/box/box.hpp
index 04eeb2df..428e0534 100644
--- a/src/phasicFlow/structuredData/box/box.hpp
+++ b/src/phasicFlow/structuredData/box/box.hpp
@@ -116,6 +116,11 @@ iIstream& operator >>(iIstream& is, box& b);
 FUNCTION_H
 iOstream& operator << (iOstream& os, const box& b);
 
+INLINE_FUNCTION_HD
+box extendBox(const box& b, const realx3& dl)
+{
+	return box(b.minPoint()-dl , b.maxPoint()+dl);
+}
 
 }
 
diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
index 6d1bf57c..6a21e808 100644
--- a/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.cpp
@@ -100,7 +100,7 @@ pFlow::uniquePtr<pFlow::int32IndexContainer>
 
 	newPoints.clear();
 	int32 numAdded = 0;
-	forAll(i, pointFlag_)
+	ForAll(i, pointFlag_)
 	{
 		if(!isActive(i))
 		{
@@ -469,7 +469,7 @@ pFlow::pointStructure::newPointsIndices(
 	// otherwise scan the points from first to the end to find empty spaces 
 	newPoints.clear();
 	int32 numAdded = 0;
-	forAll(i, pointFlag_)
+	ForAll(i, pointFlag_)
 	{
 		if(!isActive(i))
 		{
diff --git a/src/phasicFlow/structuredData/pointStructure/selectors/selectBox/selectBox.cpp b/src/phasicFlow/structuredData/pointStructure/selectors/selectBox/selectBox.cpp
index 6dd22236..89e89cfc 100644
--- a/src/phasicFlow/structuredData/pointStructure/selectors/selectBox/selectBox.cpp
+++ b/src/phasicFlow/structuredData/pointStructure/selectors/selectBox/selectBox.cpp
@@ -33,7 +33,7 @@ void pFlow::selectBox::selectAllPointsInBox()
 	selectedPoints_.clear();
 	auto pPos = pStruct().pointPosition().hostVector();
 	
-	forAll(i , pPos)
+	ForAll(i , pPos)
 	{
 		if( box_.isInside( pPos[i] )) selectedPoints_.push_back(i);
 	}
diff --git a/src/phasicFlow/structuredData/trisurfaceStructure/multiTriSurface.cpp b/src/phasicFlow/structuredData/trisurfaceStructure/multiTriSurface.cpp
index 56abdb60..f06d9ed9 100644
--- a/src/phasicFlow/structuredData/trisurfaceStructure/multiTriSurface.cpp
+++ b/src/phasicFlow/structuredData/trisurfaceStructure/multiTriSurface.cpp
@@ -58,7 +58,7 @@ void pFlow::multiTriSurface::calculateVars()
 	
 	pointsStartPos_.push_back(0);
 
-	forAll( i, surfaceNumPoints_)
+	ForAll( i, surfaceNumPoints_)
 	{
 		if(i==0)continue;
 
@@ -72,7 +72,7 @@ void pFlow::multiTriSurface::calculateVars()
 	verticesStartPos_.clear();
 
 	verticesStartPos_.push_back(0);
-	forAll(i, surfaceNumVertices_)
+	ForAll(i, surfaceNumVertices_)
 	{
 		if(i==0)continue;
 		verticesStartPos_.push_back(verticesStartPos_[i-1] + surfaceNumVertices_[i-1]);
diff --git a/src/setHelpers/initialize.hpp b/src/setHelpers/initialize.hpp
index 0b970dfe..8ee05c99 100644
--- a/src/setHelpers/initialize.hpp
+++ b/src/setHelpers/initialize.hpp
@@ -22,9 +22,9 @@ Licence:
 #define __initialize_hpp__
 
 // initilized and finalize should be placed in onc scope 
-Report(0)<<"Initializing host/device execution spaces . . . \n";
-Report(1)<<"Host execution space is "<< greenText(pFlow::DefaultHostExecutionSpace::name())<<endReport;
-Report(1)<<"Device execution space is "<<greenText(pFlow::DefaultExecutionSpace::name())<<endReport;
+REPORT(0)<<"Initializing host/device execution spaces . . . \n";
+REPORT(1)<<"Host execution space is "<< greenText(pFlow::DefaultHostExecutionSpace::name())<<endREPORT;
+REPORT(1)<<"Device execution space is "<<greenText(pFlow::DefaultExecutionSpace::name())<<endREPORT;
 
 Kokkos::initialize( argc, argv );
 {
diff --git a/src/setHelpers/initialize_Control.hpp b/src/setHelpers/initialize_Control.hpp
index a231a566..b488a6b7 100644
--- a/src/setHelpers/initialize_Control.hpp
+++ b/src/setHelpers/initialize_Control.hpp
@@ -24,7 +24,7 @@ Licence:
 // initilized and finalize should be placed in onc scope 
 #include "initialize.hpp"
 
-Report(0)<<"\nCreating Control repository . . ."<<endReport;
+REPORT(0)<<"\nCreating Control repository . . ."<<endREPORT;
 pFlow::uniquePtr<pFlow::systemControl> ControlPtr = nullptr;
 
 if(isCoupling)
diff --git a/src/setHelpers/setProperty.hpp b/src/setHelpers/setProperty.hpp
index 5c85adda..1d09d6aa 100644
--- a/src/setHelpers/setProperty.hpp
+++ b/src/setHelpers/setProperty.hpp
@@ -22,7 +22,7 @@ Licence:
 #ifndef __setProperty_hpp__
 #define __setProperty_hpp__
 
-Report(0)<<"\nReading proprties . . . "<<endReport;
+REPORT(0)<<"\nReading proprties . . . "<<endREPORT;
 
 auto proprties = pFlow::property(Control.caseSetup().path()+pFlow::propertyFile__);
 
diff --git a/src/setHelpers/setSurfaceGeometry.hpp b/src/setHelpers/setSurfaceGeometry.hpp
index 53a646ea..db6a0143 100644
--- a/src/setHelpers/setSurfaceGeometry.hpp
+++ b/src/setHelpers/setSurfaceGeometry.hpp
@@ -21,7 +21,7 @@ Licence:
 #ifndef __setSurfaceGeometry_hpp__
 #define __setSurfaceGeometry_hpp__
 
-Report(0)<< "\nCreating surface geometry . . . "<<endReport;
+REPORT(0)<< "\nCreating surface geometry . . . "<<endREPORT;
 auto surfGeometryPtr = pFlow::geometry::create(Control, proprties);
 auto& surfGeometry = surfGeometryPtr();
 
diff --git a/utilities/checkPhasicFlow/checkPhasicFlow.cpp b/utilities/checkPhasicFlow/checkPhasicFlow.cpp
index 9af45a33..3ffe4168 100755
--- a/utilities/checkPhasicFlow/checkPhasicFlow.cpp
+++ b/utilities/checkPhasicFlow/checkPhasicFlow.cpp
@@ -40,8 +40,8 @@ if(!cmds.parse(argc, argv)) return 0;
 #include "initialize.hpp"
 
 	output<<endl;
-	Report(1)<< "You are using "<<yellowText(cmds.productNameCopyright())<<endReport;
-  Report(1)<< yellowText(pFlow::floatingPointDescription())<<endReport;
+	REPORT(1)<< "You are using "<<yellowText(cmds.productNameCopyright())<<endREPORT;
+  REPORT(1)<< yellowText(pFlow::floatingPointDescription())<<endREPORT;
 	
 
 // this should be palced in each main 
diff --git a/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp b/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
index ef06a505..72559853 100755
--- a/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
+++ b/utilities/geometryPhasicFlow/geometryPhasicFlow.cpp
@@ -59,7 +59,7 @@ int main( int argc, char* argv[] )
 
 	#include "setProperty.hpp"
 
-	Report(0)<<"\nReading "<<"createGeometryDict"<<" . . ."<<endReport;
+	REPORT(0)<<"\nReading "<<"createGeometryDict"<<" . . ."<<endREPORT;
 	auto objDict = IOobject::make<dictionary>
 	(
 		objectFile
@@ -87,26 +87,26 @@ int main( int argc, char* argv[] )
 	
 	for(auto& name:wallsDictName)
 	{
-		Report(1)<<"Creating wall "<<greenText(name)<<" from file dictionary . . . "<<endReport;
+		REPORT(1)<<"Creating wall "<<greenText(name)<<" from file dictionary . . . "<<endREPORT;
 		auto wallPtr = Wall::create( surfacesDict.subDict(name));
 		auto& wall = wallPtr();
-		Report(1)<<"wall type is "<<greenText(wall.typeName())<<'\n'<<endReport;		
+		REPORT(1)<<"wall type is "<<greenText(wall.typeName())<<'\n'<<endREPORT;		
 
 		surface.addTriSurface(wall.name(), wall.triangles());
 		materials.push_back(wall.materialName());
 		motion.push_back(wall.motionName());
 	}
 
-	Report(1)<<"Selected wall materials are "<<cyanText(materials)<<'\n'<<endReport;
+	REPORT(1)<<"Selected wall materials are "<<cyanText(materials)<<'\n'<<endREPORT;
 		
-	Report(0)<< "\nCreating geometry . . ."<<endReport;
+	REPORT(0)<< "\nCreating geometry . . ."<<endREPORT;
 	auto geomPtr = geometry::create(Control, proprties, geometryDict, surface, motion, materials);
-	Report(1)<< "geometry type is "<< greenText(geomPtr().typeName())<<endReport;
+	REPORT(1)<< "geometry type is "<< greenText(geomPtr().typeName())<<endREPORT;
 	
-	Report(1)<< "Writing geometry to folder "<< geomPtr().path()<<endl;
+	REPORT(1)<< "Writing geometry to folder "<< geomPtr().path()<<endREPORT;
 	geomPtr().write();
 
-	Report(0)<< greenText("\nFinished successfully.\n");
+	REPORT(0)<< greenText("\nFinished successfully.\n");
 
 // this should be palced in each main 
 #include "finalize.hpp"
diff --git a/utilities/pFlowToVTK/geometric.hpp b/utilities/pFlowToVTK/geometric.hpp
index 2ace57c8..3c1a4313 100755
--- a/utilities/pFlowToVTK/geometric.hpp
+++ b/utilities/pFlowToVTK/geometric.hpp
@@ -45,7 +45,7 @@ bool geomObjectToVTK(IOfileHeader& header, real time, fileSystem destPath, word
 
 	if(!vtk) return false;
 
-	Report(1)<<"Converting geometry to vtk."<<endReport;
+	REPORT(1)<<"Converting geometry to vtk."<<endREPORT;
 
 	if(!dataToVTK(vtk, data) )
 	{
diff --git a/utilities/pFlowToVTK/pointFieldToVTK.hpp b/utilities/pFlowToVTK/pointFieldToVTK.hpp
index aebbf92f..ca71b028 100755
--- a/utilities/pFlowToVTK/pointFieldToVTK.hpp
+++ b/utilities/pFlowToVTK/pointFieldToVTK.hpp
@@ -105,7 +105,7 @@ bool convertIntTypesPointField(
 
 	int64* data = Field.hostVectorAll().data();
 	
-	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk.\n";
+	REPORT(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk.\n";
 
 	return addInt64PointField(
 		os,
@@ -136,7 +136,7 @@ bool convertRealTypePointField(
 
 	real* data = Field.hostVectorAll().data();
 	
-	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+	REPORT(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endREPORT;
 
 	return addRealPointField(
 		os,
@@ -166,7 +166,7 @@ bool convertRealx3TypePointField(
 
 	realx3* data = Field.hostVectorAll().data();
 	
-	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+	REPORT(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endREPORT;
 
 	return addRealx3PointField(
 		os,
@@ -324,8 +324,8 @@ bool convertTimeFolderPointFields(
 	auto posVec = std::as_const(pStruct).pointPosition().hostVectorAll();
 	auto* pos = posVec.data();
 
-	Report(1)<<"Writing pointStructure to vtk file with "<< yellowText(pStruct.numActive())
-					 <<" active particles"<<endReport;
+	REPORT(1)<<"Writing pointStructure to vtk file with "<< yellowText(pStruct.numActive())
+					 <<" active particles"<<endREPORT;
 	addUndstrcuturedGridField(
 		vtk(),
 		pStruct.numActive(),
@@ -395,8 +395,8 @@ bool convertTimeFolderPointFieldsSelected(
 	auto posVec = std::as_const(pStruct).pointPosition().hostVectorAll();
 	auto* pos = posVec.data();
 
-	Report(1)<<"Writing pointStructure to vtk file with "<< yellowText(pStruct.numActive())
-					 <<" active particles"<<endReport;
+	REPORT(1)<<"Writing pointStructure to vtk file with "<< yellowText(pStruct.numActive())
+					 <<" active particles"<<endREPORT;
 	addUndstrcuturedGridField(
 		vtk(),
 		pStruct.numActive(),
@@ -433,7 +433,7 @@ bool convertTimeFolderPointFieldsSelected(
 			}
 			else
 			{
-				Report(1)<<"Could not find "<<yellowText(fieldAddress) <<" skipping . . ."<<endReport;
+				REPORT(1)<<"Could not find "<<yellowText(fieldAddress) <<" skipping . . ."<<endREPORT;
 			}
 		}
 	}
diff --git a/utilities/pFlowToVTK/triSurfaceFieldToVTK.hpp b/utilities/pFlowToVTK/triSurfaceFieldToVTK.hpp
index 231bc2a7..d749a637 100644
--- a/utilities/pFlowToVTK/triSurfaceFieldToVTK.hpp
+++ b/utilities/pFlowToVTK/triSurfaceFieldToVTK.hpp
@@ -168,7 +168,7 @@ bool convertRealx3TypetriSurfaceField(
 
 	realx3* data = Field.hostVectorAll().data();
 	
-	Report(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endReport;
+	REPORT(2)<<"writing "<< greenColor <<header.objectName()<<defaultColor<<" field to vtk."<<endREPORT;
 
 	return addRealx3TriSurfaceField(
 		os,
@@ -209,7 +209,7 @@ bool convertTimeFolderTriSurfaceFields(
 	auto& tSurface = triSurfaceObjPtr().getObject<multiTriSurface>();
 
 	// get a list of files in this timeFolder;
-	Report(1)<<"Wrting triSurface mesh/Geometry to vtk file."<<endReport;
+	REPORT(1)<<"Wrting triSurface mesh/Geometry to vtk file."<<endREPORT;
 	if(!triDataToVTK(vtk(), tSurface))
 	{
 		fatalErrorInFunction<<
diff --git a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
index fbf834f8..0f22e484 100755
--- a/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
+++ b/utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
@@ -73,8 +73,8 @@ int main( int argc, char* argv[] )
 	
 	if(setOnly && positionOnly)
 	{
-		Err<<
-		"Options --positionParticles-only and --setFields-only cannot be used simeltanuously. \n"<<endErr;
+		ERR<<
+		"Options --positionParticles-only and --setFields-only cannot be used simeltanuously. \n"<<endERR;
 		return 1;
 	}	
 
@@ -102,7 +102,7 @@ int main( int argc, char* argv[] )
 	{
 
 		// position particles based on the dict content 
-		Report(0)<< "Positioning points . . . \n"<<endReport;
+		REPORT(0)<< "Positioning points . . . \n"<<endREPORT;
 		auto pointPosition = positionParticles::create(cpDict.subDict("positionParticles"));
 
 		fileSystem pStructPath = Control.time().path()+pointStructureFile__;
@@ -124,15 +124,15 @@ int main( int argc, char* argv[] )
 
 		auto& pSruct = pStructObj().getObject<pointStructure>();
 
-		Report(1)<< "Created pStruct with "<< pSruct.size() << " points and capacity "<<
-		pSruct.capacity()<<" . . ."<< endReport;
+		REPORT(1)<< "Created pStruct with "<< pSruct.size() << " points and capacity "<<
+		pSruct.capacity()<<" . . ."<< endREPORT;
 
-		Report(1)<< "Writing pStruct to " << pStructObj().path() << endReport<<endl<<endl;
+		REPORT(1)<< "Writing pStruct to " << pStructObj().path() << endREPORT<<endl<<endl;
 
 		if( !pStructObj().write())
 		{
 			fatalErrorInFunction<<
-			"Error in writing to file. \n ";
+			"ERRor in writing to file. \n ";
 			return 1; 
 		}
 	}else
@@ -164,7 +164,7 @@ int main( int argc, char* argv[] )
 		{
 			if( !sfEntry.setPointFieldDefaultValueNewAll(Control.time(), pStruct, true))
 			{
-				Err<< "\n error occured in setting default value fields.\n"<<endErr;
+				ERR<< "\n error occured in setting default value fields.\n"<<endERR;
 				return 1;
 			}
 		}
@@ -177,13 +177,13 @@ int main( int argc, char* argv[] )
 
 		for(auto name: selNames)
 		{
-			Report(1)<< "Applying selector " << greenText(name) <<endReport;
+			REPORT(1)<< "Applying selector " << greenText(name) <<endREPORT;
 			
 			if(
 				!applySelector(Control, pStruct, selectorsDict.subDict(name))
 			 )
 			{
-				Err<<"\n error occured in setting selector. \n"<<endErr;
+				ERR<<"\n error occured in setting selector. \n"<<endERR;
 				return 1;
 			}
 			output<<endl;
@@ -191,7 +191,7 @@ int main( int argc, char* argv[] )
 	}
 
 	Control.time().write(true);
-	Report(0)<< greenText("\nFinished successfully.\n")<<endReport;
+	REPORT(0)<< greenText("\nFinished successfully.\n")<<endREPORT;
 
 
 // this should be palced in each main 
diff --git a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
index 632851ec..0d3b8409 100755
--- a/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
+++ b/utilities/particlesPhasicFlow/positionParticles/positionParticles.cpp
@@ -54,7 +54,7 @@ pFlow::realx3Vector pFlow::positionParticles::sortByMortonCode(realx3Vector& pos
 			ind++});
 	}
 
-	Info<<"Performing morton sorting."<<endInfo;
+	INFORMATION<<"Performing morton sorting."<<endINFO;
 	std::sort(
 		indMor.begin(),
 		indMor.end(),
diff --git a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
index 63828fbc..616dbbca 100755
--- a/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
+++ b/utilities/particlesPhasicFlow/positionRandom/positionRandom.cpp
@@ -74,9 +74,9 @@ bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
 
 	ContainerType pairs(3*startNum);
 		
-	Report(1)<< "Positioning "<<
+	REPORT(1)<< "Positioning "<<
 		greenText("(Pass #"<< pass+1<<")")<< 
-		": started with "<< startNum <<" points."<<endReport;
+		": started with "<< startNum <<" points."<<endREPORT;
 	
 	fillPoints(startNum, positionHD, flagHD);
 		
@@ -86,13 +86,13 @@ bool pFlow::positionRandom::positionOnePass(int32 pass, int32 startNum)
 	int32 numCollisions = findCollisions(pairs, flagHD);
 	
 	
-	Report(2)<< "Positioned " << cyanText(startNum - numCollisions) <<
-	" without collision \n"<<endReport; 
+	REPORT(2)<< "Positioned " << cyanText(startNum - numCollisions) <<
+	" without collision \n"<<endREPORT; 
 
 	if(startNum-numCollisions >= numPoints_ )
 	{
 		
-		Report(1)<<"Selected "<< cyanText(numPoints_)<< " for the final field.\n"<<endReport;
+		REPORT(1)<<"Selected "<< cyanText(numPoints_)<< " for the final field.\n"<<endREPORT;
 
 		positionHD.syncViews();
 		position_.clear();
diff --git a/utilities/postprocessPhasicFlow/includeMask.cpp b/utilities/postprocessPhasicFlow/includeMask.cpp
index 2da4211b..0d5fe75b 100644
--- a/utilities/postprocessPhasicFlow/includeMask.cpp
+++ b/utilities/postprocessPhasicFlow/includeMask.cpp
@@ -85,7 +85,7 @@ pFlow::uniquePtr<pFlow::includeMask> pFlow::includeMask::create(
 		auto objPtr = 
 			dictionaryvCtorSelector_[method]
 			(dict, opType, timeFolder);
-		Report(2)<< dict.name()<< " with model "<<greenText(method)<<" is processed."<<endReport;
+		REPORT(2)<< dict.name()<< " with model "<<greenText(method)<<" is processed."<<endREPORT;
 		return objPtr;
 	}
 	else
diff --git a/utilities/postprocessPhasicFlow/postprocess.cpp b/utilities/postprocessPhasicFlow/postprocess.cpp
index 6e0fc199..8bfcdfb6 100644
--- a/utilities/postprocessPhasicFlow/postprocess.cpp
+++ b/utilities/postprocessPhasicFlow/postprocess.cpp
@@ -29,20 +29,20 @@ pFlow::postprocess::postprocess(systemControl& control)
 	control_(control),
 	dict_(postprocessFile__, control_.settings().path()+postprocessFile__)
 {
-	Report(1)<<"Reading numberBased dictionary ..."<<endReport;
+	REPORT(1)<<"Reading numberBased dictionary ..."<<endREPORT;
 	auto nbDict = dict_.subDictOrCreate("numberBased");
 
 	numberBasedDictNames_ = dict_.subDictOrCreate("numberBased").dictionaryKeywords();
 	if(!numberBasedDictNames_.empty())
 	{
-		Report(1)<< "numberBased dictionary contains " << yellowText(numberBasedDictNames_)<<endReport<<endl;	
+		REPORT(1)<< "numberBased dictionary contains " << yellowText(numberBasedDictNames_)<<endREPORT<<endl;	
 	}
 	
 
 	weightBasedDictNames_ = dict_.subDictOrCreate("weightBased").dictionaryKeywords();
 	if(!weightBasedDictNames_.empty())
 	{
-		Report(1)<< "numberBased dictionary contains " << yellowText(weightBasedDictNames_)<<endReport<<endl;	
+		REPORT(1)<< "numberBased dictionary contains " << yellowText(weightBasedDictNames_)<<endREPORT<<endl;	
 	}
 	
 
@@ -55,7 +55,7 @@ bool pFlow::postprocess::processTimeFolder(real time, const word& tName, const f
 	
 	time_ = time;
 
-	Report(0)<<"Working on time folder "<< cyanText(time)<<endReport; 
+	REPORT(0)<<"Working on time folder "<< cyanText(time)<<endREPORT; 
 	timeFolderReposiory_ = 
 		makeUnique<repository>
 		(
@@ -64,7 +64,7 @@ bool pFlow::postprocess::processTimeFolder(real time, const word& tName, const f
 			&control_
 		);
 	
-	Report(1)<<"Reading pointStructure"<<endReport;
+	REPORT(1)<<"Reading pointStructure"<<endREPORT;
 	timeFolderReposiory().emplaceObject<pointStructure>(
 		objectFile
 		(
@@ -75,7 +75,7 @@ bool pFlow::postprocess::processTimeFolder(real time, const word& tName, const f
 		));	
 
 	
-	Report(1)<<"Creating mesh and point to cell mapper"<<endReport;
+	REPORT(1)<<"Creating mesh and point to cell mapper"<<endREPORT;
 	pointToCell_ = makeUnique<pointRectCell>(
 		dict_.subDict("rectMesh"),
 		timeFolderReposiory().lookupObject<pointStructure>(pointStructureFile__),
@@ -125,11 +125,11 @@ bool pFlow::postprocess::writeToVTK(fileSystem destPath, word bName)const
 
 	if(!vtk) return false;
 
-	Report(1)<<"Writing processed fields to vtk file..."<<endReport;
+	REPORT(1)<<"Writing processed fields to vtk file..."<<endREPORT;
 	// mesh
 	pointToCell_->mesh().writeToVtk(vtk());
 
-	forAll( i, processedFields_)
+	ForAll( i, processedFields_)
 	{
 		
 		if( !processedFields_[i].writeToVTK(vtk()))
diff --git a/utilities/postprocessPhasicFlow/processField.cpp b/utilities/postprocessPhasicFlow/processField.cpp
index e53f2f2f..2521d535 100644
--- a/utilities/postprocessPhasicFlow/processField.cpp
+++ b/utilities/postprocessPhasicFlow/processField.cpp
@@ -115,7 +115,7 @@ pFlow::processField::create(
 		auto objPtr = 
 			dictionaryvCtorSelector_[method]
 			(dict, pToCell, rep);
-		Report(2)<<"Processing/creating " << yellowText(dict.name())<< " with model "<<greenText(method)<<"."<<endReport;
+		REPORT(2)<<"Processing/creating " << yellowText(dict.name())<< " with model "<<greenText(method)<<"."<<endREPORT;
 		return objPtr;
 	}
 	else