diff --git a/Changelog.txt b/Changelog.txt index b73f441e..d25de4f5 100644 --- a/Changelog.txt +++ b/Changelog.txt @@ -1,3 +1,6 @@ +1.8.0 + - added OBJ export + - added substepping - updated GenericParameters - bufixes - fixed several compiler warnings diff --git a/Demos/Common/DemoBase.h b/Demos/Common/DemoBase.h index c504602f..d85b6173 100644 --- a/Demos/Common/DemoBase.h +++ b/Demos/Common/DemoBase.h @@ -114,6 +114,7 @@ namespace PBD std::vector& 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; } }; diff --git a/Demos/PositionBasedElasticRodsDemo/PositionBasedElasticRodsTSC.cpp b/Demos/PositionBasedElasticRodsDemo/PositionBasedElasticRodsTSC.cpp index 48f9ccf5..632949c5 100644 --- a/Demos/PositionBasedElasticRodsDemo/PositionBasedElasticRodsTSC.cpp +++ b/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); diff --git a/Demos/SceneLoaderDemo/SceneLoaderDemo.cpp b/Demos/SceneLoaderDemo/SceneLoaderDemo.cpp index 635cd2f1..b1f390a6 100644 --- a/Demos/SceneLoaderDemo/SceneLoaderDemo.cpp +++ b/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); } + diff --git a/Demos/StiffRodsDemos/StiffRodsSceneLoader.cpp b/Demos/StiffRodsDemos/StiffRodsSceneLoader.cpp index 639d6eb6..bde42b47 100644 --- a/Demos/StiffRodsDemos/StiffRodsSceneLoader.cpp +++ b/Demos/StiffRodsDemos/StiffRodsSceneLoader.cpp @@ -524,8 +524,6 @@ void PBD::StiffRodsSceneLoader::readStiffRodsScene(const std::string &fileName, ////////////////////////////////////////////////////////////////////////// sceneData.m_timeStepSize = static_cast(0.005); sceneData.m_gravity = Vector3r(0, -static_cast(9.81), 0); - sceneData.m_maxIter = 5; - sceneData.m_maxIterVel = 5; sceneData.m_velocityUpdateMethod = 0; sceneData.m_triangleModelSimulationMethod = -1; sceneData.m_triangleModelBendingMethod = -1; diff --git a/Demos/Visualization/MiniGL.cpp b/Demos/Visualization/MiniGL.cpp index 4704fd28..141bec1c 100644 --- a/Demos/Visualization/MiniGL.cpp +++ b/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); diff --git a/README.md b/README.md index fec882ea..66517494 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,7 @@ Furthermore we use our own library: ## Build Instructions This project is based on [CMake](https://cmake.org/). Simply generate project, Makefiles, etc. using [CMake](https://cmake.org/) and compile the project with the compiler of your choice. The code was tested with the following configurations: -- Windows 10 64-bit, CMake 3.9.5, Visual Studio 2017 +- Windows 10 64-bit, CMake 3.9.5, Visual Studio 2019 - Debian 9 64-bit, CMake 3.12.3, GCC 6.3.0. Note: Please use a 64-bit target on a 64-bit operating system. 32-bit builds on a 64-bit OS are not supported. @@ -36,6 +36,8 @@ http://www.interactive-graphics.de/PositionBasedDynamics/doc/html ## Latest Important Changes +* added OBJ export +* added substepping * added DamperJoint * improved implementation of slider and hinge joints/motors * Crispin Deul added the implementation of his paper Deul, Kugelstadt, Weiler, Bender, "Direct Position-Based Solver for Stiff Rods", Computer Graphics Forum 2018 and a corresponding demo @@ -107,7 +109,7 @@ The following videos were generated using the PositionBasedDynamics library: ## Screenshots - + ![Cloth demo](http://www.interactive-graphics.de/j_images/PositionBasedDynamics.jpg "Cloth demo") ## References diff --git a/Simulation/ParticleData.h b/Simulation/ParticleData.h index d6c24715..7064249f 100644 --- a/Simulation/ParticleData.h +++ b/Simulation/ParticleData.h @@ -255,6 +255,11 @@ namespace PBD return (unsigned int) m_x.size(); } + FORCE_INLINE const std::vector* getVertices() + { + return &m_x; + } + /** Resize the array containing the particle data. */ FORCE_INLINE void resize(const unsigned int newSize) diff --git a/Simulation/Simulation.cpp b/Simulation/Simulation.cpp index 4c845c4d..7a35bf73 100644 --- a/Simulation/Simulation.cpp +++ b/Simulation/Simulation.cpp @@ -84,7 +84,6 @@ void Simulation::reset() m_timeStep->reset(); TimeManager::getCurrent()->setTime(static_cast(0.0)); - TimeManager::getCurrent()->setTimeStepSize(static_cast(0.005)); } void Simulation::setSimulationMethod(const int val) diff --git a/Simulation/TimeStepController.cpp b/Simulation/TimeStepController.cpp index fae6e856..981fc28c 100644 --- a/Simulation/TimeStepController.cpp +++ b/Simulation/TimeStepController.cpp @@ -12,6 +12,7 @@ using namespace GenParam; // int TimeStepController::SOLVER_ITERATIONS = -1; // int TimeStepController::SOLVER_ITERATIONS_V = -1; +int TimeStepController::NUM_SUB_STEPS = -1; int TimeStepController::MAX_ITERATIONS = -1; int TimeStepController::MAX_ITERATIONS_V = -1; int TimeStepController::VELOCITY_UPDATE_METHOD = -1; @@ -24,8 +25,9 @@ TimeStepController::TimeStepController() m_velocityUpdateMethod = 0; m_iterations = 0; m_iterationsV = 0; - m_maxIterations = 5; + m_maxIterations = 1; m_maxIterationsV = 5; + m_subSteps = 5; m_collisionDetection = NULL; } @@ -42,6 +44,11 @@ void TimeStepController::initParameters() // setDescription(SOLVER_ITERATIONS, "Iterations required by the solver."); // getParameter(SOLVER_ITERATIONS)->setReadOnly(true); + NUM_SUB_STEPS = createNumericParameter("subSteps", "# sub steps", &m_subSteps); + setGroup(NUM_SUB_STEPS, "PBD"); + setDescription(NUM_SUB_STEPS, "Number of sub steps of the solver."); + static_cast*>(getParameter(NUM_SUB_STEPS))->setMinValue(1); + MAX_ITERATIONS = createNumericParameter("maxIterations", "Max. iterations", &m_maxIterations); setGroup(MAX_ITERATIONS, "PBD"); setDescription(MAX_ITERATIONS, "Maximal number of iterations of the solver."); @@ -69,7 +76,7 @@ void TimeStepController::step(SimulationModel &model) { START_TIMING("simulation step"); TimeManager *tm = TimeManager::getCurrent (); - const Real h = tm->getTimeStepSize(); + const Real hOld = tm->getTimeStepSize(); ////////////////////////////////////////////////////////////////////////// // rigid body model @@ -80,46 +87,93 @@ void TimeStepController::step(SimulationModel &model) OrientationData &od = model.getOrientations(); 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); - TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i)); + #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); + TimeIntegration::semiImplicitEuler(h, pd.getMass(i), pd.getPosition(i), pd.getVelocity(i), pd.getAcceleration(i)); + } + + ////////////////////////////////////////////////////////////////////////// + // orientation model + ////////////////////////////////////////////////////////////////////////// + #pragma omp for schedule(static) + for (int i = 0; i < (int)od.size(); i++) + { + od.getLastQuaternion(i) = od.getOldQuaternion(i); + od.getOldQuaternion(i) = od.getQuaternion(i); + TimeIntegration::semiImplicitEulerRotation(h, od.getMass(i), od.getInvMass(i) * Matrix3r::Identity() ,od.getQuaternion(i), od.getVelocity(i), Vector3r(0,0,0)); + } } - ////////////////////////////////////////////////////////////////////////// - // orientation model - ////////////////////////////////////////////////////////////////////////// - #pragma omp for schedule(static) - for (int i = 0; i < (int)od.size(); i++) + START_TIMING("position constraints projection"); + positionConstraintProjection(model); + STOP_TIMING_AVG; + + #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared) { - od.getLastQuaternion(i) = od.getOldQuaternion(i); - od.getOldQuaternion(i) = od.getQuaternion(i); - TimeIntegration::semiImplicitEulerRotation(h, od.getMass(i), od.getInvMass(i) * Matrix3r::Identity() ,od.getQuaternion(i), od.getVelocity(i), Vector3r(0,0,0)); + // 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++) + { + 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)); + } + + // Update velocites of orientations + #pragma omp for schedule(static) + for (int i = 0; i < (int)od.size(); i++) + { + if (m_velocityUpdateMethod == 0) + TimeIntegration::angularVelocityUpdateFirstOrder(h, od.getMass(i), od.getQuaternion(i), od.getOldQuaternion(i), od.getVelocity(i)); + else + TimeIntegration::angularVelocityUpdateSecondOrder(h, od.getMass(i), od.getQuaternion(i), od.getOldQuaternion(i), od.getLastQuaternion(i), od.getVelocity(i)); + } } } - - START_TIMING("position constraints projection"); - positionConstraintProjection(model); - STOP_TIMING_AVG; + h = hOld; + tm->setTimeStepSize(hOld); #pragma omp parallel if(numBodies > MIN_PARALLEL_SIZE) default(shared) { @@ -127,42 +181,12 @@ void TimeStepController::step(SimulationModel &model) #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 geometry if (rb[i]->getMass() != 0.0) rb[i]->getGeometry().updateMeshTransformation(rb[i]->getPosition(), rb[i]->getRotationMatrix()); } - - // 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)); - } - - // Update velocites of orientations - #pragma omp for schedule(static) - for (int i = 0; i < (int)od.size(); i++) - { - if (m_velocityUpdateMethod == 0) - TimeIntegration::angularVelocityUpdateFirstOrder(h, od.getMass(i), od.getQuaternion(i), od.getOldQuaternion(i), od.getVelocity(i)); - else - TimeIntegration::angularVelocityUpdateSecondOrder(h, od.getMass(i), od.getQuaternion(i), od.getOldQuaternion(i), od.getLastQuaternion(i), od.getVelocity(i)); - } } + if (m_collisionDetection) { START_TIMING("collision detection"); @@ -221,8 +245,8 @@ void TimeStepController::reset() { m_iterations = 0; m_iterationsV = 0; - m_maxIterations = 5; - m_maxIterationsV = 5; + //m_maxIterations = 5; + //m_maxIterationsV = 5; } void TimeStepController::positionConstraintProjection(SimulationModel &model) diff --git a/Simulation/TimeStepController.h b/Simulation/TimeStepController.h index 7d88f3bc..dfb8b41a 100644 --- a/Simulation/TimeStepController.h +++ b/Simulation/TimeStepController.h @@ -13,6 +13,7 @@ namespace PBD public: // static int SOLVER_ITERATIONS; // static int SOLVER_ITERATIONS_V; + static int NUM_SUB_STEPS; static int MAX_ITERATIONS; static int MAX_ITERATIONS_V; static int VELOCITY_UPDATE_METHOD; @@ -24,6 +25,7 @@ namespace PBD int m_velocityUpdateMethod; unsigned int m_iterations; unsigned int m_iterationsV; + unsigned int m_subSteps; unsigned int m_maxIterations; unsigned int m_maxIterationsV; diff --git a/Utils/SceneLoader.cpp b/Utils/SceneLoader.cpp index 1577a7f3..8ca59d9f 100644 --- a/Utils/SceneLoader.cpp +++ b/Utils/SceneLoader.cpp @@ -35,8 +35,6 @@ void SceneLoader::readScene(const std::string &fileName, SceneData &sceneData) ////////////////////////////////////////////////////////////////////////// sceneData.m_timeStepSize = 0.005; sceneData.m_gravity = Vector3r(0, -9.81, 0); - sceneData.m_maxIter = 5; - sceneData.m_maxIterVel = 5; sceneData.m_velocityUpdateMethod = 0; sceneData.m_triangleModelSimulationMethod = -1; sceneData.m_triangleModelBendingMethod = -1; @@ -169,8 +167,6 @@ void SceneLoader::readSimulation(const nlohmann::json &j, const std::string &key readValue(child, "timeStepSize", sceneData.m_timeStepSize); readVector(child, "gravity", sceneData.m_gravity); - readValue(child, "maxIter", sceneData.m_maxIter); - readValue(child, "maxIterVel", sceneData.m_maxIterVel); readValue(child, "velocityUpdateMethod", sceneData.m_velocityUpdateMethod); readValue(child, "triangleModelSimulationMethod", sceneData.m_triangleModelSimulationMethod); readValue(child, "triangleModelBendingMethod", sceneData.m_triangleModelBendingMethod); diff --git a/Utils/SceneLoader.h b/Utils/SceneLoader.h index 8514d66f..de8fd6b0 100644 --- a/Utils/SceneLoader.h +++ b/Utils/SceneLoader.h @@ -196,8 +196,6 @@ namespace Utilities Real m_contactStiffnessParticleRigidBody; int m_velocityUpdateMethod; - unsigned int m_maxIter; - unsigned int m_maxIterVel; Real m_cloth_stiffness; Real m_cloth_bendingStiffness; diff --git a/data/scenes/ArmadilloCollisionScene.json b/data/scenes/ArmadilloCollisionScene.json index 8adb7e8d..22e856fb 100644 --- a/data/scenes/ArmadilloCollisionScene.json +++ b/data/scenes/ArmadilloCollisionScene.json @@ -2,9 +2,11 @@ "Name": "DeformableSolidCollisionScene1", "Simulation": { - "timeStepSize": 0.005, - "maxIter" : 5, - "maxIterVel" : 5, + "timeStepSize": 0.01, + "numberOfStepsPerRenderUpdate": 4, + "subSteps": 5, + "maxIterations" : 1, + "maxIterationsV" : 5, "velocityUpdateMethod" : 0, "contactTolerance": 0.0, "tetModelSimulationMethod": 2, diff --git a/data/scenes/CarScene.json b/data/scenes/CarScene.json index 73a61c56..927451a7 100644 --- a/data/scenes/CarScene.json +++ b/data/scenes/CarScene.json @@ -530,10 +530,11 @@ -2, 0 ], - "maxIter": 50, - "maxIterVel": 50, + "maxIterations": 1, + "maxIterationsV": 5, + "subSteps": 5, "numberOfStepsPerRenderUpdate": 10, - "timeStepSize": 0.005, + "timeStepSize": 0.01, "triangleModelBendingMethod": 2, "triangleModelSimulationMethod": 2, "velocityUpdateMethod": 0 diff --git a/data/scenes/ClothCollisionScene.json b/data/scenes/ClothCollisionScene.json index f16d1580..9a9eacaf 100644 --- a/data/scenes/ClothCollisionScene.json +++ b/data/scenes/ClothCollisionScene.json @@ -2,8 +2,9 @@ "Name": "ClothCollisionScene1", "Simulation": { - "timeStepSize": 0.0025, - "maxIter" : 5, + "timeStepSize": 0.005, + "maxIterations" : 1, + "subSteps": 2, "maxIterVel" : 5, "velocityUpdateMethod" : 0, "contactTolerance": 0.05, @@ -12,7 +13,7 @@ "contactStiffnessRigidBody" : 1.0, "contactStiffnessParticleRigidBody": 100.0, "cloth_stiffness": 1.0, - "cloth_bendingStiffness": 0.005, + "cloth_bendingStiffness": 0.01, "cloth_xxStiffness": 1.0, "cloth_yyStiffness": 1.0, "cloth_xyStiffness": 1.0,