Particle base class and sphereParticle class have been updated

src/Particles/CMakeLists.txt

@@ -1,10 +1,12 @@
 
 set(SourceFiles 
 dynamicPointStructure/dynamicPointStructure.cpp
+particles/baseShapeNames.cpp
+particles/shape.cpp
 particles/particles.cpp
 #particles/particleIdHandler.cpp
-#SphereParticles/sphereShape/sphereShape.cpp
+SphereParticles/sphereShape/sphereShape.cpp
-#SphereParticles/sphereParticles/sphereParticles.cpp
+SphereParticles/sphereParticles/sphereParticles.cpp
 #Insertion/shapeMixture/shapeMixture.cpp
 #Insertion/insertion/insertion.cpp
 #Insertion/Insertion/Insertions.cpp

src/Particles/SphereParticles/sphereParticles/sphereParticles.cpp

@@ -19,10 +19,12 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "sphereParticles.hpp"
-#include "setFieldList.hpp"
+#include "systemControl.hpp"
-#include "sphereParticlesKernels.hpp"
+#include "vocabs.hpp"
+//#include "setFieldList.hpp"
+//#include "sphereParticlesKernels.hpp"
 
-pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
+/*pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
 pFlow::sphereParticles::getFieldObjectList()const
 {
 	auto objListPtr = particles::getFieldObjectList();
@@ -78,10 +80,10 @@ bool pFlow::sphereParticles::initializeParticles()
 		static_cast<int32>(shapeName_.size()));
 	
 	return insertSphereParticles(shapeName_, indices, false);
-}
+}*/
 
 
-bool pFlow::sphereParticles::beforeIteration() 
+/*bool pFlow::sphereParticles::beforeIteration() 
 {
 	particles::beforeIteration();
 	
@@ -98,45 +100,19 @@ bool pFlow::sphereParticles::beforeIteration()
 	intPredictTimer_.end();
 
 	return true;
-}
+}*/
 
 
-bool pFlow::sphereParticles::iterate() 
-{
-
-	accelerationTimer_.start();
-	//INFO<<"before accelerationTimer_ "<<endINFO;
-		pFlow::sphereParticlesKernels::acceleration(
-			control().g(),
-			mass().deviceVectorAll(),
-			contactForce().deviceVectorAll(),
-			I().deviceVectorAll(),
-			contactTorque().deviceVectorAll(),
-			pStruct().activePointsMaskD(),
-			accelertion().deviceVectorAll(),
-			rAcceleration().deviceVectorAll()
-			);
-	accelerationTimer_.end();
-	
-	intCorrectTimer_.start();
-	
-		dynPointStruct_.correct(this->dt(), accelertion_);
-	
-		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
-	
-	intCorrectTimer_.end();
-	
-	return true;
-}
 
 
-bool pFlow::sphereParticles::afterIteration() 
+
+/*bool pFlow::sphereParticles::afterIteration() 
 {
 	return true;
-}
+}*/
 
 
-bool pFlow::sphereParticles::insertSphereParticles(
+/*bool pFlow::sphereParticles::insertSphereParticles(
 	const wordVector& names,
 	const int32IndexContainer& indices,
 	bool setId
@@ -219,6 +195,39 @@ bool pFlow::sphereParticles::insertSphereParticles(
 
 	return true;
 
+}*/
+
+bool pFlow::sphereParticles::initInertia()
+{
+	
+	using exeSpace = typename realPointField_D::execution_space;
+	using policy = Kokkos::RangePolicy<
+			exeSpace,
+			Kokkos::IndexType<uint32>>;
+	
+	auto aPointsMask = dynPointStruct().activePointsMaskDevice();
+	auto aRange = aPointsMask.activeRange();
+
+	auto field_I = I_.fieldDevice();
+	auto field_shapeIndex = shapeIndex().fieldDevice();
+	
+	const auto& shp = getShapes();
+	realVector_D I ("I", shp.Inertia());
+	auto d_I = I.deviceVector();
+	
+	Kokkos::parallel_for(
+		"particles::initInertia",
+		policy(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(aPointsMask(i))
+			{
+				uint32 index = field_shapeIndex[i];
+				field_I[i] = d_I[index];
+			}
+		});
+	Kokkos::fence();
+	return true;
 }
 
 
@@ -227,78 +236,67 @@ pFlow::sphereParticles::sphereParticles(
 	const property& prop
 )
 :
-	particles(
+	particles(control),
-		control,
+	spheres_
-		control.settingsDict().getValOrSet(
+	(
-			"integrationMethod",
+		shapeFile__,
-			word("AdamsBashforth3")
+		&control.caseSetup(),
-			)
+		prop
+	),
+	I_
+	(
+		objectFile
+		(
+			"I",
+			"",
+			objectFile::READ_NEVER,
+			objectFile::WRITE_NEVER
 		),
-	property_(prop),
+		dynPointStruct(),
-	shapes_(
+		static_cast<real>(0.0000000001)
-		control.caseSetup().emplaceObjectOrGet<sphereShape>(
+	),
-			objectFile(
+	rVelocity_
-				sphereShapeFile__,
+	(
-				"",
+		objectFile
-				objectFile::READ_ALWAYS,
+		(
-				objectFile::WRITE_NEVER
+			"rVelocity",
-				)
+			"",
-			)
+			objectFile::READ_IF_PRESENT,
-		),
+			objectFile::WRITE_ALWAYS
-	I_(
-		this->time().emplaceObject<realPointField_D>(
-			objectFile(
-				"I",
-				"",
-				objectFile::READ_NEVER,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			static_cast<real>(0.0000000001)
-			)
-		),
-	rVelocity_(
-		this->time().emplaceObject<realx3PointField_D>(
-			objectFile(
-				"rVelocity",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			zero3
-			)
-		),
-	rAcceleration_(
-		this->time().emplaceObject<realx3PointField_D>(
-			objectFile(
-				"rAcceleration",
-				"",
-				objectFile::READ_IF_PRESENT,
-				objectFile::WRITE_ALWAYS
-				),
-			pStruct(),
-			zero3
-			)
 		),
+		dynPointStruct(),
+		zero3
+	),
+	rAcceleration_
+	(
+		objectFile(
+			"rAcceleration",
+			"",
+			objectFile::READ_IF_PRESENT,
+			objectFile::WRITE_ALWAYS
+			),
+		dynPointStruct(),
+		zero3
+	),
 	accelerationTimer_(
 		"Acceleration", &this->timers() ),
 	intPredictTimer_(
 		"Integration-predict", &this->timers() ),
 	intCorrectTimer_(
 		"Integration-correct", &this->timers() )
-
 {
 
-	REPORT(1)<<"Creating integration method "<<greenText(this->integrationMethod())
+	auto intMethod = control.settingsDict().getVal<word>("integrationMethod");
-		  << " for rotational velocity."<<endREPORT;
+	REPORT(1)<<"Creating integration method "<<Green_Text(intMethod)
+		  << " for rotational velocity."<<END_REPORT;
 		  
-	rVelIntegration_ = 
+	rVelIntegration_ = integration::create
-		integration::create(
+		(
 			"rVelocity",
-			this->time().integration(),
+			dynPointStruct(),
-			this->pStruct(),
+			intMethod,
-			this->integrationMethod());
+			rVelocity_.field()
+		);
 
 	if( !rVelIntegration_ )
 	{
@@ -307,7 +305,8 @@ pFlow::sphereParticles::sphereParticles(
 		fatalExit;
 	}
 
-	if(rVelIntegration_->needSetInitialVals())
+	WARNING<<"setFields for rVelIntegration_"<<END_WARNING;
+	/*if(rVelIntegration_->needSetInitialVals())
 	{
 		
 		auto [minInd, maxInd] = pStruct().activeRange();
@@ -327,17 +326,17 @@ pFlow::sphereParticles::sphereParticles(
 		REPORT(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endREPORT;
 		rVelIntegration_->setInitialVals(indexHD, rvel);
 
-	}
+	}*/
 
 	
-	if(!initializeParticles())
+	if(!initInertia())
 	{
 		fatalExit;
 	}
 	
 } 
 
-bool pFlow::sphereParticles::update(const eventMessage& msg) 
+/*bool pFlow::sphereParticles::update(const eventMessage& msg) 
 {
 	
 	if(rVelIntegration_->needSetInitialVals())
@@ -362,9 +361,9 @@ bool pFlow::sphereParticles::update(const eventMessage& msg)
 	}	
 
 	return true;
-}
+}*/
 
-bool pFlow::sphereParticles::insertParticles
+/*bool pFlow::sphereParticles::insertParticles
 (
 	const realx3Vector& position,
  	const wordVector&  shapes,
@@ -419,4 +418,70 @@ bool pFlow::sphereParticles::insertParticles
 
 	return true;
 		
-}
+}*/
+
+bool pFlow::sphereParticles::beforeIteration()
+{
+	particles::beforeIteration();
+	intPredictTimer_.start();
+		dynPointStruct().predict(dt(), accelertion());
+		rVelIntegration_().predict(dt(),rVelocity_, rAcceleration_);
+	intPredictTimer_.end();
+	
+	WARNING<<"pFlow::sphereParticles::beforeIteration()"<<END_WARNING;
+	return true;
+}
+
+bool pFlow::sphereParticles::iterate() 
+{
+
+	particles::iterate();
+	accelerationTimer_.start();
+	WARNING<<"pFlow::sphereParticlesKernels::acceleration"<<END_WARNING;
+	//INFO<<"before accelerationTimer_ "<<endINFO;
+		/*pFlow::sphereParticlesKernels::acceleration(
+			control().g(),
+			mass().deviceVectorAll(),
+			contactForce().deviceVectorAll(),
+			I().deviceVectorAll(),
+			contactTorque().deviceVectorAll(),
+			pStruct().activePointsMaskD(),
+			accelertion().deviceVectorAll(),
+			rAcceleration().deviceVectorAll()
+			);*/
+	accelerationTimer_.end();
+	
+	intCorrectTimer_.start();
+	
+		if(!dynPointStruct().correct(dt(), accelertion()))
+		{
+			return false;
+		}
+		if(!rVelIntegration_().correct(
+			dt(), 
+			rVelocity_, 
+			rAcceleration_))
+		{
+			return false;
+		}
+	
+	intCorrectTimer_.end();
+	
+	return true;
+}
+
+bool pFlow::sphereParticles::afterIteration()
+{
+	particles::afterIteration();
+	return true;
+}
+
+pFlow::word pFlow::sphereParticles::shapeTypeName()const
+{
+	return "sphere";
+}
+
+const pFlow::shape &pFlow::sphereParticles::getShapes() const
+{
+    return spheres_;
+}

src/Particles/SphereParticles/sphereParticles/sphereParticles.hpp

@@ -44,22 +44,21 @@ class sphereParticles
 :
 	public particles 
 {
-protected:
+public:
+	using ShapeType = sphereShape;
+private:
 
-	/// reference to material properties
+		/// reference to shapes 
-	const property& 	property_;
+	ShapeType 			spheres_;
-
-	/// reference to shapes 
-	sphereShape& 		shapes_;
 
 	/// pointField of inertial of particles 
-	realPointField_D&     	I_;
+	realPointField_D     	I_;
 
 	/// pointField of rotational Velocity of particles on device 
-	realx3PointField_D& 	rVelocity_;
+	realx3PointField_D 	rVelocity_;
 
 	/// pointField of rotational acceleration of particles on device  
-	realx3PointField_D& 	rAcceleration_;	
+	realx3PointField_D 	rAcceleration_;	
 
 	/// rotational velocity integrator 
 	uniquePtr<integration>  rVelIntegration_ = nullptr;
@@ -73,9 +72,9 @@ protected:
 	/// timer for integration computations (correction step)
 	Timer 					intCorrectTimer_;
 
-	bool diameterMassInertiaPropId(const word& shName, real& diam, real& mass, real& I, int8& propIdx);
+	bool initInertia();
 
-	bool initializeParticles();
+	/*bool initializeParticles();
 
 	bool insertSphereParticles(
 		const wordVector& names, 
@@ -83,25 +82,17 @@ protected:
 		bool setId = true);
 
 	virtual uniquePtr<List<eventObserver*>> getFieldObjectList()const override;
+	*/
 
 public:
 
 	/// construct from systemControl and property 
-	sphereParticles(systemControl &control, const property& prop);
+	sphereParticles(
+		systemControl &control, 
+		const property& prop);
 
-	/// no copy constructor 
-	sphereParticles(const sphereParticles&) = delete;
 
-	/// move constructor 
+	~sphereParticles()override=default;
-	sphereParticles(sphereParticles&&) = default;
-
-	/// no copy-assignement
-	sphereParticles& operator=(const sphereParticles&) = delete;
-
-	/// move-assignement
-	sphereParticles& operator=(sphereParticles&&) = default;
-
-	virtual ~sphereParticles()=default;
 
 	/**
 	 * Insert new particles in position with specified shapes
@@ -112,17 +103,17 @@ public:
 	 * \param shape shape of new particles
 	 * \param setField initial value of the selected fields for new particles
 	 */
-	bool insertParticles
+	/*bool insertParticles
 	(
 		const realx3Vector& position,
 	 	const wordVector&  shapes,
 	 	const setFieldList& setField
-	) override ;
+	) override ;*/
 
 	/// const reference to shapes object 
-	const auto& shapes()const
+	const auto& spheres()const
 	{
-		return shapes_;
+		return spheres_;
 	}
 
 	/// const reference to inertia pointField
@@ -137,41 +128,29 @@ public:
 		return I_;
 	}
 	
-	const realx3Vector_D rVelocity()const 
+	const auto& rVelocity()const 
 	{	
 		return rVelocity_;
 	}
 
-	const realVector_D& boundingSphere()const override
+	auto& rVelocity()
-	{
+	{	
-		return this->diameter();
+		return rVelocity_;
 	}
 
-	word shapeTypeName() const override
+	bool hearChanges
+	(
+		real t,
+		real dt,
+		uint32 iter,
+		const message& msg, 
+    	const anyList& varList
+	) override
 	{
-		return "sphere";
+		notImplementedFunction;
-	}
+		return false;
-
-	void boundingSphereMinMax(
-		real& minDiam,
-		real& maxDiam )const override
-	{
-		shapes_.diameterMinMax(minDiam, maxDiam);
-	}
-
-	bool update(const eventMessage& msg) override;
-	
-	
-	realx3PointField_D& rAcceleration() override
-	{
-		return rAcceleration_;
 	}
 	
-	const realx3PointField_D& rAcceleration() const override
-	{
-		return rAcceleration_;
-	}
-
 	/// before iteration step 
 	bool beforeIteration() override;
 
@@ -180,8 +159,25 @@ public:
 
 	/// after iteration step
 	bool afterIteration() override;
-	
 
+	realx3PointField_D& rAcceleration() override
+	{
+		return rAcceleration_;
+	}
+
+	const realx3PointField_D& rAcceleration() const override
+	{
+		return rAcceleration_;
+	}
+	
+	const realPointField_D& boundingSphere()const override
+	{
+		return diameter();
+	}
+
+	word shapeTypeName()const override;
+	
+	const shape& getShapes()const override;
 	
 }; //sphereParticles
 

src/Particles/SphereParticles/sphereShape/sphereShape.cpp

@@ -19,186 +19,177 @@ Licence:
 -----------------------------------------------------------------------------*/
 
 #include "sphereShape.hpp"
-#include "dictionary.hpp"
-#include "vocabs.hpp"
-#include "streams.hpp"
 
 
-bool pFlow::sphereShape::insertNames
+bool pFlow::sphereShape::readDictionary()
-(
-	const wordVector& names
-)
 {
-	names_.clear();
-	uint32 i=0;
-	for(const auto& nm:names)
-	{
-		if(!names_.insertIf(nm, i))
-		{
-			fatalErrorInFunction<<
-			"  repeated name in the list of sphere names: " << names;
-			return false;
-		}
-		i++;
-	}
-	names_.rehash(0);
 	
-	numShapes_ = names_.size();
+	diameters_ = getVal<realVector>("diameters");
-	
-	return true;
-}
 
-
+	if(diameters_.size() != numShapes() )
-bool pFlow::sphereShape::readDictionary
-(
-	const dictionary& dict
-)
-{
-	diameters_ = dict.getVal<realVector>("diameters");
-	materials_ = dict.getVal<wordVector>("materials");
-	auto names = dict.getVal<wordVector>("names");
-
-	if(diameters_.size() != materials_.size() )
 	{
-		fatalErrorInFunction<<
+		fatalErrorInFunction<< 
-		"  number of elements in diameters and properties are not the same in "<< dict.globalName()<<endl;
+		" number of elements in diameters in "<< globalName()<<" is not consistent"<<endl;
-		return false;
-	}
-	else if(diameters_.size() != names.size() )
-	{
-		fatalErrorInFunction<<
-		"  number of elements in diameters and names are not the same in "<< dict.globalName()<<endl;
-		return false;
-	}
-
-	if( !insertNames(names) )
-	{
-		fatalErrorInFunction<<
-		"  error in reading dictionary "<< dict.globalName();
 		return false;
 	}
 
 	return true;
 }
 
-bool pFlow::sphereShape::writeDictionary
+bool pFlow::sphereShape::writeToDict(dictionary& dict)const 
-(
-	dictionary& dict
-)const
 {
+	
+	if(!shape::writeToDict(dict))return false;
 
-	if( !dict.add("diamters", diameters_) )
+	if( !dict.add("diameters", diameters_) )
 	{
-		fatalErrorInFunction<<  
+		fatalErrorInFunction<<
 		"  Error in writing diameters to dictionary "<< dict.globalName()<<endl;
 		return false;
 	}
 
-	if( !dict.add("properties", materials_) )
-	{
-		fatalErrorInFunction<<
-		"  Error in writing properties to dictionary "<< dict.globalName()<<endl;
-		return false;
-	}
-
-	size_t n = names_.size();
-	wordVector names(n);
-	names.clear();
-	word nm;
-	
-	for(label i=0; i<n; i++)
-	{
-		indexToName(i, nm);
-		names.push_back(nm);
-	}
-
-	if( !dict.add("names", names) )
-	{
-		fatalErrorInFunction<<
-		"  Error in writing names to dictionary "<< dict.globalName()<<endl;
-		return false;
-	}
-
 	return true;
 }
 
-
-
 pFlow::sphereShape::sphereShape
 (
-	const realVector& diameter,
+	const word& fileName,
-	const wordVector& property,
+	repository* owner,
-	const wordVector& name
+	const property& prop
 )
 :
-	diameters_(diameter),
+	shape(fileName, owner, prop)
-	materials_(property)
-{
-	if( !insertNames( name) )
-	{
-		fatalExit;
-	}
-}
-
-bool pFlow::sphereShape::shapeToDiameter
-(
-	wordVector& names,
-	realVector& diams
-)const
-{
-	diams.clear();
-	uint32 idx;
-	for(const auto& nm:names)
-	{
-		if(!nameToIndex(nm, idx))
-		{
-			fatalErrorInFunction<<
-			"  invalid shape name requested "<< nm <<endl;
-			return false;
-		}
-		diams.push_back(diameters_[idx]);
-	}
-
-	return true;
-}
-
-
-bool pFlow::sphereShape::read(iIstream& is)
 {
 
-	dictionary sphereDict(sphereShapeFile__, true);
-
-	if( !sphereDict.read(is) )
-	{
-		ioErrorInFile(is.name(), is.lineNumber()) <<
-		"  error in reading dictionray " << sphereShapeFile__ <<" from file. \n";
-		return false;
-	}
-
-	if( !readDictionary(sphereDict) )
-	{
-		return false;
-	}
-
-	return true;
 }
 
-bool pFlow::sphereShape::write(iOstream& os)const
+pFlow::real pFlow::sphereShape::maxBoundingSphere() const
 {
-
+    return max(diameters_);
-	dictionary sphereDict(sphereShapeFile__, true);
+}
-
+
-	if( !writeDictionary(sphereDict))
+pFlow::real pFlow::sphereShape::minBoundingSphere() const
-	{
+{
-		return false;
+    return min(diameters_);
-	}
+}
-
+
-	if( !sphereDict.write(os) )
+bool pFlow::sphereShape::boundingDiameter(uint32 index, real &bDiam) const
-	{
+{
-		ioErrorInFile( os.name(), os.lineNumber() )<<
+	if( indexValid(index))
-		"  error in writing dictionray to file. \n";
+	{
-		return false;
+		bDiam = diameters_[index];
-	}
+		return true;
-
+	}
-	return true;
+    return false;
+}
+
+pFlow::real pFlow::sphereShape::boundingDiameter(uint32 index) const
+{
+    if(indexValid(index))
+	{
+		return diameters_[index];
+	}
+	fatalErrorInFunction<<"Invalid index for diameter "<<
+	index<<endl;
+	fatalExit;
+	return 0.0;
+}
+
+pFlow::realVector pFlow::sphereShape::boundingDiameter() const
+{
+    return diameters_;
+}
+
+bool pFlow::sphereShape::mass(uint32 index, real &m) const
+{
+    if( indexValid(index) )
+	{
+		real d = diameters_[index];
+		real rho = indexToDensity(index);
+		m = Pi/6.0*pow(d,3)*rho;
+		return true;
+	}
+	return false;
+}
+
+pFlow::real pFlow::sphereShape::mass(uint32 index) const
+{
+    if(real m; mass(index, m))
+	{
+		return m;
+	}
+	fatalErrorInFunction<<"bad index for mass "<< index<<endl;
+	fatalExit;
+	return 0;
+}
+
+pFlow::realVector pFlow::sphereShape::mass() const
+{
+    return realVector ("mass", Pi/6*pow(diameters_,(real)3.0)*density());
+}
+
+pFlow::realVector pFlow::sphereShape::density()const 
+{
+	auto pids = shapePropertyIds();
+	realVector rho("rho", numShapes());
+	ForAll(i, pids)
+	{
+		rho[i] = properties().density(pids[i]);
+	}
+	return rho;
+}
+
+bool pFlow::sphereShape::Inertia(uint32 index, real &I) const
+{
+	if( indexValid(index) )
+	{
+		I = 0.4 * mass(index) * pow(diameters_[index]/2.0,2.0);
+		return true;
+	}
+    return false;
+}
+
+pFlow::real pFlow::sphereShape::Inertia(uint32 index) const
+{
+	if(real I; Inertia(index, I))
+	{
+		return I;
+	}
+    fatalExit;
+	return 0;
+}
+
+pFlow::realVector pFlow::sphereShape::Inertia() const
+{
+    return realVector("I", 0.4*mass()*pow(0.5*diameters_,(real)2.0));
+}
+
+bool pFlow::sphereShape::Inertia_xx(uint32 index, real &Ixx) const
+{
+    return Inertia(index,Ixx);
+}
+
+pFlow::real pFlow::sphereShape::Inertial_xx(uint32 index) const
+{
+    return Inertia(index);
+}
+
+bool pFlow::sphereShape::Inertia_yy(uint32 index, real &Iyy) const
+{
+    return Inertia(index,Iyy);
+}
+
+pFlow::real pFlow::sphereShape::Inertial_yy(uint32 index) const
+{
+    return Inertia(index);
+}
+
+bool pFlow::sphereShape::Inertia_zz(uint32 index, real &Izz) const
+{
+    return Inertia(index,Izz);
+}
+
+pFlow::real pFlow::sphereShape::Inertial_zz(uint32 index) const
+{
+    return Inertia(index);
 }

src/Particles/SphereParticles/sphereShape/sphereShape.hpp

@@ -21,154 +21,78 @@ Licence:
 #ifndef __sphereShape_hpp__
 #define __sphereShape_hpp__
 
-#include "Vectors.hpp"
+#include "shape.hpp"
-#include "hashMap.hpp"
 
 namespace pFlow
 {
 
-class dictionary;
+
 
 class sphereShape
+:
+	public shape
 {
-protected:
+private:
 	
 	// - diameter of spheres
 	realVector 			diameters_;
 	
-	// - property name of spheres
+	bool readDictionary();
-	wordVector 			materials_;
 
-	// - hashed list of spheres names
+protected:
-	wordHashMap<uint32> names_;
 
-	size_t 				numShapes_;
+	bool writeToDict(dictionary& dict)const override;
-
-
-	bool insertNames(const wordVector& names);
-
-	bool readDictionary(const dictionary& dict);
-
-	bool writeDictionary(dictionary& dict)const;
 
 public:
 
 	// - type info
-	TypeInfoNV("sphereShape");
+	TypeInfo("shape<sphere>");
 
+	sphereShape(
+		const word& fileName,
+		repository* owner,
+		const property& prop);
 
-		sphereShape(){}
+	
-
+	~sphereShape() override = default;
-		sphereShape(
-			const realVector& diameter,
-			const wordVector& property,
-			const wordVector& name
-			);
-
-		sphereShape(const sphereShape&) = default;
-
-		sphereShape(sphereShape&&) = default;
-
-		sphereShape& operator=(const sphereShape&) = default;
-
-		sphereShape& operator=(sphereShape&&) = default;
-
-		auto clone()const
-		{
-			return makeUnique<sphereShape>(*this);
-		}
-
-		sphereShape* clonePtr()const
-		{
-			return new sphereShape(*this);
-		}
-
-		~sphereShape() = default;
 
 	//// - Methods
-		const auto& names()const{
-			return names_;
-		}
 
-		const auto& diameters()const{
+	real maxBoundingSphere()const override;
-			return diameters_;
-		}
 
-		const auto& materials()const{
+    real minBoundingSphere()const override;
-			return materials_;
-		}
 
-		const auto diameter(label i)const{
+	bool boundingDiameter(uint32 index, real& bDiam)const override;
-			return diameters_[i];
-		}
 
-		const auto material(label i)const{
+	real boundingDiameter(uint32 index)const override;
-			return materials_[i];
-		}
 
+	realVector boundingDiameter()const override;
 
-		// name to index
+	bool mass(uint32 index, real& m)const override;
-		bool nameToIndex(const word& name, uint32& index)const
-		{
-			if(auto[iter, found] = names_.findIf(name); found )
-			{
-				index = iter->second;
-				return true;
-			}
-			else
-			{
-				index = 0;
-				return false;
-			}
-		}
 
-		uint32 nameToIndex(const word& name)const
+	real mass(uint32 index) const override;
-		{
-			return names_.at(name);
-		}
 
-		bool indexToName(uint32 i, word& name)const
+	realVector mass()const override;
-		{
-			for(auto& nm: names_)
-			{
-				if(nm.second == i)
-				{
-					name = nm.first;
-					return true;
-				}
-			}
-			name = "";
-			return false;
-		}
 
-		bool shapeToDiameter(wordVector& names, realVector& diams)const;
+    realVector density() const override;
 
-		void diameterMinMax(real& minD, real& maxD)const
+	bool Inertia(uint32 index, real& I)const override;
-		{
-			minD = min(diameters_);
-			maxD = max(diameters_);
-		}
 
-	//// - IO operatoin 
+	real Inertia(uint32 index)const override;
 
-		// - read from stream/file
+	realVector Inertia()const override;
-		bool read(iIstream& is);
+	 
+    bool Inertia_xx(uint32 index, real& Ixx)const override;
 
-		// - write to stream/file
+    real Inertial_xx(uint32 index)const override;
-		bool write(iOstream& os)const;
+	
+    bool Inertia_yy(uint32 index, real& Iyy)const override;
 
-		// - read from dictionary 
+    real Inertial_yy(uint32 index)const override;
-		bool read(const dictionary& dict)
-		{
-			return readDictionary(dict);
-		}
 
-		// - write to dictionary 
+    bool Inertia_zz(uint32 index, real& Izz)const override;
-		bool write(dictionary& dict)const
-		{
-			return writeDictionary(dict);
-		}
 
+    real Inertial_zz(uint32 index)const override;
 
 };
 

src/Particles/particles/baseShapeNames.cpp

@@ -0,0 +1,114 @@
+/*------------------------------- 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 "baseShapeNames.hpp"
+
+bool pFlow::baseShapeNames::createHashNames()
+{
+	hashNames_.clear();
+	uint32 i=0;
+	for(const auto& nm:shapeNames_)
+	{
+		if(!hashNames_.insertIf(nm, i))
+		{
+			fatalErrorInFunction<<
+			"  repeated name in the list of shape names: " << shapeNames_;
+			return false;
+		}
+		i++;
+	}
+	hashNames_.rehash(0);
+	
+	return true;
+}
+
+
+bool pFlow::baseShapeNames::readFromDictionary()
+{
+	
+	shapeNames_ = getVal<wordVector>("names");
+	numShapes_ = shapeNames_.size();
+
+	return true;
+}
+
+
+pFlow::baseShapeNames::baseShapeNames
+(
+    const word &fileName, 
+    repository *owner
+)
+:
+    fileDictionary
+    (
+        objectFile
+        (
+            fileName,
+            "",
+            objectFile::READ_ALWAYS,
+            objectFile::WRITE_ALWAYS
+        ),
+        owner
+    )
+{
+	
+	if( !readFromDictionary() )
+	{
+		fatalExit;
+	}
+
+	if( !createHashNames())
+	{
+		fatalExit;
+	}
+
+}
+
+bool pFlow::baseShapeNames::writeToDict(dictionary &dict)const
+{
+	
+	if( !dict.add("names", shapeNames_) )
+	{
+		fatalErrorInFunction<<
+		"  Error in writing names to dictionary "<< dict.globalName()<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+bool pFlow::baseShapeNames::write(iOstream &os) const
+{
+	dictionary newDict(fileDictionary::dictionary::name(), true);
+
+	if( !writeToDict(newDict) )
+	{
+		return false;
+	}
+
+	if( !newDict.write(os) )
+	{
+		fatalErrorInFunction<<"Error in writing dictionary "<< globalName()<<
+		"to stream "<< os.name()<<endl; 
+		return false;
+	}
+
+	return true;
+}

src/Particles/particles/baseShapeNames.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 __baseShapeNames_hpp__
+#define __baseShapeNames_hpp__
+
+#include "fileDictionary.hpp"
+#include "Vectors.hpp"
+#include "hashMap.hpp"
+
+
+namespace pFlow
+{
+
+class repository;
+
+class baseShapeNames
+:
+	public fileDictionary
+{
+private:
+
+	size_t 				numShapes_;
+
+	// - hashed list of spheres names
+	wordHashMap<uint32> hashNames_;
+
+	/// list of shape names 
+	wordVector 			shapeNames_;
+	
+	bool createHashNames();
+
+	bool readFromDictionary();	
+
+protected:
+
+	virtual
+	bool writeToDict(dictionary& dict)const;
+
+public:
+    
+    TypeInfo("baseShapeNames");
+
+	baseShapeNames(
+		const word& fileName,
+		repository* owner);
+
+    ~baseShapeNames() override=default;
+	
+
+    inline
+    const wordVector& shapeNames()const
+	{
+		return shapeNames_;
+	}
+
+	inline 
+	wordList shapeNameList()const
+	{
+		wordList wl;
+		wl.insert(wl.begin(), shapeNames_.begin(), shapeNames_.end());
+		return wl;
+	}
+
+    inline
+    size_t numShapes()const
+	{
+		return numShapes_;
+	}
+
+	// name to index
+	inline
+	bool shapeNameToIndex(const word& name, uint32& index)const
+	{
+		if(auto[iter, found] = hashNames_.findIf(name); found )
+		{
+			index = iter->second;
+			return true;
+		}
+		else
+		{
+			index = 0;
+			return false;
+		}
+	}
+
+	inline
+	bool indexToShapeName(uint32 i, word& name)const
+	{
+		if( i < numShapes_)
+		{
+			name = shapeNames_[i];
+			return true;
+		}
+		return false;
+	}
+
+	inline
+	bool indexValid(uint32 index)const
+	{
+		return index < numShapes_;
+	}
+
+	// - IO
+
+	bool write(iOstream& os)const override;
+
+};
+
+}
+
+#endif //__shapeNames_hpp__

src/Particles/particles/particles.cpp

@@ -21,17 +21,71 @@ Licence:
 
 #include "particles.hpp"
 
+bool pFlow::particles::initMassDiameter()const 
+{
+
+	auto [minIndex, maxIndex] = minMax(shapeIndex_.internal());
+
+	const auto& shp = getShapes();
+
+	if( !shp.indexValid(maxIndex) )
+	{
+		fatalErrorInFunction<< 
+		"the maximum value of shape index is "<< maxIndex << 
+		" which is not valid."<<endl;
+		return false;
+	}
+
+	realVector_D d("diameter", shp.boundingDiameter());
+	realVector_D mass("mass",shp.mass());
+	uint8Vector_D propId("propId", shp.shapePropertyIds());
+
+	
+	auto aPointsMask = dynPointStruct_.activePointsMaskDevice();
+	auto aRange = aPointsMask.activeRange();
+
+	using exeSpace = typename realPointField_D::execution_space;
+	using policy = Kokkos::RangePolicy<
+			exeSpace,
+			Kokkos::IndexType<uint32>>;
+
+	auto field_diameter = diameter_.fieldDevice();
+	auto field_mass = mass_.fieldDevice();
+	auto field_propId = propertyId_.fieldDevice();
+	auto field_shapeIndex = shapeIndex_.fieldDevice();
+
+	auto d_d = d.deviceVector();
+	auto d_mass = mass.deviceVector();
+	auto d_propId = propId.deviceVector();
+
+	Kokkos::parallel_for(
+		"particles::initMassDiameter",
+		policy(aRange.start(), aRange.end()),
+		LAMBDA_HD(uint32 i)
+		{
+			if(aPointsMask(i))
+			{
+				uint32 index = field_shapeIndex[i];
+				field_diameter[i] = d_d[index];
+				field_mass[i] = d_mass[index];
+				field_propId[i] = d_propId[index];
+			}
+		});
+	Kokkos::fence();
+
+    return true;
+}
 
 pFlow::particles::particles
 (
-	systemControl& control
+    systemControl &control
 )
-:
+: 
-	observer(),
+	observer(defaultMessage_),
 	demComponent("particles", control),
 	dynPointStruct_(control),
 	id_
-	(	
+	(
 		objectFile
 		(
 			"id",
@@ -48,127 +102,108 @@ pFlow::particles::particles
 		(
 			"propertyId",
 			"",
-			objectFile::READ_NEVER,
+			objectFile::READ_ALWAYS,
 			objectFile::WRITE_NEVER
 		),
 		dynPointStruct_,
 		static_cast<int8>(0)
 	),
-	diameter_
+	shapeIndex_
 	(
 		objectFile
 		(
+			"shapeIndex",
+			"",
+			objectFile::READ_ALWAYS,
+			objectFile::WRITE_ALWAYS
+		),
+		dynPointStruct_,
+		0
+	),
+	diameter_
+	(
+		objectFile(
 			"diameter",
 			"",
 			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
+			objectFile::WRITE_NEVER),
-		),
 		dynPointStruct_,
 		0.00000000001
 	),
 	mass_
 	(
-		objectFile
+		objectFile(
-		(
 			"mass",
 			"",
 			objectFile::READ_NEVER,
-			objectFile::WRITE_ALWAYS
+			objectFile::WRITE_NEVER),
-		),
 		dynPointStruct_,
 		0.0000000001
 	),
 	accelertion_
 	(
-		objectFile
+		objectFile(
-		(
 			"accelertion",
 			"",
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			objectFile::WRITE_ALWAYS),
-			),
 		dynPointStruct_,
 		zero3
 	),
-	contactForce_
+	contactForce_(
-	(
+		objectFile(
-		objectFile
-		(
 			"contactForce",
 			"",
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			objectFile::WRITE_ALWAYS),
-		),
 		dynPointStruct_,
-		zero3
+		zero3),
-	),
+	contactTorque_(
-	contactTorque_
+		objectFile(
-	(
-	
-		objectFile
-		(
 			"contactTorque",
 			"",
 			objectFile::READ_IF_PRESENT,
-			objectFile::WRITE_ALWAYS
+			objectFile::WRITE_ALWAYS),
-		),
 		dynPointStruct_,
-		zero3
+		zero3)
-	)
 {
 	
-	WARNING<<"Subscribe particles"<<END_WARNING;
+	this->addToSubscriber(dynPointStruct_);
-	//this->subscribe(pStruct());
+
+	if(!initMassDiameter())
+	{
+		fatalExit;
+	}
 
 }
 
 bool pFlow::particles::beforeIteration() 
 {
-	/*auto domain = this->control().domain();
-
-	auto numMarked = dynPointStruct_.markDeleteOutOfBox(domain);
 	
-	if(time_.sortTime())
+	dynPointStruct_.beforeIteration();
-	{
+	zeroForce();
-		real min_dx, max_dx;
+	zeroTorque();
-		boundingSphereMinMax(min_dx, max_dx);
-		Timer t;
-		t.start();
-		REPORT(0)<<"Performing morton sorting on particles ...."<<endREPORT;
-		if(!pStruct().mortonSortPoints(domain, max_dx))
-		{
-			fatalErrorInFunction<<"Morton sorting was not successful"<<endl;
-			return false;
-		}
-		t.end();
-		REPORT(1)<<"Active range is "<< pStruct().activeRange()<<endREPORT;
-		REPORT(1)<<"It took "<< yellowText(t.totalTime())<<" seconds."<<endREPORT;
-	}
-
-	this->zeroForce();
-	this->zeroTorque();*/
 
 	return true;
 }
 
-/*pFlow::uniquePtr<pFlow::List<pFlow::eventObserver*>> 
+bool pFlow::particles::iterate()
-pFlow::particles::getFieldObjectList()const
 {
-	auto objListPtr = makeUnique<pFlow::List<pFlow::eventObserver*>>();
+	return dynPointStruct_.iterate();
-	objListPtr().push_back(
+}
-		static_cast<eventObserver*>(&id_) );
-		
-	objListPtr().push_back(
-		static_cast<eventObserver*>(&propertyId_) );
 
-	objListPtr().push_back(
+bool pFlow::particles::afterIteration() 
-		static_cast<eventObserver*>(&diameter_) );
+{
+	return dynPointStruct_.afterIteration();
+}
 
-	objListPtr().push_back(
+void pFlow::particles::boundingSphereMinMax
-		static_cast<eventObserver*>(&mass_) );
+(
-	
+	real &minDiam, 
-	objListPtr().push_back(
+	real &maxDiam
-		static_cast<eventObserver*>(&shapeName_) );
+) const
-	
+{
-	return objListPtr;
+	auto& shp = getShapes();
-}*/
+	minDiam = shp.minBoundingSphere();
+	maxDiam = shp.maxBoundingSphere();
+}

src/Particles/particles/particles.hpp

@@ -24,6 +24,7 @@ Licence:
 
 #include "dynamicPointStructure.hpp"
 #include "demComponent.hpp"
+#include "shape.hpp"
 
 namespace pFlow
 {
@@ -39,43 +40,51 @@ private:
 	/// dynamic point structure for particles center mass
 	dynamicPointStructure 	dynPointStruct_;
 
-	// - name of shapes - this is managed by particles 
+	/// id of particles on host
-	//wordPointField&  		shapeName_;
+	uint32PointField_D		id_;
 
-	// id of particles on host
+	/// property id on device
-	int32PointField_D		id_;
+	uint8PointField_D 		propertyId_;
 
-	// property id on device
+	/// property id on device
-	int8PointField_D 		propertyId_;
+	uint32PointField_D 		shapeIndex_;
 
-	// diameter / boundig sphere size of particles on device
+	/// diameter / boundig sphere size of particles on device
 	realPointField_D		diameter_;
 
-	// mass of particles field 
+	/// mass of particles field 
 	realPointField_D 		mass_;
 
-	// - acceleration on device
+	/// acceleration on device
-	realx3PointField_D 	accelertion_;
+	realx3PointField_D 		accelertion_;
 
 	/// contact force field
 	realx3PointField_D     contactForce_;
 
 	/// contact torque field
-	realx3PointField_D     contactTorque_;	
+	realx3PointField_D     contactTorque_;
 
+	static inline
+	const message defaultMessage_ = 
+	(
+		message::CAP_CHANGED+
+		message::SIZE_CHANGED+
+		message::ITEM_INSERT+
+		message::ITEM_DELETE
+	);
 	
+	bool initMassDiameter()const;
 	
 	void zeroForce()
 	{
-		WARNING<<"fill contactTorque"<<END_WARNING;
+		contactForce_.fill(zero3);
-		//contactForce_.fill(zero3);
 	}
 
 	void zeroTorque()
 	{
-		WARNING<<"fill contactTorque"<<END_WARNING;
+		contactTorque_.fill(zero3);
-		//contactTorque_.fill(zero3);
 	}
+
 protected:
 
 	inline auto& dynPointStruct()
@@ -94,6 +103,12 @@ protected:
 		return dynPointStruct_.velocity();
 	}
 
+	inline 
+	auto& shapeIndex()
+	{
+		return shapeIndex_;
+	}
+
 public:
 
 	// type info
@@ -208,16 +223,12 @@ public:
 		return propertyId_;
 	}
 
-	/*inline const auto& shapeName()const{
-		return shapeName_;
-	}*/
-
-	/*inline auto& shapName(){
-		return shapeName_;
-	}*/
-
 	bool beforeIteration() override;
 
+	bool iterate()override;
+
+	bool afterIteration() override;
+
 	/*virtual
 	bool insertParticles
 	(
@@ -234,14 +245,17 @@ public:
 	virtual
 	const realx3PointField_D& rAcceleration() const = 0;
 
-	/*virtual
+	virtual
-	const realVector_D& boundingSphere()const = 0;*/
+	const realPointField_D& boundingSphere()const = 0;
 
 	virtual 
 	word shapeTypeName()const = 0;
 
 	virtual
-	void boundingSphereMinMax(real & minDiam, real& maxDiam)const = 0;
+	const shape& getShapes()const = 0;
+
+	virtual
+	void boundingSphereMinMax(real & minDiam, real& maxDiam)const;
 
 	
 

src/Particles/particles/shape.cpp

@@ -0,0 +1,76 @@
+#include "shape.hpp"
+
+
+bool pFlow::shape::findPropertyIds()
+{	
+	shapePropertyIds_.resize(numShapes());
+
+	ForAll( i, materialNames_)
+	{	
+		
+		if(uint32 propId; property_.nameToIndex(materialNames_[i], propId) )
+		{
+			shapePropertyIds_[i] = static_cast<uint8>(propId);
+		}
+		else
+		{
+			fatalErrorInFunction<<"Material name "<< materialNames_[i]<< 
+			"is not valid in dictionary "<<globalName()<<endl;
+			return false;
+		}
+	}
+    return false;
+}
+
+bool pFlow::shape::readFromDictionary()
+{
+	
+	materialNames_ = getVal<wordVector>("materials");
+
+	if(materialNames_.size() != numShapes() )
+	{
+		fatalErrorInFunction<<
+		"  number of elements in materials and names are not the same in "<< globalName()<<endl;
+		return false;
+	}
+
+	return true;
+}
+
+
+pFlow::shape::shape
+(
+    const word &fileName, 
+    repository *owner, 
+    const property &prop
+)
+:
+   baseShapeNames(fileName,owner),
+    property_(prop)
+{
+	
+	if( !readFromDictionary() )
+	{
+		fatalExit;
+	}
+
+	if(!findPropertyIds())
+	{
+		fatalExit;
+	}
+
+}
+
+bool pFlow::shape::writeToDict(dictionary &dict)const
+{
+	if(!baseShapeNames::writeToDict(dict))return false;
+
+	if( !dict.add("materials", materialNames_) )
+	{
+		fatalErrorInFunction<<
+		"  Error in writing materials to dictionary "<< dict.globalName()<<endl;
+		return false;
+	}
+
+	return true;
+}

src/Particles/particles/shape.hpp

@@ -0,0 +1,183 @@
+/*------------------------------- 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 __shape_hpp__
+#define __shape_hpp__
+
+#include "baseShapeNames.hpp"
+#include "property.hpp"
+
+namespace pFlow
+{
+
+
+class shape
+:
+	public baseShapeNames
+{
+private:
+
+	/// property of materials 
+  	const property& 	property_; 
+
+	/// list of property ids of shapes (index)
+	uint8Vector			shapePropertyIds_;
+
+	/// list of material names 
+	wordVector 			materialNames_;
+		
+	
+	bool findPropertyIds();
+
+	bool readFromDictionary();	
+
+protected:
+
+	
+	bool writeToDict(dictionary& dict)const override;
+
+public:
+    
+    TypeInfo("shape");
+
+	shape(
+		const word& fileName,
+		repository* owner,
+		const property& prop);
+
+    ~shape() override=default;
+
+	inline 
+	const auto& properties()const
+	{
+		return property_;
+	}
+
+    
+	inline
+    const wordVector& materialNames()const
+	{
+		return materialNames_;
+	}
+
+	inline
+	const auto& shapePropertyIds()const
+	{
+		return shapePropertyIds_;
+	}
+
+    
+	
+    inline
+    bool shapeNameToPropId(const word& name, int8& propId)const
+	{	
+		if(uint32 index; shapeNameToIndex(name, index))
+		{
+			propId = shapePropertyIds_[index];
+			return true;
+		}
+		return false;
+	}
+
+	inline
+    bool propIdValid(int8 propId)const
+	{
+		return static_cast<uint32>(propId) < property_.numMaterials();
+	}
+
+
+	inline
+    bool indexToDensity(uint32 index, real& rho)const
+	{
+		if(indexValid(index))
+			return property_.density(shapePropertyIds_[index], rho);
+		return false;
+	}
+
+    inline
+    real indexToDensity(uint32 index)const
+	{
+		if(indexValid(index))
+		{
+			return property_.density(shapePropertyIds_[index]);
+		}
+		
+		fatalExit;
+		return 0.0;
+	}
+
+	virtual 
+    real maxBoundingSphere()const = 0;
+
+    virtual 
+    real minBoundingSphere()const = 0;
+
+    virtual 
+    bool boundingDiameter(uint32 index, real& bDiam)const =0;
+
+    virtual
+    real boundingDiameter(uint32 index)const = 0;
+
+	virtual
+	realVector boundingDiameter()const = 0;
+
+    virtual 
+    bool mass(uint32 index, real& m)const = 0;
+
+    virtual 
+    real mass(uint32 index) const =0;
+
+	virtual
+	realVector mass()const =0;
+
+	virtual
+	realVector density()const =0;
+
+    virtual 
+    bool Inertia(uint32 index, real& I)const = 0;
+
+    virtual
+    real Inertia(uint32 index)const = 0;
+
+	virtual 
+	realVector Inertia()const=0;
+
+    virtual 
+    bool Inertia_xx(uint32 index, real& Ixx)const = 0;
+
+    virtual
+    real Inertial_xx(uint32 index)const =0;
+	
+    virtual 
+    bool Inertia_yy(uint32 index, real& Iyy)const = 0;
+
+    virtual 
+    real Inertial_yy(uint32 index)const = 0;
+
+    virtual 
+    bool Inertia_zz(uint32 index, real& Izz)const = 0;
+
+    virtual 
+    real Inertial_zz(uint32 index)const = 0;
+
+};
+
+}
+
+#endif //__shape_hpp__