- added substepping

- added OBJ export

Author: janbender

Changelog.txt

@@ -1,3 +1,6 @@
+1.8.0
+	- added OBJ export
+	- added substepping
 	- updated GenericParameters
 	- bufixes
 	- fixed several compiler warnings

Demos/Common/DemoBase.h

@@ -114,6 +114,7 @@ namespace PBD
 		std::vector<unsigned int>& getSelectedRigidBodies() { return m_selectedBodies; }
 		bool getUseCache() const { return m_useCache; }
 		void setUseCache(bool val) { m_useCache = val; }
+		std::string getOutputPath() const { return m_outputPath; }
 
 		Utilities::SceneLoader::SceneData& getSceneData() { return m_scene; }
 	};

Demos/PositionBasedElasticRodsDemo/PositionBasedElasticRodsTSC.cpp

@@ -21,7 +21,7 @@ PositionBasedElasticRodsTSC::~PositionBasedElasticRodsTSC()
 void PositionBasedElasticRodsTSC::step(SimulationModel &model)
 {
  	TimeManager *tm = TimeManager::getCurrent ();
- 	const Real h = tm->getTimeStepSize();
+	const Real hOld = tm->getTimeStepSize();
 	PositionBasedElasticRodsModel &ermodel = (PositionBasedElasticRodsModel&)model;
 
 	//////////////////////////////////////////////////////////////////////////
@@ -33,86 +33,102 @@ void PositionBasedElasticRodsTSC::step(SimulationModel &model)
 	ParticleData &pg = ermodel.getGhostParticles();
 
 	const int numBodies = (int)rb.size();
-	#pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
-	{
-		#pragma omp for schedule(static) nowait
-		for (int i = 0; i < numBodies; i++)
- 		{ 
-			rb[i]->getLastPosition() = rb[i]->getOldPosition();
-			rb[i]->getOldPosition() = rb[i]->getPosition();
-			TimeIntegration::semiImplicitEuler(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getVelocity(), rb[i]->getAcceleration());
-			rb[i]->getLastRotation() = rb[i]->getOldRotation();
-			rb[i]->getOldRotation() = rb[i]->getRotation();
-			TimeIntegration::semiImplicitEulerRotation(h, rb[i]->getMass(), rb[i]->getInertiaTensorInverseW(), rb[i]->getRotation(), rb[i]->getAngularVelocity(), rb[i]->getTorque());
-			rb[i]->rotationUpdated();
- 		}
 
-		//////////////////////////////////////////////////////////////////////////
-		// particle model
-		//////////////////////////////////////////////////////////////////////////
-		#pragma omp for schedule(static) 
-		for (int i = 0; i < (int) pd.size(); i++)
+	Real h = hOld / (Real)m_subSteps;
+	tm->setTimeStepSize(h);
+	for (unsigned int step = 0; step < m_subSteps; step++)
+	{
+		#pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
 		{
-			pd.getLastPosition(i) = pd.getOldPosition(i);
-			pd.getOldPosition(i) = pd.getPosition(i);
-			pd.getVelocity(i) *= (1.0f - m_damping);
+			#pragma omp for schedule(static) nowait
+			for (int i = 0; i < numBodies; i++)
+ 			{ 
+				rb[i]->getLastPosition() = rb[i]->getOldPosition();
+				rb[i]->getOldPosition() = rb[i]->getPosition();
+				TimeIntegration::semiImplicitEuler(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getVelocity(), rb[i]->getAcceleration());
+				rb[i]->getLastRotation() = rb[i]->getOldRotation();
+				rb[i]->getOldRotation() = rb[i]->getRotation();
+				TimeIntegration::semiImplicitEulerRotation(h, rb[i]->getMass(), rb[i]->getInertiaTensorInverseW(), rb[i]->getRotation(), rb[i]->getAngularVelocity(), rb[i]->getTorque());
+				rb[i]->rotationUpdated();
+ 			}
+
+			//////////////////////////////////////////////////////////////////////////
+			// particle model
+			//////////////////////////////////////////////////////////////////////////
+			#pragma omp for schedule(static) 
+			for (int i = 0; i < (int) pd.size(); i++)
+			{
+				pd.getLastPosition(i) = pd.getOldPosition(i);
+				pd.getOldPosition(i) = pd.getPosition(i);
+				pd.getVelocity(i) *= (1.0f - m_damping);
 
-			TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i));
-		}
+				TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i));
+			}
 
-		for (int i = 0; i < (int)pg.size(); i++)
-		{
-			pg.getLastPosition(i) = pg.getOldPosition(i);
-			pg.getOldPosition(i) = pg.getPosition(i);
+			for (int i = 0; i < (int)pg.size(); i++)
+			{
+				pg.getLastPosition(i) = pg.getOldPosition(i);
+				pg.getOldPosition(i) = pg.getPosition(i);
 
-			pg.getVelocity(i) *= (1.0f - m_damping);
+				pg.getVelocity(i) *= (1.0f - m_damping);
 
-			TimeIntegration::semiImplicitEuler(h, pg.getMass(i), pg.getPosition(i), pg.getVelocity(i), pg.getAcceleration(i));
+				TimeIntegration::semiImplicitEuler(h, pg.getMass(i), pg.getPosition(i), pg.getVelocity(i), pg.getAcceleration(i));
+			}
 		}
-	}
 
-	positionConstraintProjection(model);
+		positionConstraintProjection(model);
 
-	#pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
-	{
- 		// Update velocities	
-		#pragma omp for schedule(static) nowait
-		for (int i = 0; i < numBodies; i++)
- 		{
-			if (m_velocityUpdateMethod == 0)
+		#pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
+		{
+ 			// Update velocities	
+			#pragma omp for schedule(static) nowait
+			for (int i = 0; i < numBodies; i++)
+ 			{
+				if (m_velocityUpdateMethod == 0)
+				{
+					TimeIntegration::velocityUpdateFirstOrder(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getOldPosition(), rb[i]->getVelocity());
+					TimeIntegration::angularVelocityUpdateFirstOrder(h, rb[i]->getMass(), rb[i]->getRotation(), rb[i]->getOldRotation(), rb[i]->getAngularVelocity());
+				}
+				else
+				{
+					TimeIntegration::velocityUpdateSecondOrder(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getOldPosition(), rb[i]->getLastPosition(), rb[i]->getVelocity());
+					TimeIntegration::angularVelocityUpdateSecondOrder(h, rb[i]->getMass(), rb[i]->getRotation(), rb[i]->getOldRotation(), rb[i]->getLastRotation(), rb[i]->getAngularVelocity());
+				}
+ 			}
+
+			// Update velocities	
+			#pragma omp for schedule(static) 
+			for (int i = 0; i < (int) pd.size(); i++)
 			{
-				TimeIntegration::velocityUpdateFirstOrder(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getOldPosition(), rb[i]->getVelocity());
-				TimeIntegration::angularVelocityUpdateFirstOrder(h, rb[i]->getMass(), rb[i]->getRotation(), rb[i]->getOldRotation(), rb[i]->getAngularVelocity());
+				if (m_velocityUpdateMethod == 0)
+					TimeIntegration::velocityUpdateFirstOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getVelocity(i));
+				else
+					TimeIntegration::velocityUpdateSecondOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getLastPosition(i), pd.getVelocity(i));
 			}
-			else
+	
+			#pragma omp for schedule(static) 
+			for (int i = 0; i < (int)pg.size(); i++)
 			{
-				TimeIntegration::velocityUpdateSecondOrder(h, rb[i]->getMass(), rb[i]->getPosition(), rb[i]->getOldPosition(), rb[i]->getLastPosition(), rb[i]->getVelocity());
-				TimeIntegration::angularVelocityUpdateSecondOrder(h, rb[i]->getMass(), rb[i]->getRotation(), rb[i]->getOldRotation(), rb[i]->getLastRotation(), rb[i]->getAngularVelocity());
+				if (m_velocityUpdateMethod == 0)
+					TimeIntegration::velocityUpdateFirstOrder(h, pg.getMass(i), pg.getPosition(i), pg.getOldPosition(i), pg.getVelocity(i));
+				else
+					TimeIntegration::velocityUpdateSecondOrder(h, pg.getMass(i), pg.getPosition(i), pg.getOldPosition(i), pg.getLastPosition(i), pg.getVelocity(i));
 			}
-			// update geometry
-			rb[i]->getGeometry().updateMeshTransformation(rb[i]->getPosition(), rb[i]->getRotationMatrix());
- 		}
+		}
+	}
+	h = hOld;
+	tm->setTimeStepSize(hOld);
 
+	#pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared)
+	{
 		// Update velocities	
-		#pragma omp for schedule(static) 
-		for (int i = 0; i < (int) pd.size(); i++)
-		{
-			if (m_velocityUpdateMethod == 0)
-				TimeIntegration::velocityUpdateFirstOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getVelocity(i));
-			else
-				TimeIntegration::velocityUpdateSecondOrder(h, pd.getMass(i), pd.getPosition(i), pd.getOldPosition(i), pd.getLastPosition(i), pd.getVelocity(i));
-		}
-	
-		#pragma omp for schedule(static) 
-		for (int i = 0; i < (int)pg.size(); i++)
+		#pragma omp for schedule(static) nowait
+		for (int i = 0; i < numBodies; i++)
 		{
-			if (m_velocityUpdateMethod == 0)
-				TimeIntegration::velocityUpdateFirstOrder(h, pg.getMass(i), pg.getPosition(i), pg.getOldPosition(i), pg.getVelocity(i));
-			else
-				TimeIntegration::velocityUpdateSecondOrder(h, pg.getMass(i), pg.getPosition(i), pg.getOldPosition(i), pg.getLastPosition(i), pg.getVelocity(i));
+			if (rb[i]->getMass() != 0.0)
+				rb[i]->getGeometry().updateMeshTransformation(rb[i]->getPosition(), rb[i]->getRotationMatrix());
 		}
-
-}
+	}
 
 	if (m_collisionDetection)
 		m_collisionDetection->collisionDetection(model);

Demos/SceneLoaderDemo/SceneLoaderDemo.cpp

@@ -39,6 +39,7 @@ void buildModel ();
 void readScene(const bool readFile);
 void render ();
 void reset();
+void exportOBJ();
 void TW_CALL setContactTolerance(const void *value, void *clientData);
 void TW_CALL getContactTolerance(void *value, void *clientData);
 void TW_CALL setContactStiffnessRigidBody(const void *value, void *clientData);
@@ -60,6 +61,10 @@ CubicSDFCollisionDetection cd;
 short clothSimulationMethod = 2;
 short solidSimulationMethod = 2;
 short bendingMethod = 2;
+bool enableExportOBJ = false;
+unsigned int exportFPS = 25;
+Real nextFrameTime = 0.0;
+unsigned int frameCounter = 1;
 
 
 // main 
@@ -129,6 +134,9 @@ void initParameters()
 
 	MiniGL::initTweakBarParameters();
 
+	TwAddVarRW(MiniGL::getTweakBar(), "ExportOBJ", TW_TYPE_BOOL32, &enableExportOBJ, " label='Export OBJ'");
+	TwAddVarRW(MiniGL::getTweakBar(), "ExportFPS", TW_TYPE_UINT32, &exportFPS, " label='Export FPS'");
+
 	TweakBarParameters::createParameterGUI();
 	TweakBarParameters::createParameterObjectGUI(base);
 	TweakBarParameters::createParameterObjectGUI(Simulation::getCurrent());
@@ -138,6 +146,7 @@ void initParameters()
 
 void reset()
 {
+	const Real h = TimeManager::getCurrent()->getTimeStepSize();
 	Utilities::Timing::printAverageTimes();
 	Utilities::Timing::reset();
 
@@ -147,7 +156,11 @@ void reset()
 	Simulation::getCurrent()->getModel()->cleanup();
 	Simulation::getCurrent()->getTimeStep()->getCollisionDetection()->cleanup();
 
+	nextFrameTime = 0.0;
+	frameCounter = 1;
+
 	readScene(false);
+	TimeManager::getCurrent()->setTimeStepSize(h);
 }
 
 void timeStep ()
@@ -167,6 +180,8 @@ void timeStep ()
 		START_TIMING("SimStep");
 		Simulation::getCurrent()->getTimeStep()->step(*model);
 		STOP_TIMING_AVG;
+
+		exportOBJ();
 	}
 
 	// Update visualization models
@@ -1089,6 +1104,101 @@ void readScene(const bool readFile)
 	}
 }
 
+void exportMeshOBJ(const std::string &exportFileName, const unsigned int nVert, const Vector3r* pos, const unsigned int nTri, const unsigned int* faces)
+{
+	// Open the file
+	std::ofstream outfile(exportFileName);
+	if (!outfile)
+	{
+		LOG_WARN << "Cannot open a file to save OBJ mesh.";
+		return;
+	}
+
+	// Header
+	outfile << "# Created by the PositionBasedDynamics library\n";
+	outfile << "g default\n";
+
+	// Vertices
+	{
+		for (auto j = 0u; j < nVert; j++)
+		{
+			const Vector3r& x = pos[j];
+			outfile << "v " << x[0] << " " << x[1] << " " << x[2] << "\n";
+		}
+	}
+
+	// faces
+	{
+		for (auto j = 0; j < nTri; j++)
+		{
+			outfile << "f " << faces[3 * j + 0] + 1 << " " << faces[3 * j + 1] + 1 << " " << faces[3 * j + 2] + 1 << "\n";
+		}
+	}
+	outfile.close();
+}
+
+
+void exportOBJ()
+{
+	if (!enableExportOBJ)
+		return;
+
+	if (TimeManager::getCurrent()->getTime() < nextFrameTime)
+		return;
+
+	nextFrameTime += 1.0 / (Real)exportFPS;
+
+	//////////////////////////////////////////////////////////////////////////
+	// rigid bodies
+	//////////////////////////////////////////////////////////////////////////
+
+	std::string exportPath = base->getOutputPath() + "/export";
+	FileSystem::makeDirs(exportPath);
+
+	SimulationModel* model = Simulation::getCurrent()->getModel();
+	const ParticleData& pd = model->getParticles();
+	for (unsigned int i = 0; i < model->getTriangleModels().size(); i++)
+	{
+		const IndexedFaceMesh& mesh = model->getTriangleModels()[i]->getParticleMesh();
+		const unsigned int offset = model->getTriangleModels()[i]->getIndexOffset();
+		const Vector3r* x = model->getParticles().getVertices()->data();
+		
+		std::string fileName = "triangle_model";
+		fileName = fileName + std::to_string(i) + "_" + std::to_string(frameCounter) + ".obj";
+		std::string exportFileName = FileSystem::normalizePath(exportPath + "/" + fileName);
+
+		exportMeshOBJ(exportFileName, mesh.numVertices(), &x[offset], mesh.numFaces(), mesh.getFaces().data());
+	}
+
+	for (unsigned int i = 0; i < model->getTetModels().size(); i++)
+	{
+		const IndexedFaceMesh& mesh = model->getTetModels()[i]->getVisMesh();
+		const unsigned int offset = model->getTetModels()[i]->getIndexOffset();
+		const Vector3r* x = model->getTetModels()[i]->getVisVertices().getVertices()->data();
+
+		std::string fileName = "tet_model";
+		fileName = fileName + std::to_string(i) + "_" + std::to_string(frameCounter) + ".obj";
+		std::string exportFileName = FileSystem::normalizePath(exportPath + "/" + fileName);
+
+		exportMeshOBJ(exportFileName, mesh.numVertices(), x, mesh.numFaces(), mesh.getFaces().data());
+	}
+
+	for (unsigned int i = 0; i < model->getRigidBodies().size(); i++)
+	{
+		const IndexedFaceMesh& mesh = model->getRigidBodies()[i]->getGeometry().getMesh();
+		const Vector3r *x = model->getRigidBodies()[i]->getGeometry().getVertexData().getVertices()->data();
+
+		std::string fileName = "rigid_body";
+		fileName = fileName + std::to_string(i) + "_" + std::to_string(frameCounter) + ".obj";
+		std::string exportFileName = FileSystem::normalizePath(exportPath + "/" + fileName);
+
+		exportMeshOBJ(exportFileName, mesh.numVertices(), x, mesh.numFaces(), mesh.getFaces().data());
+	}
+
+
+	frameCounter++;
+}
+
 void TW_CALL setContactTolerance(const void *value, void *clientData)
 {
 	const Real val = *(const Real *)(value);
@@ -1158,3 +1268,4 @@ void TW_CALL getSolidSimulationMethod(void *value, void *clientData)
 {
 	*(short *)(value) = *((short*)clientData);
 }
+

Demos/StiffRodsDemos/StiffRodsSceneLoader.cpp

@@ -524,8 +524,6 @@ void PBD::StiffRodsSceneLoader::readStiffRodsScene(const std::string &fileName,
 		//////////////////////////////////////////////////////////////////////////
 		sceneData.m_timeStepSize = static_cast<Real>(0.005);
 		sceneData.m_gravity = Vector3r(0, -static_cast<Real>(9.81), 0);
-		sceneData.m_maxIter = 5;
-		sceneData.m_maxIterVel = 5;
 		sceneData.m_velocityUpdateMethod = 0;
 		sceneData.m_triangleModelSimulationMethod = -1;
 		sceneData.m_triangleModelBendingMethod = -1;

Demos/Visualization/MiniGL.cpp

@@ -466,7 +466,8 @@ void MiniGL::display ()
 		scenefunc ();
 
 	TwDraw();  // draw the tweak bar(s)
-	glutSwapBuffers();
+	//glutSwapBuffers();
+	glFlush();
 }
 
 void MiniGL::init(int argc, char **argv, const int width, const int height, const int posx, const int posy, const char *name)
@@ -478,7 +479,7 @@ void MiniGL::init(int argc, char **argv, const int width, const int height, cons
 	scenefunc = NULL;
 
 	glutInit (&argc, argv);
-	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
+	glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB | GLUT_DEPTH);
 
 	atexit(destroy);