change for predictor-corrector and adams-moulton3

src/Particles/SphereParticles/sphereParticles/sphereParticles.C

@@ -83,8 +83,15 @@ bool pFlow::sphereParticles::beforeIteration()
 	particles::beforeIteration();
 	
 	intPredictTimer_.start();
+	
+	//Info<<"before dyn predict"<<endInfo;
 		dynPointStruct_.predict(this->dt(), accelertion_);
+	//Info<<"after dyn predict"<<endInfo;
+
+	//Info<<"before revel predict"<<endInfo;
 		rVelIntegration_().predict(this->dt(),rVelocity_, rAcceleration_);
+	//Info<<"after rvel predict"<<endInfo;
+
 	intPredictTimer_.end();
 
 	return true;
@@ -93,7 +100,9 @@ bool pFlow::sphereParticles::beforeIteration()
 
 bool pFlow::sphereParticles::iterate() 
 {
+
 	accelerationTimer_.start();
+	//Info<<"before accelerationTimer_ "<<endInfo;
 		pFlow::sphereParticlesKernels::acceleration(
 			control().g(),
 			mass().deviceVectorAll(),
@@ -107,8 +116,11 @@ bool pFlow::sphereParticles::iterate()
 	accelerationTimer_.end();
 	
 	intCorrectTimer_.start();
+	//Info<<"before correct dyn "<<endInfo;
 		dynPointStruct_.correct(this->dt(), accelertion_);
+	//Info<<"after correct dyn "<<endInfo;
 		rVelIntegration_().correct(this->dt(), rVelocity_, rAcceleration_);
+	//Info<<"after correct rvel "<<endInfo;
 	intCorrectTimer_.end();
 	
 	return true;
@@ -285,6 +297,28 @@ pFlow::sphereParticles::sphereParticles(
 		fatalExit;
 	}
 
+	
+	if(rVelIntegration_->needSetInitialVals())
+	{
+		
+		auto [minInd, maxInd] = pStruct().activeRange();
+		int32IndexContainer indexHD(minInd, maxInd);
+	
+		auto n = indexHD.size();
+		auto index = indexHD.indicesHost();
+
+		realx3Vector rvel(n,RESERVE());
+		const auto hrVel = rVelocity_.hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			rvel.push_back( hrVel[index(i)]);
+		}
+		
+		Report(2)<< "Initializing the required vectors for rotational velocity integratoin\n "<<endReport;
+		rVelIntegration_->setInitialVals(indexHD, rvel);
+
+	}
 
 	if(!initializeParticles())
 	{
@@ -293,6 +327,32 @@ pFlow::sphereParticles::sphereParticles(
 	
 } 
 
+bool pFlow::sphereParticles::update(const eventMessage& msg) 
+{
+	
+	if(rVelIntegration_->needSetInitialVals())
+	{
+		
+		
+		auto indexHD = pStruct().insertedPointIndex();
+	
+		auto n = indexHD.size();
+		auto index = indexHD.indicesHost();
+
+		realx3Vector rvel(n,RESERVE());
+		const auto hrVel = rVelocity_.hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			rvel.push_back( hrVel[index(i)]);
+		}
+		
+		rVelIntegration_->setInitialVals(indexHD, rvel);
+
+	}	
+
+	return true;
+}
 
 bool pFlow::sphereParticles::insertParticles
 (
@@ -312,7 +372,7 @@ bool pFlow::sphereParticles::insertParticles
 
 	auto exclusionListAllPtr = getFieldObjectList(); 
 
-	auto newInsertedPtr = pStruct().insertPoints( position, exclusionListAllPtr());
+	auto newInsertedPtr = pStruct().insertPoints( position, setField, time(), exclusionListAllPtr());
 	
 	if(!newInsertedPtr)
 	{
@@ -337,10 +397,6 @@ bool pFlow::sphereParticles::insertParticles
 		return false;
 	}
 	
-	for(auto sfEntry:setField)
-	{
-		if(!sfEntry.setPointFieldSelectedAll( time(), newInserted, false ))return false;
-	}
 
 	auto activeR = this->activeRange();
 

src/Particles/SphereParticles/sphereParticles/sphereParticles.H

@@ -157,11 +157,8 @@ public:
 		shapes_.diameterMinMax(minDiam, maxDiam);
 	}
 
-	bool update(const eventMessage& msg) override
-	{
-		CONSUME_PARAM(msg);
-		return true;
-	}
+	bool update(const eventMessage& msg) override;
+	
 	
 	realx3PointField_D& rAcceleration() override
 	{

src/Particles/dynamicPointStructure/dynamicPointStructure.C

@@ -47,26 +47,28 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 				objectFile::READ_ALWAYS,
 				objectFile::WRITE_ALWAYS
 				),
-			pStruct_,
+			pStruct(),
 			zero3
 			)
 		)
 
 {
 
+	this->subscribe(pStruct());
+
 	Report(1)<< "Creating integration method "<<
 		greenText(integrationMethod_)<<" for dynamicPointStructure."<<endReport;
 
 	integrationPos_ = integration::create(
 		"pStructPosition",
 		time_.integration(),
-		pStruct_,
+		pStruct(),
 		integrationMethod_);
 	
 	integrationVel_ = integration::create(
 		"pStructVelocity",
 		time_.integration(),
-		pStruct_,
+		pStruct(),
 		integrationMethod_);
 
 	if( !integrationPos_ )
@@ -83,6 +85,36 @@ pFlow::dynamicPointStructure::dynamicPointStructure
 		fatalExit;
 	}
 
+	
+	if(!integrationPos_->needSetInitialVals()) return;
+
+	
+		
+	auto [minInd, maxInd] = pStruct().activeRange();
+	int32IndexContainer indexHD(minInd, maxInd);
+	
+	auto n = indexHD.size();
+	auto index = indexHD.indicesHost();
+
+	realx3Vector pos(n,RESERVE());
+	realx3Vector vel(n,RESERVE());
+	const auto hVel = velocity().hostVector();
+	const auto hPos = pStruct().pointPosition().hostVector();
+
+	for(auto i=0; i<n; i++)
+	{
+		pos.push_back( hPos[index(i)]);
+		vel.push_back( hVel[index(i)]);
+	}
+
+	//output<< "pos "<< pos<<endl;
+	//output<< "vel "<< vel<<endl;
+
+	Report(2)<< "Initializing the required vectors for position integratoin "<<endReport;
+	integrationPos_->setInitialVals(indexHD, pos);
+
+	Report(2)<< "Initializing the required vectors for velocity integratoin\n "<<endReport;
+	integrationVel_->setInitialVals(indexHD, vel);
 }
 
 bool pFlow::dynamicPointStructure::predict
@@ -91,7 +123,7 @@ bool pFlow::dynamicPointStructure::predict
 	realx3PointField_D& acceleration
 )
 {
-	auto& pos = pStruct_.pointPosition();
+	auto& pos = pStruct().pointPosition();
 
 	if(!integrationPos_().predict(dt, pos.VectorField(), velocity_.VectorField() ))return false;
 	if(!integrationVel_().predict(dt, velocity_.VectorField(), acceleration.VectorField()))return false;
@@ -105,11 +137,82 @@ bool pFlow::dynamicPointStructure::correct
 	realx3PointField_D& acceleration
 )
 {
-	auto& pos = pStruct_.pointPosition();
+	auto& pos = pStruct().pointPosition();
 	
 	if(!integrationPos_().correct(dt, pos.VectorField(), velocity_.VectorField() ))return false;
 	
 	if(!integrationVel_().correct(dt, velocity_.VectorField(), acceleration.VectorField()))return false;
 
 	return true;	
+}
+
+/*FUNCTION_H
+pFlow::uniquePtr<pFlow::int32IndexContainer> pFlow::dynamicPointStructure::insertPoints
+(
+	const realx3Vector& pos,
+	const List<eventObserver*>& exclusionList
+)
+{
+	auto newIndicesPtr = pointStructure::insertPoints(pos, exclusionList);
+
+	// no new point is inserted
+	if( !newIndicesPtr ) return newIndicesPtr;
+
+	if(!integrationPos_().needSetInitialVals()) return newIndicesPtr;
+
+	auto hVel = velocity_.hostVector();
+	auto n = newIndicesPtr().size();
+	auto index = newIndicesPtr().indicesHost();
+
+	realx3Vector velVec(n, RESERVE());
+	for(auto i=0; i<n; i++)
+	{
+		velVec.push_back( hVel[ index(i) ]);
+	}
+
+	integrationPos_->setInitialVals(newIndicesPtr(), pos );
+	integrationVel_->setInitialVals(newIndicesPtr(), velVec );
+	
+	return newIndicesPtr;
+
+}*/
+
+
+bool pFlow::dynamicPointStructure::update(
+	const eventMessage& msg) 
+{
+	if( msg.isInsert())
+	{
+		
+		if(!integrationPos_->needSetInitialVals())return true;
+
+		const auto indexHD = pStruct().insertedPointIndex();
+
+	
+		auto n = indexHD.size();
+
+		if(n==0) return true;
+
+		auto index = indexHD.indicesHost();
+
+		realx3Vector pos(n,RESERVE());
+		realx3Vector vel(n,RESERVE());
+		const auto hVel = velocity().hostVector();
+		const auto hPos = pStruct().pointPosition().hostVector();
+
+		for(auto i=0; i<n; i++)
+		{
+			pos.push_back( hPos[index(i)]);
+			vel.push_back( hVel[index(i)]);
+		}
+
+
+		
+		integrationPos_->setInitialVals(indexHD, pos);
+
+		integrationVel_->setInitialVals(indexHD, vel);		
+
+	}
+
+	return true;
 }

src/Particles/dynamicPointStructure/dynamicPointStructure.H

@@ -31,12 +31,15 @@ namespace pFlow
 {
 
 class dynamicPointStructure
+:
+	//public pointStructure
+public eventObserver
 {
 protected:
 
 	Time& time_;
 
-	const word 			integrationMethod_;
+	word 			integrationMethod_;
 
 	pointStructure& 	pStruct_;
 
@@ -52,8 +55,36 @@ public:
 
 	dynamicPointStructure(Time& time, const word& integrationMethod);
 
+	dynamicPointStructure(const dynamicPointStructure& ps) = default;
+
+	/*dynamicPointStructure(const dynamicPointStructure& ps):
+		pointStructure(ps),
+		time_(ps.time_),
+		integrationMethod_(ps.integrationMethod_),
+		velocity_(ps.velocity_),
+		integrationPos_(ps.integrationPos_->clone()),
+		integrationVel_(ps.integrationVel_->clone())
+	{
+
+	}*/
+
+		
+		// - no move construct 
+	dynamicPointStructure(dynamicPointStructure&&) = delete;
+
+	// - copy assignment 
+	//
+	//   should be changed, may causs undefined behavior 
+	//////////////////////////////////////////////////////////////
+	dynamicPointStructure& operator=(const dynamicPointStructure&) = default;
+
+	// - no move assignment 
+	dynamicPointStructure& operator=(dynamicPointStructure&&) = delete;
+
+	// - destructor 
 	virtual ~dynamicPointStructure() = default;
 
+
 	inline pointStructure& pStruct() 
 	{
 		return pStruct_;
@@ -79,6 +110,19 @@ public:
 	bool correct(real dt, realx3PointField_D& acceleration);
 
 
+	// - update data structure by inserting/setting new points 
+	//   Notifies all the fields in the registered list of data structure
+	//   and exclude the fields that re in the exclusionList
+	//   retrun nullptr if it fails 
+	/*FUNCTION_H
+	virtual uniquePtr<int32IndexContainer> insertPoints(
+		const realx3Vector& pos,
+		const List<eventObserver*>& exclusionList={nullptr}
+	)override;*/
+
+	bool update(const eventMessage& msg) override;
+	
+
 };
 
 }

src/Particles/particles/particles.C

@@ -31,10 +31,21 @@ pFlow::particles::particles
 	demParticles(control),
 	time_(control.time()),
 	integrationMethod_(integrationMethod),
+	/*dynPointStruct_(
+		time_.emplaceObject<dynamicPointStructure>(
+			objectFile(
+				pointStructureFile__,
+				"",
+				objectFile::READ_ALWAYS,
+				objectFile::WRITE_ALWAYS			
+				),
+			control.time(),
+			integrationMethod
+			)
+		),*/
 	dynPointStruct_(
 		control.time(),
-		integrationMethod
-		),
+		integrationMethod),
 	shapeName_(
 		control.time().emplaceObject<wordPointField>(
 			objectFile(