diff --git a/Base/DynamicalModels/ObjectModels/imstkDynamicalModel.h b/Base/DynamicalModels/ObjectModels/imstkDynamicalModel.h
index 55f08614647b1dd142386110657a2901b42faefb..43873c2f67dd2910fdbdf91ea74d2d70c9705ade 100644
--- a/Base/DynamicalModels/ObjectModels/imstkDynamicalModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkDynamicalModel.h
@@ -119,6 +119,11 @@ public:
///
virtual void updatePhysicsGeometry(){};
+ ///
+ /// \brief Initialize the dynamical model
+ ///
+ virtual bool initialize() = 0;
+
protected:
DynamicalModelType m_type; ///> Mathematical model type
diff --git a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
index 80671c7665cc52e6989e3dc2cbbb99069326cb5a..5c90e10f9dbe24c85d64bdb015ff634c26f6ee4b 100644
--- a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
+++ b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.cpp
@@ -89,28 +89,55 @@ FEMDeformableBodyModel::configure(const std::string& configFileName)
m_forceModelConfiguration = std::make_shared<ForceModelConfig>(configFileName);
}
-void
-FEMDeformableBodyModel::initialize(std::shared_ptr<VolumetricMesh> physicsMesh)
+bool
+FEMDeformableBodyModel::initialize()
{
- this->setModelGeometry(physicsMesh);
-
// prerequisite of for successfully initializing
if (!m_forceModelGeometry || !m_forceModelConfiguration)
{
LOG(WARNING) << "DeformableBodyModel::initialize: Physics mesh or force model configuration not set yet!";
- return;
+ return false;
}
+ auto physicsMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(this->getModelGeometry());
m_vegaPhysicsMesh = physicsMesh->getAttachedVegaMesh();
- this->initializeForceModel();
- this->initializeMassMatrix();
- this->initializeDampingMatrix();
- this->initializeTangentStiffness();
+ if (!this->initializeForceModel())
+ {
+ return false;
+ }
+ ;
+ if (!this->initializeMassMatrix())
+ {
+ return false;
+ }
+ ;
+ if (!this->initializeDampingMatrix())
+ {
+ return false;
+ }
+ ;
+ if (!this->initializeTangentStiffness())
+ {
+ return false;
+ }
+ ;
this->loadInitialStates();
- this->loadBoundaryConditions();
- this->initializeGravityForce();
- this->initializeExplicitExternalForces();
+ if (!this->loadBoundaryConditions())
+ {
+ return false;
+ }
+ ;
+ if (!this->initializeGravityForce())
+ {
+ return false;
+ }
+ ;
+ if (!this->initializeExplicitExternalForces())
+ {
+ return false;
+ }
+ ;
m_Feff.resize(m_numDOF);
m_Finternal.resize(m_numDOF);
@@ -119,18 +146,25 @@ FEMDeformableBodyModel::initialize(std::shared_ptr<VolumetricMesh> physicsMesh)
m_Fcontact.setConstant(0.0);
m_qSol.resize(m_numDOF);
m_qSol.setConstant(0.0);
+
+ return true;
}
void
FEMDeformableBodyModel::loadInitialStates()
{
+ if (m_numDOF == 0)
+ {
+ LOG(WARNING) << "FEMDeformableBodyModel::loadInitialStates() - Num. of degree of freedom is zero!";
+ }
+
// For now the initial states are set to zero
m_initialState = std::make_shared<kinematicState>(m_numDOF);
m_previousState = std::make_shared<kinematicState>(m_numDOF);
m_currentState = std::make_shared<kinematicState>(m_numDOF);
}
-void
+bool
FEMDeformableBodyModel::loadBoundaryConditions()
{
auto fileName = m_forceModelConfiguration->getStringOptionsMap().at("fixedDOFFilename");
@@ -138,7 +172,7 @@ FEMDeformableBodyModel::loadBoundaryConditions()
if (fileName.empty())
{
LOG(WARNING) << "DeformableBodyModel::loadBoundaryConditions: The external boundary conditions file name is empty";
- return;
+ return false;
}
else
{
@@ -147,7 +181,6 @@ FEMDeformableBodyModel::loadBoundaryConditions()
if (file.peek() == std::ifstream::traits_type::eof())
{
LOG(INFO) << "DeformableBodyModel::loadBoundaryConditions: The external boundary conditions file is empty";
- return;
}
if (file.is_open())
@@ -165,6 +198,7 @@ FEMDeformableBodyModel::loadBoundaryConditions()
{
LOG(WARNING) << "FEMDeformableBodyModel::loadBoundaryConditions(): " <<
"The boundary condition node id provided is greater than number of nodes and hence excluded!!";
+ return false;
}
}
@@ -174,11 +208,13 @@ FEMDeformableBodyModel::loadBoundaryConditions()
else
{
LOG(WARNING) << "DeformableBodyModel::loadBoundaryConditions: Could not open boundary conditions file!";
+ return false;
}
}
+ return true;
}
-void
+bool
FEMDeformableBodyModel::initializeForceModel()
{
const double g = m_forceModelConfiguration->getFloatsOptionsMap().at("gravity");
@@ -215,16 +251,19 @@ FEMDeformableBodyModel::initializeForceModel()
default:
LOG(WARNING) << "DeformableBodyModel::initializeForceModel: Unknown force model type";
+ return false;
} //switch
+
+ return true;
}
-void
+bool
FEMDeformableBodyModel::initializeMassMatrix(const bool saveToDisk /*= false*/)
{
if (!m_forceModelGeometry)
{
LOG(WARNING) << "DeformableBodyModel::initializeMassMatrix Force model geometry not set!";
- return;
+ return false;
}
vega::SparseMatrix *vegaMatrix;
@@ -235,9 +274,11 @@ FEMDeformableBodyModel::initializeMassMatrix(const bool saveToDisk /*= false*/)
this->initializeEigenMatrixFromVegaMatrix(*vegaMatrix, m_M);
// TODO: Add option to save mass matrix to file
+
+ return true;
}
-void
+bool
FEMDeformableBodyModel::initializeDampingMatrix()
{
auto dampingLaplacianCoefficient = m_forceModelConfiguration->getFloatsOptionsMap().at("dampingLaplacianCoefficient");
@@ -248,13 +289,13 @@ FEMDeformableBodyModel::initializeDampingMatrix()
if (!m_damped)
{
- return;
+ return false;
}
- if (dampingLaplacianCoefficient <= 0.0)
+ if (dampingLaplacianCoefficient < 0.0)
{
LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix Damping coefficient is negative!";
- return;
+ return false;
}
auto imstkVolMesh = std::static_pointer_cast<VolumetricMesh>(m_forceModelGeometry);
@@ -265,7 +306,7 @@ FEMDeformableBodyModel::initializeDampingMatrix()
if (!meshGraph)
{
LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix: Mesh graph not avaliable!";
- return;
+ return false;
}
vega::SparseMatrix *matrix;
@@ -274,7 +315,7 @@ FEMDeformableBodyModel::initializeDampingMatrix()
if (!matrix)
{
LOG(WARNING) << "DeformableBodyModel::initializeDampingMatrix: Mesh Laplacian not avaliable!";
- return;
+ return false;
}
matrix->ScalarMultiply(dampingLaplacianCoefficient);
@@ -282,15 +323,17 @@ FEMDeformableBodyModel::initializeDampingMatrix()
m_vegaDampingMatrix.reset(matrix);
this->initializeEigenMatrixFromVegaMatrix(*matrix, m_C);
+
+ return true;
}
-void
+bool
FEMDeformableBodyModel::initializeTangentStiffness()
{
if (!m_internalForceModel)
{
LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness: Tangent stiffness cannot be initialized without force model";
- return;
+ return false;
}
vega::SparseMatrix *matrix;
@@ -299,13 +342,13 @@ FEMDeformableBodyModel::initializeTangentStiffness()
if (!matrix)
{
LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness - Tangent stiffness matrix topology not avaliable!";
- return;
+ return false;
}
if (!m_vegaMassMatrix)
{
LOG(WARNING) << "DeformableBodyModel::initializeTangentStiffness - Vega mass matrix doesn't exist!";
- return;
+ return false;
}
matrix->BuildSubMatrixIndices(*m_vegaMassMatrix.get());
@@ -329,9 +372,11 @@ FEMDeformableBodyModel::initializeTangentStiffness()
}
m_internalForceModel->setTangentStiffness(m_vegaTangentStiffnessMatrix);
+
+ return true;
}
-void
+bool
FEMDeformableBodyModel::initializeGravityForce()
{
m_Fgravity.resize(m_numDOF);
@@ -339,6 +384,8 @@ FEMDeformableBodyModel::initializeGravityForce()
const double gravity = m_forceModelConfiguration->getFloatsOptionsMap().at("gravity");
m_vegaPhysicsMesh->computeGravity(m_Fgravity.data(), gravity);
+
+ return true;
}
void
@@ -448,11 +495,13 @@ FEMDeformableBodyModel::computeImplicitSystemLHS(const kinematicState& stateAtT,
}
}
-void
+bool
FEMDeformableBodyModel::initializeExplicitExternalForces()
{
m_FexplicitExternal.resize(m_numDOF);
m_FexplicitExternal.setZero();
+
+ return true;
}
void
diff --git a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
index e0963e86bbf416a2f04c069970752023e87d94e3..b6ac2850ff78859c69095c066732724710fec8e3 100644
--- a/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkFEMDeformableBodyModel.h
@@ -100,7 +100,7 @@ public:
///
/// \brief Initialize the deformable body model
///
- void initialize(std::shared_ptr<VolumetricMesh> physicsMesh);
+ bool initialize() override;
///
/// \brief Load the initial conditions of the deformable object
@@ -110,37 +110,37 @@ public:
///
/// \brief Load the boundary conditions from external file
///
- void loadBoundaryConditions();
+ bool loadBoundaryConditions();
///
/// \brief Initialize the force model
///
- void initializeForceModel();
+ bool initializeForceModel();
///
/// \brief Initialize the mass matrix from the mesh
///
- void initializeMassMatrix(const bool saveToDisk = false);
+ bool initializeMassMatrix(const bool saveToDisk = false);
///
/// \brief Initialize the damping (combines structural and viscous damping) matrix
///
- void initializeDampingMatrix();
+ bool initializeDampingMatrix();
///
/// \brief Initialize the tangent stiffness matrix
///
- void initializeTangentStiffness();
+ bool initializeTangentStiffness();
///
/// \brief Initialize the gravity force
///
- void initializeGravityForce();
+ bool initializeGravityForce();
///
/// \brief Initialize explicit external forces
///
- void initializeExplicitExternalForces();
+ bool initializeExplicitExternalForces();
///
/// \brief Initialize the Eigen matrix with data inside vega sparse matrix
diff --git a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
index 161459f8ce5a7970816003c952b993e0615c50e5..361d7402095aeec59f45eefeca13e8897c408e01 100644
--- a/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
+++ b/Base/DynamicalModels/ObjectModels/imstkPbdModel.cpp
@@ -37,9 +37,6 @@ namespace imstk
PbdModel::PbdModel() :
DynamicalModel(DynamicalModelType::positionBasedDynamics)
{
- m_initialState = std::make_shared<PbdState>();
- m_previousState = std::make_shared<PbdState>();
- m_currentState = std::make_shared<PbdState>();
}
void
@@ -48,11 +45,60 @@ PbdModel::setModelGeometry(std::shared_ptr<PointSet> m)
m_mesh = m;
}
+bool
+PbdModel::configure(const int nCons, ...)
+{
+ va_list args;
+ va_start(args, nCons);
+ for (int i = 0; i < nCons; ++i)
+ {
+ m_constraintConfig.push_back(std::string(va_arg(args, char*)));
+ }
+
+ m_uniformMassValue = va_arg(args, double);
+
+ if (nCons > 0)
+ {
+ char *sgrav = va_arg(args, char*);
+ std::stringstream gStream(sgrav);
+ float x, y, z;
+ gStream >> x;
+ gStream >> y;
+ gStream >> z;
+ Vec3d g(x,y,z);
+ this->setGravity(g);
+
+ this->setTimeStep(va_arg(args, double));
+
+ char *s = va_arg(args, char*);
+ if (strlen(s) > 0)
+ {
+ std::stringstream stream(s);
+ size_t n;
+ while (stream >> n)
+ {
+ m_fixedNodeIds.push_back(n-1);
+ }
+ }
+ this->setMaxNumIterations(va_arg(args, int));
+ }
+ this->setProximity(va_arg(args, double));
+ this->setContactStiffness(va_arg(args, double));
+ this->setNumDegreeOfFreedom(this->getModelGeometry()->getNumVertices() * 3);
+ va_end(args);
+
+ return true;
+}
+
bool
PbdModel::initialize()
{
if (m_mesh)
{
+ m_initialState = std::make_shared<PbdState>();
+ m_previousState = std::make_shared<PbdState>();
+ m_currentState = std::make_shared<PbdState>();
+
bool option[3] = { 1, 0, 0 };
m_initialState->initialize(m_mesh, option);
m_previousState->initialize(m_mesh, option);
@@ -64,19 +110,139 @@ PbdModel::initialize()
m_currentState->setPositions(m_mesh->getVertexPositions());
auto nP = m_mesh->getNumVertices();
- m_invMass.resize(nP, 0);
+
m_mass.resize(nP, 0);
+ this->setUniformMass(m_uniformMassValue);
- return true;
+ m_invMass.resize(nP, 0);
+ for (auto i :m_fixedNodeIds)
+ {
+ setFixedPoint(i);
+ }
}
else
{
LOG(WARNING) << "Model geometry is not yet set! Cannot initialize without model geometry.";
return false;
}
+
+ auto nCons = m_constraintConfig.size();
+ for (int i = 0; i < nCons; ++i)
+ {
+ auto s = m_constraintConfig.at(i).c_str();
+ int len = 0;
+ while (s[len] != ' ' && s[len] != '\0')
+ {
+ ++len;
+ }
+
+ if (strncmp("FEM", &s[0], len) == 0)
+ {
+ int pos = len + 1;
+ len = 0;
+ while (s[pos + len] != ' ' && s[pos + len] != '\0')
+ {
+ ++len;
+ }
+
+ if (strncmp("Corotation", &s[pos], len) == 0)
+ {
+ if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("NeoHookean", &s[pos], len) == 0)
+ {
+ if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("Stvk", &s[pos], len) == 0)
+ {
+ if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK))
+ {
+ return false;
+ }
+ ;
+ }
+ else
+ { // default
+ if (!this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK))
+ {
+ return false;
+ }
+ ;
+ }
+
+ float YoungModulus, PoissonRatio;
+ sscanf(&s[pos + len + 1], "%f %f", &YoungModulus, &PoissonRatio);
+ this->computeLameConstants(YoungModulus, PoissonRatio);
+ }
+ else if (strncmp("Volume", &s[0], len) == 0)
+ {
+ float stiffness;
+ sscanf(&s[len + 1], "%f", &stiffness);
+ if (!this->initializeVolumeConstraints(stiffness))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("Distance", &s[0], len) == 0)
+ {
+ float stiffness;
+ sscanf(&s[len + 1], "%f", &stiffness);
+ if (!this->initializeDistanceConstraints(stiffness))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("Area", &s[0], len) == 0)
+ {
+ float stiffness;
+ sscanf(&s[len + 1], "%f", &stiffness);
+ if (!this->initializeAreaConstraints(stiffness))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("Dihedral", &s[0], len) == 0)
+ {
+ float stiffness;
+ sscanf(&s[len + 1], "%f", &stiffness);
+ if (!this->initializeDihedralConstraints(stiffness))
+ {
+ return false;
+ }
+ ;
+ }
+ else if (strncmp("ConstantDensity", &s[0], len) == 0)
+ {
+ float stiffness;
+ sscanf(&s[len + 1], "%f", &stiffness);
+ if (!this->initializeConstantDensityConstraint(stiffness))
+ {
+ return false;
+ }
+ ;
+ }
+ else
+ {
+ LOG(WARNING) << "PbdModel::initialize() - Type of PBD constraints not identified!";
+ return false;
+ }
+ }
+ return true;
}
-void PbdModel::computeLameConstants(const double& E, const double nu)
+void
+PbdModel::computeLameConstants(const double E, const double nu)
{
m_mu = E/(2*(1+nu));
m_lambda = E*nu/((1-2*nu)*(1+nu));
@@ -108,7 +274,7 @@ PbdModel::initializeFEMConstraints(PbdFEMConstraint::MaterialType type)
}
bool
-PbdModel::initializeVolumeConstraints(const double& stiffness)
+PbdModel::initializeVolumeConstraints(const double stiffness)
{
// Check if constraint type matches the mesh type
if (m_mesh->getType() != Geometry::Type::TetrahedralMesh)
@@ -133,7 +299,7 @@ PbdModel::initializeVolumeConstraints(const double& stiffness)
}
bool
-PbdModel::initializeDistanceConstraints(const double& stiffness)
+PbdModel::initializeDistanceConstraints(const double stiffness)
{
if (m_mesh->getType() == Geometry::Type::TetrahedralMesh)
{
@@ -256,7 +422,7 @@ PbdModel::initializeDistanceConstraints(const double& stiffness)
}
bool
-PbdModel::initializeAreaConstraints(const double& stiffness)
+PbdModel::initializeAreaConstraints(const double stiffness)
{
// check if constraint type matches the mesh type
if (m_mesh->getType() != Geometry::Type::SurfaceMesh)
@@ -281,7 +447,7 @@ PbdModel::initializeAreaConstraints(const double& stiffness)
}
bool
-PbdModel::initializeDihedralConstraints(const double& stiffness)
+PbdModel::initializeDihedralConstraints(const double stiffness)
{
if (m_mesh->getType() != Geometry::Type::SurfaceMesh)
{
@@ -303,7 +469,10 @@ PbdModel::initializeDihedralConstraints(const double& stiffness)
onering[tri[2]].push_back(k);
}
- std::vector<std::vector<bool>> E(triMesh->getNumVertices(), std::vector<bool>(triMesh->getNumVertices(), 1));
+ std::vector<std::vector<bool>> E(triMesh->getNumVertices(),
+ std::vector<bool>(triMesh->getNumVertices(),
+ 1));
+
for (size_t k = 0; k < elements.size(); ++k)
{
auto& tri = elements[k];
@@ -397,11 +566,14 @@ PbdModel::initializeDihedralConstraints(const double& stiffness)
}
bool
-PbdModel::initializeConstantDensityConstraint(const double& stiffness)
+PbdModel::initializeConstantDensityConstraint(const double stiffness)
{
// check if constraint type matches the mesh type
- if (m_mesh->getType() != Geometry::Type::SurfaceMesh && m_mesh->getType() != Geometry::Type::TetrahedralMesh && m_mesh->getType() != Geometry::Type::LineMesh &&
- m_mesh->getType() != Geometry::Type::HexahedralMesh && m_mesh->getType() != Geometry::Type::PointSet)
+ if (m_mesh->getType() != Geometry::Type::SurfaceMesh &&
+ m_mesh->getType() != Geometry::Type::TetrahedralMesh &&
+ m_mesh->getType() != Geometry::Type::LineMesh &&
+ m_mesh->getType() != Geometry::Type::HexahedralMesh &&
+ m_mesh->getType() != Geometry::Type::PointSet)
{
LOG(WARNING) << "Constant constraint should come with a mesh"; //TODO: Really only need a point cloud, so may need to change this.
return false;
@@ -448,12 +620,13 @@ PbdModel::setViscousDamping(const double damping)
}
else
{
- LOG(WARNING) << "WARNING - PbdModel::setViscousDamping: Viscous damping coefficients is out of bounds [0, 1]";
+ LOG(WARNING) << "WARNING - PbdModel::setViscousDamping: " <<
+ "Viscous damping coefficients is out of bounds [0, 1]";
}
}
void
-PbdModel::setUniformMass(const double& val)
+PbdModel::setUniformMass(const double val)
{
if (val != 0.0)
{
@@ -468,7 +641,7 @@ PbdModel::setUniformMass(const double& val)
}
void
-PbdModel::setParticleMass(const double& val, const size_t& idx)
+PbdModel::setParticleMass(const double val, const size_t idx)
{
if (idx < m_mesh->getNumVertices())
{
@@ -478,7 +651,7 @@ PbdModel::setParticleMass(const double& val, const size_t& idx)
}
void
-PbdModel::setFixedPoint(const size_t& idx)
+PbdModel::setFixedPoint(const size_t idx)
{
if (idx < m_mesh->getNumVertices())
{
@@ -487,7 +660,7 @@ PbdModel::setFixedPoint(const size_t& idx)
}
double
-PbdModel::getInvMass(const size_t& idx) const
+PbdModel::getInvMass(const size_t idx) const
{
return m_invMass[idx];
}
diff --git a/Base/DynamicalModels/ObjectModels/imstkPbdModel.h b/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
index 32265183886201852dead09b1b059c957a184079..3982146a39e9b00b3eafaa142d32fe1fe9049a4c 100644
--- a/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
+++ b/Base/DynamicalModels/ObjectModels/imstkPbdModel.h
@@ -51,49 +51,58 @@ public:
///
~PbdModel() = default;
- ///
- /// \brief Initialize the states
- ///
- bool initialize();
-
///
/// \brief Set/Get the geometry (mesh in this case) used by the pbd model
///
void setModelGeometry(std::shared_ptr<PointSet> m);
std::shared_ptr<PointSet> getModelGeometry() const { return m_mesh; }
+ ///
+ /// \brief Configure the PBD model. Arguments should be in the following order
+ /// 1. Number of Constraints (eg: 1)
+ /// 2. Constraint configuration (eg: "FEM NeoHookean 1.0 0.3")
+ /// 3. Mass (eg: 1.0)
+ /// 4. Gravity (eg: "0 -9.8 0")
+ /// 5. TimeStep (eg: 0.001)
+ /// 6. FixedPoint (eg: "10, 21")
+ /// 7. NumberOfIterationInConstraintSolver (eg: 2)
+ /// 8. Proximity (eg: 0.1)
+ /// 9. Contact stiffness (eg: 0.01)
+ ///
+ bool configure(const int nCons, ...);
+
///
/// \brief setElasticModulus
/// \param E Young's modulus
/// \param nu Poisson's ratio
///
- void computeLameConstants(const double& E, const double nu);
+ void computeLameConstants(const double E, const double nu);
///
/// \brief Returns the first Lame constant
///
- const double& getFirstLame() const { return m_mu; }
+ const double getFirstLame() const { return m_mu; }
///
/// \brief Returns the second Lame constant
///
- const double& getSecondLame() const { return m_lambda; }
+ const double getSecondLame() const { return m_lambda; }
///
/// \brief Set the maximum number of iterations for the pbd solver
///
- void setMaxNumIterations(const unsigned int& n) { m_maxIter = n; }
+ void setMaxNumIterations(const unsigned int n) { m_maxIter = n; }
///
/// \brief Get/Set proximity used for collision
///
- void setProximity(const double& prox) { m_proximity = prox; }
+ void setProximity(const double prox) { m_proximity = prox; }
double getProximity() const { return m_proximity; }
///
/// \brief Get/Set contact stiffness that is used for collision constraints
///
- void setContactStiffness(const double& stiffness) { m_contactStiffness = stiffness;}
+ void setContactStiffness(const double stiffness) { m_contactStiffness = stiffness;}
double getContactStiffness() const { return m_contactStiffness; }
///
@@ -104,28 +113,28 @@ public:
///
/// \brief Initialize volume constraints
///
- bool initializeVolumeConstraints(const double& stiffness);
+ bool initializeVolumeConstraints(const double stiffness);
///
/// \brief Initialize distance constraints
///
- bool initializeDistanceConstraints(const double& stiffness);
+ bool initializeDistanceConstraints(const double stiffness);
///
/// \brief Initialize area constraints
///
- bool initializeAreaConstraints(const double& stiffness);
+ bool initializeAreaConstraints(const double stiffness);
///
/// \brief Initialize dihedral constraints
///
- bool initializeDihedralConstraints(const double& stiffness);
+ bool initializeDihedralConstraints(const double stiffness);
///
/// \brief addConstraint add elastic constraint
/// \param constraint
///
- bool initializeConstantDensityConstraint(const double& stiffness);
+ bool initializeConstantDensityConstraint(const double stiffness);
///
/// \todo: add the initialization parameters for the constraint
@@ -135,29 +144,29 @@ public:
inline void addConstraint(std::shared_ptr<PbdConstraint> constraint) { m_constraints.push_back(constraint); }
///
- /// \brief compute delta x and update position
+ /// \brief compute delta x (position) and update position
///
void projectConstraints();
///
- /// \brief Update the model geometry from the newest pbd state
+ /// \brief Update the model geometry from the newest PBD state
///
void updatePhysicsGeometry() override;
///
- /// \brief Update the pbd state from the model geometry
+ /// \brief Update the PBD state from the model geometry
///
void updatePbdStateFromPhysicsGeometry();
///
- /// \brief Returns true if there is atleast one constraint
+ /// \brief Returns true if there is at least one constraint
///
inline bool hasConstraints() const { return !m_constraints.empty(); }
///
/// \brief Set the time step size
///
- void setTimeStep(const double& timeStep) { m_dt = timeStep; };
+ void setTimeStep(const double timeStep) { m_dt = timeStep; };
///
/// \brief Returns the time step size
@@ -178,22 +187,22 @@ public:
///
/// \brief Set uniform mass to all the nodes
///
- void setUniformMass(const double& val);
+ void setUniformMass(const double val);
///
/// \brief Set mass to particular node
///
- void setParticleMass(const double& val, const size_t& idx);
+ void setParticleMass(const double val, const size_t idx);
///
/// \brief Se the node as fixed
///
- void setFixedPoint(const size_t& idx);
+ void setFixedPoint(const size_t idx);
///
/// \brief Get the inverse of mass of a certain node
///
- double getInvMass(const size_t& idx) const;
+ double getInvMass(const size_t idx) const;
///
/// \brief Time integrate the position
@@ -206,31 +215,41 @@ public:
void updateVelocity();
///
- /// \brief
+ /// \brief Update body states given the newest update and the type of update
///
- void updateBodyStates(const Vectord& q, const stateUpdateType updateType = stateUpdateType::displacement) override {};
+ void updateBodyStates(const Vectord& q,
+ const stateUpdateType updateType = stateUpdateType::displacement) override {};
+
+ ///
+ /// \brief Initialize the PBD model
+ ///
+ bool initialize() override;
protected:
- std::shared_ptr<PointSet> m_mesh; ///> PointSet on which the pbd model operates on
- std::vector<std::shared_ptr<PbdConstraint>> m_constraints; ///> List of pbd constraints
+ std::shared_ptr<PointSet> m_mesh; ///> PointSet on which the pbd model operates on
+ std::vector<std::shared_ptr<PbdConstraint>> m_constraints; ///> List of pbd constraints
+
+ std::vector<std::size_t> m_fixedNodeIds; ///> Nodal IDs of the nodes that are fixed
+ std::vector<std::string> m_constraintConfig;
// Lame's constants
- double m_mu; ///> Lame constant
- double m_lambda; ///> Lame constant
+ double m_mu; ///> Lame constant
+ double m_lambda; ///> Lame constant
// Mass properties
- std::vector<double> m_mass; ///> Mass of nodes
- std::vector<double> m_invMass; ///> Inverse of mass of nodes
+ double m_uniformMassValue = 1.0;
+ std::vector<double> m_mass; ///> Mass of nodes
+ std::vector<double> m_invMass; ///> Inverse of mass of nodes
- double m_contactStiffness = 1.; ///> Contact stiffness for collisions
- Vec3d m_gravity; ///> Gravity
+ double m_contactStiffness = 1.; ///> Contact stiffness for collisions
+ Vec3d m_gravity; ///> Gravity
- double m_viscousDampingCoeff = 0.01; ///> Viscous damping coefficient [0, 1]
+ double m_viscousDampingCoeff = 0.01; ///> Viscous damping coefficient [0, 1]
- unsigned int m_maxIter; ///> Max. pbd iterations
- double m_proximity; ///> Proximity for collisions
+ unsigned int m_maxIter; ///> Max. pbd iterations
+ double m_proximity; ///> Proximity for collisions
- double m_dt; ///> Time step size
+ double m_dt; ///> Time step size
};
} // imstk
diff --git a/Base/Scene/imstkScene.cpp b/Base/Scene/imstkScene.cpp
index 7585394fe8626007ac35c09107078f296459cd70..fb798927db53bd75a690942feb705e18df0aecff 100644
--- a/Base/Scene/imstkScene.cpp
+++ b/Base/Scene/imstkScene.cpp
@@ -27,6 +27,21 @@
namespace imstk
{
+bool
+Scene::initialize()
+{
+ for (auto const& it : m_sceneObjectsMap)
+ {
+ auto sceneObject = it.second;
+ if (!sceneObject->initialize())
+ {
+ LOG(WARNING) << "Error initializing scene object: " << sceneObject->getName();
+ return false;
+ }
+ }
+ return true;
+}
+
bool
Scene::isObjectRegistered(std::string sceneObjectName) const
{
diff --git a/Base/Scene/imstkScene.h b/Base/Scene/imstkScene.h
index 426d8e9783a28a7ef72f1b3428ee06af7ce60149..4443c8bdc49a8b4a3a7e450b9a0a35ffa45341b8 100644
--- a/Base/Scene/imstkScene.h
+++ b/Base/Scene/imstkScene.h
@@ -57,6 +57,11 @@ public:
///
~Scene() = default;
+ ///
+ /// \brief Initialize the scene
+ ///
+ bool initialize();
+
///
/// \brief Returns true if the object with a given name is registered, else false
///
diff --git a/Base/SceneElements/Objects/imstkCollidingObject.h b/Base/SceneElements/Objects/imstkCollidingObject.h
index 5a1fed3aafab53769eaf12349bff762edbbae256..1a4c343c45f5bf4ce649eed6e236bbf573ac4322 100644
--- a/Base/SceneElements/Objects/imstkCollidingObject.h
+++ b/Base/SceneElements/Objects/imstkCollidingObject.h
@@ -78,6 +78,21 @@ public:
///
virtual void updateGeometries() override;
+ ///
+ /// \brief Initialize the scene object
+ ///
+ virtual bool initialize()
+ {
+ if (SceneObject::initialize())
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
protected:
std::shared_ptr<Geometry> m_collidingGeometry; ///> Geometry for collisions
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.cpp b/Base/SceneElements/Objects/imstkDeformableObject.cpp
index 0e8ca53ead70d55d34e1261edd3ba8574c7c05cb..0b2ad8125b80e9d9a40178a21d34a6a031fea1ac 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.cpp
+++ b/Base/SceneElements/Objects/imstkDeformableObject.cpp
@@ -26,11 +26,9 @@ namespace imstk
Vectord&
DeformableObject::getContactForce()
{
- //m_defModel = std::dynamic_pointer_cast<imstk::DeformableBodyModel>(m_dynamicalModel);
-
if (!m_defModel)
{
- LOG(WARNING) << "Dynamics pointer cast failure in DeformableObject::getContactForce()";
+ LOG(WARNING) << "deformation model pointer not valid DeformableObject::getContactForce()";
}
return m_defModel->getContactForce();
@@ -40,12 +38,15 @@ bool
DeformableObject::initialize()
{
m_defModel = std::dynamic_pointer_cast<FEMDeformableBodyModel>(m_dynamicalModel);
- if (m_defModel == nullptr)
+ if (m_defModel)
+ {
+ return DynamicObject::initialize();
+ }
+ else
{
- LOG(WARNING) << "Dynamic model set is not of expected type (DeformableBodyModel)!";
+ LOG(WARNING) << "Dynamics pointer cast failure in DeformableObject::initialize()";
return false;
}
- return true;
}
const Vectord&
diff --git a/Base/SceneElements/Objects/imstkDeformableObject.h b/Base/SceneElements/Objects/imstkDeformableObject.h
index 24d3468043c46e39c17ce5dd976a1701ca28fb1f..4797896c6533b9f4ae45d4225d1b8efce7703b13 100644
--- a/Base/SceneElements/Objects/imstkDeformableObject.h
+++ b/Base/SceneElements/Objects/imstkDeformableObject.h
@@ -56,9 +56,9 @@ public:
~DeformableObject() = default;
///
- /// \brief
+ /// \brief Initialize the deformable object
///
- bool initialize();
+ bool initialize() override;
///
/// \brief Initialize the kinematic state of the body
diff --git a/Base/SceneElements/Objects/imstkDynamicObject.h b/Base/SceneElements/Objects/imstkDynamicObject.h
index 6f6aeb561ea28bb0315e0912ae20d4d636514440..e6010bde8dae62c5797beb5466d1f604d41e3863 100644
--- a/Base/SceneElements/Objects/imstkDynamicObject.h
+++ b/Base/SceneElements/Objects/imstkDynamicObject.h
@@ -114,6 +114,21 @@ public:
}
}
+ ///
+ /// \brief Initialize the scene object
+ ///
+ virtual bool initialize() override
+ {
+ if (CollidingObject::initialize())
+ {
+ return m_dynamicalModel->initialize();
+ }
+ else
+ {
+ return false;
+ }
+ }
+
protected:
///
@@ -122,12 +137,12 @@ protected:
DynamicObject(std::string name) : CollidingObject(name){}
std::shared_ptr<DynamicalModel<StateType>> m_dynamicalModel; ///> Dynamical model
- std::shared_ptr<Geometry> m_physicsGeometry; ///> Geometry used for Physics
+ std::shared_ptr<Geometry> m_physicsGeometry; ///> Geometry used for Physics
//Maps
- std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap; ///> Maps from Physics to collision geometry
- std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap; ///> Maps from Physics to visual geometry
- bool m_updateVisualFromPhysicsGeometry = true; ///> Defines if visual is updated from colliding mapping or physics mapping
+ std::shared_ptr<GeometryMap> m_physicsToCollidingGeomMap; ///> Maps from Physics to collision geometry
+ std::shared_ptr<GeometryMap> m_physicsToVisualGeomMap; ///> Maps from Physics to visual geometry
+ bool m_updateVisualFromPhysicsGeometry = true; ///> Defines if visual is updated from colliding mapping or physics mapping
};
} // imstk
diff --git a/Base/SceneElements/Objects/imstkPbdObject.cpp b/Base/SceneElements/Objects/imstkPbdObject.cpp
index fbbdc70d0273ca533b71c69ca9b492711d0b6420..61537799d503cae7a14e12f517073a48aa2f0b33 100644
--- a/Base/SceneElements/Objects/imstkPbdObject.cpp
+++ b/Base/SceneElements/Objects/imstkPbdObject.cpp
@@ -27,155 +27,18 @@
namespace imstk
{
bool
-PbdObject::initialize(int nCons, ...)
+PbdObject::initialize()
{
- std::shared_ptr<PointSet> mesh = std::dynamic_pointer_cast<PointSet>(m_physicsGeometry);
- if (mesh == nullptr)
- {
- LOG(WARNING) << "Physics geometry is not a mesh!";
- return false;
- }
-
m_pbdModel = std::dynamic_pointer_cast<PbdModel>(m_dynamicalModel);
- if (m_pbdModel == nullptr)
- {
- LOG(WARNING) << "Dynamic model set is not of expected type (PbdModel)!";
- return false;
- }
-
- m_pbdModel->setModelGeometry(mesh);
- m_pbdModel->initialize();
-
- va_list args;
- va_start(args, nCons);
- for(int i = 0; i < nCons; ++i)
- {
- char* s = va_arg(args,char*);
- int len = 0;
- while (s[len] != ' ' && s[len] != '\0')
- {
- ++len;
- }
-
- if (strncmp("FEM",&s[0],len) == 0)
- {
- int pos = len+1;
- len = 0;
- while (s[pos+len] != ' ' && s[pos+len] != '\0')
- {
- ++len;
- }
-
- if (strncmp("Corotation",&s[pos],len) == 0)
- {
- LOG(INFO) << "Creating Corotation constraints";
- m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation);
- }
- else if (strncmp("NeoHookean",&s[pos],len) == 0)
- {
- LOG(INFO) << "Creating Neohookean constraints";
- m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean);
- }
- else if (strncmp("Stvk",&s[pos],len) == 0)
- {
- LOG(INFO) << "Creating StVenant-Kirchhoff constraints";
- m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
- }
- else
- { // default
- m_pbdModel->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
- }
-
- float YoungModulus, PoissonRatio;
- sscanf(&s[pos+len+1], "%f %f", &YoungModulus, &PoissonRatio);
- m_pbdModel->computeLameConstants(YoungModulus, PoissonRatio);
- }
- else if (strncmp("Volume",&s[0],len) == 0)
- {
- float stiffness;
- sscanf(&s[len+1], "%f", &stiffness);
- LOG(INFO) << "Creating Volume constraints " << stiffness;
- m_pbdModel->initializeVolumeConstraints(stiffness);
- }
- else if (strncmp("Distance",&s[0],len) == 0)
- {
- float stiffness;
- sscanf(&s[len+1], "%f", &stiffness);
- LOG(INFO) << "Creating Distance constraints " << stiffness;
- m_pbdModel->initializeDistanceConstraints(stiffness);
- }
- else if (strncmp("Area",&s[0],len) == 0)
- {
- float stiffness;
- sscanf(&s[len+1], "%f", &stiffness);
- LOG(INFO) << "Creating Area constraints " << stiffness;
- m_pbdModel->initializeAreaConstraints(stiffness);
- }
- else if (strncmp("Dihedral",&s[0],len) == 0)
- {
- float stiffness;
- sscanf(&s[len+1], "%f", &stiffness);
- LOG(INFO) << "Creating Dihedral constraints " << stiffness;
- m_pbdModel->initializeDihedralConstraints(stiffness);
- }
- else if (strncmp("ConstantDensity", &s[0], len) == 0)
- {
- float stiffness;
- sscanf(&s[len + 1], "%f", &stiffness);
- LOG(INFO) << "Creating Constant Density constraints ";
- m_pbdModel->initializeConstantDensityConstraint(stiffness);
- }
- else
- {
- exit(0);
- }
- }
-
- m_pbdModel->setUniformMass(va_arg(args, double));
-
- if (nCons > 0)
+ if (m_pbdModel)
{
- char* gstring = va_arg(args,char*);
- float x,y,z;
- sscanf(gstring,"%f %f %f", &x, &y, &z);
- Vec3d g(x,y,z);
- m_pbdModel->setGravity(g);
- m_pbdModel->setTimeStep(va_arg(args, double));
- char *s = va_arg(args,char*);
-
- if (strlen(s) > 0)
- {
- while (1)
- {
- int idx = atoi(s);
- m_pbdModel->setFixedPoint(idx - 1);
- while (*s != ' ' && *s != '\0')
- {
- ++s;
- }
-
- if (*s == '\0')
- {
- break;
- }
- else
- {
- ++s;
- }
- }
- }
- m_pbdModel->setMaxNumIterations(va_arg(args,int));
+ return DynamicObject::initialize();
}
-
- if (m_collidingGeometry)
+ else
{
- m_pbdModel->setProximity(va_arg(args,double));
- m_pbdModel->setContactStiffness(va_arg(args,double));
+ LOG(WARNING) << "Dynamics pointer cast failure in PbdObject::initialize()";
+ return false;
}
-
- m_pbdModel->setNumDegreeOfFreedom(m_pbdModel->getModelGeometry()->getNumVertices() * 3);
-
- return true;
}
void
diff --git a/Base/SceneElements/Objects/imstkPbdObject.h b/Base/SceneElements/Objects/imstkPbdObject.h
index b8f06e0d2f76a3771e58e8b98ba5434e2ad6ed1c..554f3fd6bba673ad70fe19c44a68ef0723182be3 100644
--- a/Base/SceneElements/Objects/imstkPbdObject.h
+++ b/Base/SceneElements/Objects/imstkPbdObject.h
@@ -56,10 +56,9 @@ public:
virtual ~PbdObject() = default;
///
- /// \brief Initialize the pbd configuration
- /// TODO: Parse from config file
+ /// \brief Initialize the pbd scene object
///
- bool initialize(int nCons, ...);
+ bool initialize() override;
///
/// \brief Update the position based on Verlet time stepping rule
diff --git a/Base/SceneElements/Objects/imstkSceneObject.h b/Base/SceneElements/Objects/imstkSceneObject.h
index 6b9655d007bd4fb445cacf584d44f4f2e17086ce..4e9b24a1967e390fea28ce79c295770c1ef5ef66 100644
--- a/Base/SceneElements/Objects/imstkSceneObject.h
+++ b/Base/SceneElements/Objects/imstkSceneObject.h
@@ -88,6 +88,11 @@ public:
///
virtual void updateGeometries() {};
+ ///
+ /// \brief Initialize the scene object
+ ///
+ virtual bool initialize(){ return true; };
+
protected:
///
/// \brief Assigns the type of the object
diff --git a/Base/SimulationManager/imstkSimulationManager.cpp b/Base/SimulationManager/imstkSimulationManager.cpp
index bd87151927fd7f9fdff1c6c2c5a3aa64af2beaaf..1eb11613cc84316ff4509a9fdf977a2259c5bee0 100644
--- a/Base/SimulationManager/imstkSimulationManager.cpp
+++ b/Base/SimulationManager/imstkSimulationManager.cpp
@@ -194,7 +194,7 @@ SimulationManager::setCurrentScene(std::string newSceneName, bool unloadCurrentS
return;
}
- std::shared_ptr<Scene> newScene = this->getScene(newSceneName);
+ auto newScene = this->getScene(newSceneName);
if (!newScene)
{
LOG(WARNING) << "Can not find scene";
@@ -253,6 +253,7 @@ SimulationManager::setCurrentScene(std::string newSceneName, bool unloadCurrentS
void
SimulationManager::startSimulation(bool debug)
{
+ this->getCurrentScene()->initialize();
if (m_status != SimulationStatus::INACTIVE)
{
LOG(WARNING) << "Simulation already active";
diff --git a/Examples/Sandbox/main.cpp b/Examples/Sandbox/main.cpp
index 9301c0fb6ca409f1ef45fbd74ed62d80917237d0..ce15e5f9986dfe10a16ecd5bbe07874808374174 100644
--- a/Examples/Sandbox/main.cpp
+++ b/Examples/Sandbox/main.cpp
@@ -1202,7 +1202,7 @@ void testDeformableBody()
dynaModel->configure(iMSTK_DATA_ROOT "/oneTet/oneTet.config");
//dynaModel->configure(iMSTK_DATA_ROOT"/asianDragon/asianDragon.config");
//dynaModel->configure(iMSTK_DATA_ROOT"/liver/liver.config");
- dynaModel->initialize(volTetMesh);
+ dynaModel->setModelGeometry(volTetMesh);
auto timeIntegrator = std::make_shared<BackwardEuler>(0.001);// Create and add Backward Euler time integrator
dynaModel->setTimeIntegrator(timeIntegrator);
@@ -1349,13 +1349,8 @@ void testPbdVolume()
auto deformableObj = std::make_shared<PbdObject>("Beam");
auto pbdModel = std::make_shared<PbdModel>();
-
- deformableObj->setDynamicalModel(pbdModel);
- deformableObj->setVisualGeometry(surfMesh);
- deformableObj->setPhysicsGeometry(volTetMesh);
- deformableObj->setPhysicsToVisualMap(oneToOneNodalMap); //assign the computed map
-
- deformableObj->initialize(/*Number of Constraints*/ 1,
+ pbdModel->setModelGeometry(volTetMesh);
+ pbdModel->configure(/*Number of Constraints*/ 1,
/*Constraint configuration*/ "FEM NeoHookean 100.0 0.3",
/*Mass*/ 1.0,
/*Gravity*/ "0 -9.8 0",
@@ -1364,13 +1359,17 @@ void testPbdVolume()
/*NumberOfIterationInConstraintSolver*/ 5
);
+ deformableObj->setDynamicalModel(pbdModel);
+ deformableObj->setVisualGeometry(surfMesh);
+ deformableObj->setPhysicsGeometry(volTetMesh);
+ deformableObj->setPhysicsToVisualMap(oneToOneNodalMap); //assign the computed map
+
auto pbdSolver = std::make_shared<PbdSolver>();
pbdSolver->setPbdObject(deformableObj);
scene->addNonlinearSolver(pbdSolver);
scene->addSceneObject(deformableObj);
-
auto planeGeom = std::make_shared<Plane>();
planeGeom->setWidth(40);
planeGeom->setTranslation(0, -6, 0);
@@ -1435,21 +1434,21 @@ void testPbdCloth()
surfMesh->setTrianglesVertices(triangles);
- // Object & Model
+ // Create Object & Model
auto deformableObj = std::make_shared<PbdObject>("Cloth");
auto pbdModel = std::make_shared<PbdModel>();
- deformableObj->setDynamicalModel(pbdModel);
- deformableObj->setVisualGeometry(surfMesh);
- deformableObj->setPhysicsGeometry(surfMesh);
- deformableObj->initialize(/*Number of constraints*/ 2,
+ pbdModel->setModelGeometry(surfMesh);
+ pbdModel->configure(/*Number of constraints*/ 2,
/*Constraint configuration*/ "Distance 0.1",
/*Constraint configuration*/ "Dihedral 0.001",
/*Mass*/ 1.0,
/*Gravity*/ "0 -9.8 0",
/*TimeStep*/ 0.01,
/*FixedPoint*/ "1 2 3 4 5 6 7 8 9 10 11",
- /*NumberOfIterationInConstraintSolver*/ 5
- );
+ /*NumberOfIterationInConstraintSolver*/ 5);
+ deformableObj->setDynamicalModel(pbdModel);
+ deformableObj->setVisualGeometry(surfMesh);
+ deformableObj->setPhysicsGeometry(surfMesh);
// Solver
auto pbdSolver = std::make_shared<PbdSolver>();
@@ -1552,9 +1551,7 @@ void testPbdCollision()
deformableObj->setCollidingToVisualMap(deformMapC2V);
auto pbdModel = std::make_shared<PbdModel>();
- deformableObj->setDynamicalModel(pbdModel);
-
- deformableObj->initialize(/*Number of Constraints*/ 1,
+ pbdModel->configure(/*Number of Constraints*/ 1,
/*Constraint configuration*/ "FEM NeoHookean 1.0 0.3",
/*Mass*/ 1.0,
/*Gravity*/ "0 -9.8 0",
@@ -1563,6 +1560,7 @@ void testPbdCollision()
/*NumberOfIterationInConstraintSolver*/ 2,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 0.01);
+ deformableObj->setDynamicalModel(pbdModel);
auto pbdSolver = std::make_shared<PbdSolver>();
pbdSolver->setPbdObject(deformableObj);
@@ -1630,16 +1628,16 @@ void testPbdCollision()
floor->setPhysicsToCollidingMap(oneToOneFloor);
floor->setPhysicsToVisualMap(oneToOneFloor);
//floor->setCollidingToVisualMap(oneToOneFloor);
- floor->initialize(/*Number of constraints*/ 2,
- /*Constraint configuration*/ "Distance 0.1",
- /*Constraint configuration*/ "Dihedral 0.001",
- /*Mass*/ 0.1,
- /*Gravity*/ "0 9.8 0",
- /*TimeStep*/ 0.002,
- /*FixedPoint*/ fixed_corner.c_str(),
- /*NumberOfIterationInConstraintSolver*/ 5,
- /*Proximity*/ 0.1,
- /*Contact stiffness*/ 0.95);
+ //floor->initialize(/*Number of constraints*/ 2,
+ // /*Constraint configuration*/ "Distance 0.1",
+ // /*Constraint configuration*/ "Dihedral 0.001",
+ // /*Mass*/ 0.1,
+ // /*Gravity*/ "0 9.8 0",
+ // /*TimeStep*/ 0.002,
+ // /*FixedPoint*/ fixed_corner.c_str(),
+ // /*NumberOfIterationInConstraintSolver*/ 5,
+ // /*Proximity*/ 0.1,
+ // /*Contact stiffness*/ 0.95);
scene->addSceneObject(floor);
std::cout << "nbr of vertices in cloth mesh" << clothMesh->getNumVertices() << std::endl;
@@ -1707,15 +1705,15 @@ void testPbdCollision()
deformableObj1->setPhysicsToCollidingMap(deformMapP2C1);
deformableObj1->setPhysicsToVisualMap(deformMapP2V1);
deformableObj1->setCollidingToVisualMap(deformMapC2V1);
- deformableObj1->initialize(/*Number of Constraints*/ 1,
- /*Constraint configuration*/ "FEM NeoHookean 10.0 0.5",
- /*Mass*/ 0.0,
- /*Gravity*/ "0 -9.8 0",
- /*TimeStep*/ 0.002,
- /*FixedPoint*/ "",
- /*NumberOfIterationInConstraintSolver*/ 2,
- /*Proximity*/ 0.1,
- /*Contact stiffness*/ 0.01);
+ //deformableObj1->initialize(/*Number of Constraints*/ 1,
+ // /*Constraint configuration*/ "FEM NeoHookean 10.0 0.5",
+ // /*Mass*/ 0.0,
+ // /*Gravity*/ "0 -9.8 0",
+ // /*TimeStep*/ 0.002,
+ // /*FixedPoint*/ "",
+ // /*NumberOfIterationInConstraintSolver*/ 2,
+ // /*Proximity*/ 0.1,
+ // /*Contact stiffness*/ 0.01);
scene->addSceneObject(deformableObj1);
@@ -1794,12 +1792,11 @@ void testPbdCollision()
floor->setCollidingToVisualMap(floorMapC2V);
auto pbdModel2 = std::make_shared<PbdModel>();
- floor->setDynamicalModel(pbdModel2);
-
- floor->initialize(/*Number of Constraints*/ 0,
+ pbdModel2->configure(/*Number of Constraints*/ 0,
/*Mass*/ 0.0,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 1.0);
+ floor->setDynamicalModel(pbdModel2);
auto pbdSolverfloor = std::make_shared<PbdSolver>();
pbdSolverfloor->setPbdObject(floor);
@@ -1865,9 +1862,8 @@ void testPbdFluidBenchmarking()
cube->setPhysicsGeometry(cubeMesh);
auto pbdModel = std::make_shared<PbdModel>();
- cube->setDynamicalModel(pbdModel);
-
- cube->initialize(/*Number of Constraints*/ 1,
+ pbdModel->setModelGeometry(cubeMesh);
+ pbdModel->configure(/*Number of Constraints*/ 1,
/*Constraint configuration*/ "ConstantDensity 1.0 0.3",
/*Mass*/ 1.0,
/*Gravity*/ "0 -9.8 0",
@@ -1876,6 +1872,7 @@ void testPbdFluidBenchmarking()
/*NumberOfIterationInConstraintSolver*/ 2,
/*Proximity*/ 0.2,
/*Contact stiffness*/ 1.0);
+ cube->setDynamicalModel(pbdModel);
auto pbdSolver = std::make_shared<PbdSolver>();
pbdSolver->setPbdObject(cube);
@@ -1946,12 +1943,12 @@ void testPbdFluidBenchmarking()
floor->setCollidingToVisualMap(floorMapC2V);
auto pbdModel2 = std::make_shared<PbdModel>();
- floor->setDynamicalModel(pbdModel2);
-
- floor->initialize(/*Number of Constraints*/ 0,
+ pbdModel2->setModelGeometry(floorMeshPhysics);
+ pbdModel2->configure(/*Number of Constraints*/ 0,
/*Mass*/ 0.0,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 1.0);
+ floor->setDynamicalModel(pbdModel2);
auto pbdSolverfloor = std::make_shared<PbdSolver>();
pbdSolverfloor->setPbdObject(floor);
@@ -2031,9 +2028,8 @@ void testPbdFluid()
deformableObj->setPhysicsGeometry(fluidMesh);
auto pbdModel = std::make_shared<PbdModel>();
- deformableObj->setDynamicalModel(pbdModel);
-
- deformableObj->initialize(/*Number of Constraints*/ 1,
+ pbdModel->setModelGeometry(fluidMesh);
+ pbdModel->configure(/*Number of Constraints*/ 1,
/*Constraint configuration*/ "ConstantDensity 1.0 0.3",
/*Mass*/ 1.0,
/*Gravity*/ "0 -9.8 0",
@@ -2042,6 +2038,7 @@ void testPbdFluid()
/*NumberOfIterationInConstraintSolver*/ 2,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 1.0);
+ deformableObj->setDynamicalModel(pbdModel);
auto pbdSolver = std::make_shared<PbdSolver>();
pbdSolver->setPbdObject(deformableObj);
@@ -2187,12 +2184,12 @@ void testPbdFluid()
floor->setCollidingToVisualMap(floorMapC2V);
auto pbdModel2 = std::make_shared<PbdModel>();
- floor->setDynamicalModel(pbdModel2);
-
- floor->initialize(/*Number of Constraints*/ 0,
+ pbdModel2->setModelGeometry(floorMeshPhysics);
+ pbdModel2->configure(/*Number of Constraints*/ 0,
/*Mass*/ 0.0,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 1.0);
+ floor->setDynamicalModel(pbdModel2);
auto pbdSolverfloor = std::make_shared<PbdSolver>();
pbdSolverfloor->setPbdObject(floor);
@@ -2246,11 +2243,7 @@ void testLineMesh()
if (line)
{
- // Read LineMesh
- auto lineMeshColliding = std::make_shared<imstk::LineMesh>();
- auto lineMeshVisual = std::make_shared<imstk::LineMesh>();
- auto lineMeshPhysics = std::make_shared<imstk::LineMesh>();
-
+ // Make LineMesh
imstk::StdVectorOfVec3d vertList;
vertList.resize(3);
vertList[0] = Vec3d(0.0, -10.0, -10.0);
@@ -2267,101 +2260,43 @@ void testLineMesh()
connectivity.push_back(line);
}
- lineMeshColliding->setInitialVertexPositions(vertList);
- lineMeshColliding->setVertexPositions(vertList);
- lineMeshColliding->setConnectivity(connectivity);
-
- lineMeshPhysics->setInitialVertexPositions(vertList);
- lineMeshPhysics->setVertexPositions(vertList);
- lineMeshPhysics->setConnectivity(connectivity);
-
- lineMeshVisual->setInitialVertexPositions(vertList);
- lineMeshVisual->setVertexPositions(vertList);
- lineMeshVisual->setConnectivity(connectivity);
-
-
- auto mapC2P = std::make_shared<imstk::OneToOneMap>();
- mapC2P->setMaster(lineMeshColliding);
- mapC2P->setSlave(lineMeshPhysics);
- mapC2P->compute();
-
- auto mapC2V = std::make_shared<imstk::OneToOneMap>();
- mapC2V->setMaster(lineMeshColliding);
- mapC2V->setSlave(lineMeshVisual);
- mapC2V->compute();
-
- auto mapP2C = std::make_shared<imstk::OneToOneMap>();
- mapP2C->setMaster(lineMeshPhysics);
- mapP2C->setSlave(lineMeshColliding);
- mapP2C->compute();
-
- auto mapP2V = std::make_shared<imstk::OneToOneMap>();
- mapP2V->setMaster(lineMeshPhysics);
- mapP2V->setSlave(lineMeshVisual);
- mapP2V->compute();
-
- linesTool->setCollidingGeometry(lineMeshColliding);
- linesTool->setVisualGeometry(lineMeshVisual);
- linesTool->setPhysicsGeometry(lineMeshPhysics);
- linesTool->setPhysicsToCollidingMap(mapP2C);
- linesTool->setCollidingToVisualMap(mapC2V);
- linesTool->setPhysicsToVisualMap(mapP2V);
- linesTool->setColldingToPhysicsMap(mapC2P);
- linesTool->initialize(/*Number of constraints*/ 1,
- /*Constraint configuration*/ "Distance 100",
- /*Mass*/ 0.0,
- /*Gravity*/ "0 -9.8 0",
- /*TimeStep*/ 0.002,
- /*FixedPoint*/ "0 1 2",
- /*NumberOfIterationInConstraintSolver*/ 5,
- /*Proximity*/ 0.1,
- /*Contact stiffness*/ 0.1);
+ auto lineMesh = std::make_shared<imstk::LineMesh>();
+ lineMesh->setInitialVertexPositions(vertList);
+ lineMesh->setVertexPositions(vertList);
+ lineMesh->setConnectivity(connectivity);
+
+ linesTool->setCollidingGeometry(lineMesh);
+ linesTool->setVisualGeometry(lineMesh);
+ linesTool->setPhysicsGeometry(lineMesh);
+ //linesTool->initialize(/*Number of constraints*/ 1,
+ // /*Constraint configuration*/ "Distance 100",
+ // /*Mass*/ 0.0,
+ // /*Gravity*/ "0 -9.8 0",
+ // /*TimeStep*/ 0.002,
+ // /*FixedPoint*/ "0 1 2",
+ // /*NumberOfIterationInConstraintSolver*/ 5,
+ // /*Proximity*/ 0.1,
+ // /*Contact stiffness*/ 0.1);
scene->addSceneObject(linesTool);
}
else
{
std::string path2obj = iMSTK_DATA_ROOT "/ETI/resources/Tools/blade2.obj";
- auto collidingMesh = imstk::MeshIO::read(path2obj);
- auto viusalMesh = imstk::MeshIO::read(path2obj);
- auto physicsMesh = imstk::MeshIO::read(path2obj);
-
- auto bladeMapP2V = std::make_shared<imstk::OneToOneMap>();
- bladeMapP2V->setMaster(physicsMesh);
- bladeMapP2V->setSlave(viusalMesh);
- bladeMapP2V->compute();
-
- auto bladeMapP2C = std::make_shared<imstk::OneToOneMap>();
- bladeMapP2C->setMaster(physicsMesh);
- bladeMapP2C->setSlave(collidingMesh);
- bladeMapP2C->compute();
-
- auto bladeMapC2V = std::make_shared<imstk::OneToOneMap>();
- bladeMapC2V->setMaster(collidingMesh);
- bladeMapC2V->setSlave(viusalMesh);
- bladeMapC2V->compute();
-
- auto bladeMapC2P = std::make_shared<imstk::OneToOneMap>();
- bladeMapC2P->setMaster(collidingMesh);
- bladeMapC2P->setSlave(physicsMesh);
- bladeMapC2P->compute();
-
- blade->setCollidingGeometry(collidingMesh);
- blade->setVisualGeometry(viusalMesh);
- blade->setPhysicsGeometry(physicsMesh);
- blade->setPhysicsToCollidingMap(bladeMapP2C);
- blade->setCollidingToVisualMap(bladeMapC2V);
- blade->setPhysicsToVisualMap(bladeMapP2V);
- blade->setColldingToPhysicsMap(bladeMapC2P);
- blade->initialize(/*Number of constraints*/ 1,
- /*Constraint configuration*/ "Distance 0.1",
- /*Mass*/ 0.0,
- /*Gravity*/ "0 0 0",
- /*TimeStep*/ 0.001,
- /*FixedPoint*/ "",
- /*NumberOfIterationInConstraintSolver*/ 5,
- /*Proximity*/ 0.1,
- /*Contact stiffness*/ 0.01);
+ auto bladeMesh = imstk::MeshIO::read(path2obj);
+
+ blade->setCollidingGeometry(bladeMesh);
+ blade->setVisualGeometry(bladeMesh);
+ blade->setPhysicsGeometry(bladeMesh);
+ //blade->initialize(/*Number of constraints*/ 1,
+ // /*Constraint configuration*/ "Distance 0.1",
+ // /*Mass*/ 0.0,
+ // /*Gravity*/ "0 0 0",
+ // /*TimeStep*/ 0.001,
+ // /*FixedPoint*/ "",
+ // /*NumberOfIterationInConstraintSolver*/ 5,
+ // /*Proximity*/ 0.1,
+ // /*Contact stiffness*/ 0.01);
scene->addSceneObject(blade);
}
@@ -2410,35 +2345,13 @@ void testLineMesh()
auto clothMeshVisual = std::make_shared<imstk::SurfaceMesh>();
clothMeshVisual->initialize(vertList, triangles);
- auto clothMeshColliding = std::make_shared<imstk::SurfaceMesh>();
- clothMeshColliding->initialize(vertList, triangles);
- auto clothMeshPhysics = std::make_shared<imstk::SurfaceMesh>();
- clothMeshPhysics->initialize(vertList, triangles);
-
- auto clothMapP2V = std::make_shared<imstk::OneToOneMap>();
- clothMapP2V->setMaster(clothMeshPhysics);
- clothMapP2V->setSlave(clothMeshVisual);
- clothMapP2V->compute();
-
- auto clothMapC2V = std::make_shared<imstk::OneToOneMap>();
- clothMapC2V->setMaster(clothMeshColliding);
- clothMapC2V->setSlave(clothMeshVisual);
- clothMapC2V->compute();
-
- auto clothMapP2C = std::make_shared<imstk::OneToOneMap>();
- clothMapP2C->setMaster(clothMeshPhysics);
- clothMapP2C->setSlave(clothMeshColliding);
- clothMapP2C->compute();
-
-
- auto floor = std::make_shared<PbdObject>("cloth");
- floor->setCollidingGeometry(clothMeshColliding);
- floor->setVisualGeometry(clothMeshVisual);
- floor->setPhysicsGeometry(clothMeshPhysics);
- floor->setPhysicsToCollidingMap(clothMapP2C);
- floor->setPhysicsToVisualMap(clothMapP2V);
- floor->setCollidingToVisualMap(clothMapC2V);
- floor->initialize(/*Number of constraints*/ 2,
+
+ auto clothObject = std::make_shared<PbdObject>("cloth");
+ clothObject->setCollidingGeometry(clothMeshVisual);
+ clothObject->setVisualGeometry(clothMeshVisual);
+ clothObject->setPhysicsGeometry(clothMeshVisual);
+ auto clothModel = std::make_shared<PbdModel>();
+ clothModel->configure(/*Number of constraints*/ 2,
/*Constraint configuration*/ "Distance 0.1",
/*Constraint configuration*/ "Dihedral 0.001",
/*Mass*/ 0.1,
@@ -2448,7 +2361,8 @@ void testLineMesh()
/*NumberOfIterationInConstraintSolver*/ 5,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 0.1);
- scene->addSceneObject(floor);
+ clothObject->setDynamicalModel(clothModel);
+ scene->addSceneObject(clothObject);
std::cout << "nbr of vertices in cloth mesh" << clothMeshVisual->getNumVertices() << std::endl;
@@ -2463,7 +2377,7 @@ void testLineMesh()
tool = blade;
}
- auto pair1 = std::make_shared<PbdInteractionPair>(PbdInteractionPair(tool, floor));
+ auto pair1 = std::make_shared<PbdInteractionPair>(PbdInteractionPair(tool, clothObject));
pair1->setNumberOfInterations(5);
clothTestcolGraph->addInteractionPair(pair1);
@@ -2526,7 +2440,9 @@ void testLineMesh()
deformableObj->setPhysicsToCollidingMap(dragonMapP2C);
deformableObj->setPhysicsToVisualMap(dragonMapP2V);
deformableObj->setCollidingToVisualMap(dragonMapC2V);
- deformableObj->initialize(/*Number of Constraints*/ 1,
+
+ auto pbdDefModel = std::make_shared<PbdModel>();
+ pbdDefModel->configure(/*Number of Constraints*/ 1,
/*Constraint configuration*/ "FEM NeoHookean 10.0 0.3",
/*Mass*/ 0.1,
/*Gravity*/ "0 0 0",
@@ -2535,7 +2451,7 @@ void testLineMesh()
/*NumberOfIterationInConstraintSolver*/ 5,
/*Proximity*/ 0.1,
/*Contact stiffness*/ 0.01);
-
+ deformableObj->setDynamicalModel(pbdDefModel);
scene->addSceneObject(deformableObj);
std::cout << "nbr of vertices in tongue mesh = " << surfMesh->getNumVertices() << std::endl;
@@ -2674,7 +2590,7 @@ void testDeformableBodyCollision()
// Configure the dynamic model
auto dynaModel = std::make_shared<FEMDeformableBodyModel>();
dynaModel->configure(iMSTK_DATA_ROOT "/oneTet/oneTet.config");
- dynaModel->initialize(volTetMesh);
+ dynaModel->setModelGeometry(volTetMesh);
// Create and add Backward Euler time integrator
auto timeIntegrator = std::make_shared<BackwardEuler>(0.001);
@@ -2778,7 +2694,7 @@ void liverToolInteraction()
auto dynaModel = std::make_shared<FEMDeformableBodyModel>();
//dynaModel->configure(iMSTK_DATA_ROOT"/liver/liver.config");
dynaModel->configure(iMSTK_DATA_ROOT "/oneTet/oneTet.config");
- dynaModel->initialize(volTetMesh);
+ dynaModel->setModelGeometry(volTetMesh);
// Create and add Backward Euler time integrator
auto timeIntegrator = std::make_shared<BackwardEuler>(0.001);
@@ -3067,7 +2983,7 @@ void testPicking()
// Configure the dynamic model
auto dynaModel = std::make_shared<FEMDeformableBodyModel>();
dynaModel->configure(iMSTK_DATA_ROOT "/oneTet/oneTet.config");
- dynaModel->initialize(volTetMesh);
+ dynaModel->setModelGeometry(volTetMesh);
// Create and add Backward Euler time integrator
auto timeIntegrator = std::make_shared<BackwardEuler>(0.01);
@@ -3347,7 +3263,7 @@ int main()
/*------------------
Test physics
------------------*/
- //testPbdVolume();
+ testPbdVolume();
//testPbdCloth();
//testPbdCollision();
testPbdFluidBenchmarking();