yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3005
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Thu 2012-01-26 09:46:25 +0100
message:
  - fix a bracket mistake triggering ghosts definition after each insertion in periodic::triangulate()
  - add tests on boundary ids (<0?) before doing anything on them.
modified:
  pkg/dem/FlowEngine.cpp
  pkg/dem/FlowEngine.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp	2012-01-24 17:14:09 +0000
+++ pkg/dem/FlowEngine.cpp	2012-01-26 08:46:25 +0000
@@ -37,17 +37,6 @@
 void FlowEngine::action()
 {
 	if (!isActivated) return;
-	if (!solver) {
-		solver = shared_ptr<FlowSolver> (new FlowSolver);
-		first=true;
-		cerr <<"first = true"<<endl;
-		Update_Triangulation=false;
-		eps_vol_max=0.f;
-		Eps_Vol_Cumulative=0.f;
-		ReTrg=1;
-		retriangulationLastIter=0;
-	}
-
 	timingDeltas->start();
 
 	if (first) Build_Triangulation(P_zero,solver);
@@ -87,8 +76,8 @@
 		///End Compute flow and forces
 		timingDeltas->checkpoint("Applying Forces");
 	}
-// 	if ( scene->iter % PermuteInterval == 0 )
-// 	{ Update_Triangulation = true; }
+	if ( scene->iter % PermuteInterval == 0 )
+	{ Update_Triangulation = true; }
 
 	if (Update_Triangulation && !first) {
 		Build_Triangulation(solver);
@@ -106,26 +95,36 @@
 template<class Solver>
 void FlowEngine::BoundaryConditions(Solver& flow)
 {
-	flow->boundary ( flow->y_min_id ).flowCondition=Flow_imposed_BOTTOM_Boundary;
-	flow->boundary ( flow->y_max_id ).flowCondition=Flow_imposed_TOP_Boundary;
-	flow->boundary ( flow->x_max_id ).flowCondition=Flow_imposed_RIGHT_Boundary;
-	flow->boundary ( flow->x_min_id ).flowCondition=Flow_imposed_LEFT_Boundary;
-	flow->boundary ( flow->z_max_id ).flowCondition=Flow_imposed_FRONT_Boundary;
-	flow->boundary ( flow->z_min_id ).flowCondition=Flow_imposed_BACK_Boundary;
-
-	flow->boundary ( flow->y_min_id ).value=Pressure_BOTTOM_Boundary;
-	flow->boundary ( flow->y_max_id ).value=Pressure_TOP_Boundary;
-	flow->boundary ( flow->x_max_id ).value=Pressure_RIGHT_Boundary;
-	flow->boundary ( flow->x_min_id ).value=Pressure_LEFT_Boundary;
-	flow->boundary ( flow->z_max_id ).value=Pressure_FRONT_Boundary;
-	flow->boundary ( flow->z_min_id ).value=Pressure_BACK_Boundary;
-	
-	flow->boundary ( flow->y_min_id ).velocity = Vector3r::Zero();
-	flow->boundary ( flow->y_max_id ).velocity = topBoundaryVelocity;
-	flow->boundary ( flow->x_max_id ).velocity = Vector3r::Zero();
-	flow->boundary ( flow->x_min_id ).velocity = Vector3r::Zero();
-	flow->boundary ( flow->z_max_id ).velocity = Vector3r::Zero();
-	flow->boundary ( flow->z_min_id ).velocity = Vector3r::Zero();
+    if (flow->y_min_id>=0) {
+        flow->boundary ( flow->y_min_id ).flowCondition=Flow_imposed_BOTTOM_Boundary;
+        flow->boundary ( flow->y_min_id ).value=Pressure_BOTTOM_Boundary;
+        flow->boundary ( flow->y_min_id ).velocity = Vector3r::Zero();
+    }
+    if (flow->y_max_id>=0) {
+        flow->boundary ( flow->y_max_id ).flowCondition=Flow_imposed_TOP_Boundary;
+        flow->boundary ( flow->y_max_id ).value=Pressure_TOP_Boundary;
+        flow->boundary ( flow->y_max_id ).velocity = topBoundaryVelocity;
+    }
+    if (flow->x_max_id>=0) {
+        flow->boundary ( flow->x_max_id ).flowCondition=Flow_imposed_RIGHT_Boundary;
+        flow->boundary ( flow->x_max_id ).value=Pressure_RIGHT_Boundary;
+        flow->boundary ( flow->x_max_id ).velocity = Vector3r::Zero();
+    }
+    if (flow->x_min_id>=0) {
+        flow->boundary ( flow->x_min_id ).flowCondition=Flow_imposed_LEFT_Boundary;
+        flow->boundary ( flow->x_min_id ).value=Pressure_LEFT_Boundary;
+        flow->boundary ( flow->x_min_id ).velocity = Vector3r::Zero();
+    }
+    if (flow->z_max_id>=0) {
+        flow->boundary ( flow->z_max_id ).flowCondition=Flow_imposed_FRONT_Boundary;
+        flow->boundary ( flow->z_max_id ).value=Pressure_FRONT_Boundary;
+        flow->boundary ( flow->z_max_id ).velocity = Vector3r::Zero();
+    }
+    if (flow->z_min_id>=0) {
+        flow->boundary ( flow->z_min_id ).flowCondition=Flow_imposed_BACK_Boundary;
+        flow->boundary ( flow->z_min_id ).value=Pressure_BACK_Boundary;
+        flow->boundary ( flow->z_min_id ).velocity = Vector3r::Zero();
+    }
 }
 
 template<class Solver>
@@ -153,13 +152,14 @@
 	Update_Triangulation=true;
 }
 template<class Solver>
-void FlowEngine::clearImposedPressure(Solver& flow) { flow->imposedP.clear();}
+void FlowEngine::clearImposedPressure(Solver& flow) { flow->imposedP.clear(); flow->IPCells.clear();}
+
 template<class Solver>
 Real FlowEngine::getFlux(unsigned int cond,Solver& flow) {
 	if (cond>=flow->imposedP.size()) LOG_ERROR("Getting flux with cond higher than imposedP size.");
 	RTriangulation& Tri = flow->T[flow->currentTes].Triangulation();
 	double flux=0;
-	Cell_handle cell= Tri.locate(flow->imposedP[cond].first);
+	Cell_handle& cell= flow->IPCells[cond];
 	for (int ngb=0;ngb<4;ngb++) {
 		/*if (!cell->neighbor(ngb)->info().Pcondition)*/ flux+= cell->info().k_norm()[ngb]*(cell->info().p()-cell->neighbor(ngb)->info().p());}
 	return flux+cell->info().dv();
@@ -265,12 +265,12 @@
 	flow->z_min_id=wallBackId;
 	flow->z_max_id=wallFrontId;
 
-	flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
-	flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
-	flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
-	flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
-	flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
-	flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
+	if (flow->y_min_id>=0) flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
+	if (flow->y_max_id>=0) flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
+	if (flow->x_max_id>=0) flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
+	if (flow->x_min_id>=0) flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
+	if (flow->z_max_id>=0) flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
+	if (flow->z_min_id>=0) flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
 
 	//FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
 	flow->boundsIds[0]= &flow->x_min_id;
@@ -592,17 +592,6 @@
 	CGT::PeriodicCellInfo::hSize[1] = makeCgVect(scene->cell->hSize.col(1));
 	CGT::PeriodicCellInfo::hSize[2] = makeCgVect(scene->cell->hSize.col(2));
 	if (!isActivated) return;
-	if (!solver) {
-		solver = shared_ptr<CGT::PeriodicFlow> (new CGT::PeriodicFlow);
-		first=true;
-		cerr <<"first = true"<<endl;
-		Update_Triangulation=false;
-		eps_vol_max=0.f;
-		Eps_Vol_Cumulative=0.f;
-		ReTrg=1;
-		retriangulationLastIter=0;
-	}
-
 // 	timingDeltas->start();
 
 	if (first) Build_Triangulation(P_zero);
@@ -644,7 +633,7 @@
 // 	{ Update_Triangulation = true; }
 
 	if (Update_Triangulation && !first) {
-		Build_Triangulation(0);
+		Build_Triangulation(P_zero);
 		Update_Triangulation = false;}
 		
 // 	if (velocity_profile) /*flow->FluidVelocityProfile();*/solver->Average_Fluid_Velocity();
@@ -676,10 +665,6 @@
 		Real y = wpos[1];
 		Real z = wpos[2];
 		solver->T[solver->currentTes].insert(x, y, z, rad, id);
-#ifdef XVIEW
-		Vue1.SetSpheresColor(0.8,0.6,0.6,1);
-		Vue1.Dessine_Sphere(x,y,z, rad, 15);
-#endif
 		const Vector3r& cellSize(scene->cell->getSize());
 		wpos=scene->cell->unshearPt(wpos);
 		// traverse all periodic cells around the body, to see if any of them touches
@@ -697,23 +682,16 @@
 						Vector3r pt=scene->cell->shearPt(pos2);
 						Vertex_handle vh=solver->T[solver->currentTes].insert(pt[0],pt[1],pt[2],rad,id,false,id);
 						for (int k=0;k<3;k++) vh->info().period[k]=i[k];    
-#ifdef XVIEW
-                            Vue1.SetSpheresColor(0.1,0.5,0.8,1);
-                            Vue1.Dessine_Sphere(pt[0],pt[1],pt[2], rad, 15);
-#endif
 					}
 				}
-		cerr <<"number of vtx"<< solver->T[solver -> currentTes].Triangulation().number_of_vertices()<<endl;
-		// Define the ghost cells
-		Finite_cells_iterator cellend=solver->T[solver->currentTes].Triangulation().finite_cells_end();
-		for (Finite_cells_iterator cell=solver -> T[solver -> currentTes].Triangulation().finite_cells_begin(); cell!=cellend; cell++)
-			locateCell(cell);
-		solver -> viscousShearForces.resize(solver -> T[solver -> currentTes].max_id+1);
-#ifdef XVIEW
-    Vue1.Affiche();
-#endif
+// 		cerr <<"number of vtx"<< solver->T[solver -> currentTes].Triangulation().number_of_vertices()<<endl;
 	}
-     }
+    }
+	// Define the ghost cells
+	Finite_cells_iterator cellend=solver->T[solver->currentTes].Triangulation().finite_cells_end();
+	for (Finite_cells_iterator cell=solver -> T[solver -> currentTes].Triangulation().finite_cells_begin(); cell!=cellend; cell++)
+		locateCell(cell);
+	solver -> viscousShearForces.resize(solver -> T[solver -> currentTes].max_id+1);
 }
 
 
@@ -901,6 +879,7 @@
 	solver->x_min = 1000.0, solver->x_max = -10000.0, solver->y_min = 1000.0, solver->y_max = -10000.0, solver->z_min = 1000.0, solver->z_max = -10000.0;
 
 	AddBoundary (solver);
+	if (Debug) cout << endl << "Added boundaries------" << endl << endl;
 	Triangulate ();
 	if (Debug) cout << endl << "Tesselating------" << endl << endl;
 	solver->T[solver->currentTes].Compute();

=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp	2012-01-25 10:02:04 +0000
+++ pkg/dem/FlowEngine.hpp	2012-01-26 08:46:25 +0000
@@ -33,7 +33,7 @@
 	typedef FlowSolver::Cell_handle						Cell_handle;
 	typedef RTriangulation::Finite_edges_iterator				Finite_edges_iterator;
 	
-	private:
+	protected:
 		shared_ptr<FlowSolver> solver;
 
 	public :
@@ -263,7 +263,11 @@
 			((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
 			,,
 			wallTopId=wallBottomId=wallFrontId=wallBackId=wallLeftId=wallRightId=-1;
+// 			FlowEngine::solver=shared_ptr<FlowSolver>;
 			solver = shared_ptr<FlowSolver> (new FlowSolver);
+			Update_Triangulation=false;
+			eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
+			ReTrg=1;
 			,
 			.def("MeanVelocity",&PeriodicFlowEngine::MeanVelocity,"measure the mean velocity in the period")
 // 			.def("imposeFlux",&FlowEngine::_imposeFlux,(python::arg("pos"),python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")