Commit 5f88f46c authored by Sreekanth Arikatla's avatar Sreekanth Arikatla
Browse files

ENH: Add configure function to dynamic models

parent a7158781
Pipeline #66588 passed with stage
......@@ -119,6 +119,11 @@ public:
///
virtual void updatePhysicsGeometry(){};
///
/// \brief Initialize the dynamical model
///
virtual bool initialize() = 0;
protected:
DynamicalModelType m_type; ///> Mathematical model type
......
......@@ -89,16 +89,16 @@ 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);
auto physicsMesh = std::dynamic_pointer_cast<imstk::VolumetricMesh>(this->getModelGeometry());
// 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;
}
m_vegaPhysicsMesh = physicsMesh->getAttachedVegaMesh();
......@@ -119,6 +119,8 @@ 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
......
......@@ -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
......
......@@ -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,16 +110,106 @@ 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)
{
LOG(INFO) << "Creating Corotation constraints";
this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::Corotation);
}
else if (strncmp("NeoHookean", &s[pos], len) == 0)
{
LOG(INFO) << "Creating Neohookean constraints";
this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::NeoHookean);
}
else if (strncmp("Stvk", &s[pos], len) == 0)
{
LOG(INFO) << "Creating StVenant-Kirchhoff constraints";
this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
}
else
{ // default
this->initializeFEMConstraints(PbdFEMConstraint::MaterialType::StVK);
}
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);
LOG(INFO) << "Creating Volume constraints " << stiffness;
this->initializeVolumeConstraints(stiffness);
}
else if (strncmp("Distance", &s[0], len) == 0)
{
float stiffness;
sscanf(&s[len + 1], "%f", &stiffness);
LOG(INFO) << "Creating Distance constraints " << stiffness;
this->initializeDistanceConstraints(stiffness);
}
else if (strncmp("Area", &s[0], len) == 0)
{
float stiffness;
sscanf(&s[len + 1], "%f", &stiffness);
LOG(INFO) << "Creating Area constraints " << stiffness;
this->initializeAreaConstraints(stiffness);
}
else if (strncmp("Dihedral", &s[0], len) == 0)
{
float stiffness;
sscanf(&s[len + 1], "%f", &stiffness);
LOG(INFO) << "Creating Dihedral constraints " << stiffness;
this->initializeDihedralConstraints(stiffness);
}
else if (strncmp("ConstantDensity", &s[0], len) == 0)
{
float stiffness;
sscanf(&s[len + 1], "%f", &stiffness);
LOG(INFO) << "Creating Constant Density constraints ";
this->initializeConstantDensityConstraint(stiffness);
}
else
{
exit(0);
}
}
return true;
}
void PbdModel::computeLameConstants(const double& E, const double nu)
......
......@@ -51,17 +51,26 @@ 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
......@@ -210,17 +219,26 @@ public:
///
void updateBodyStates(const Vectord& q, const stateUpdateType updateType = stateUpdateType::displacement) override {};
///
/// \brief Initialize the scene object
///
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::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
// 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
......
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment