Commit 49bb04d0 authored by Sreekanth Arikatla's avatar Sreekanth Arikatla
Browse files

Merge branch 'DeformableBodybugFixes' into 'master'

Bug fixes and refactoring

See merge request !292
parents fd9f91f2 837d2d03
Pipeline #134626 failed with stage
in 0 seconds
......@@ -32,6 +32,7 @@
#include "imstkMeshIO.h"
#include "imstkOneToOneMap.h"
#include "imstkAPIUtilities.h"
#include "imstkConjugateGradient.h"
using namespace imstk;
......@@ -120,24 +121,21 @@ int main()
dynaModel->getFunction(),
dynaModel->getFunctionGradient());
std::vector<LinearProjectionConstraint> projList;
for (auto i : dynaModel->getFixNodeIds())
{
auto s = LinearProjectionConstraint(i, false);
s.setProjectorToDirichlet(i, Vec3d(0.001, 0, 0));
projList.push_back(s);
}
nlSystem->setUnknownVector(dynaModel->getUnknownVec());
nlSystem->setUpdateFunction(dynaModel->getUpdateFunction());
nlSystem->setUpdatePreviousStatesFunction(dynaModel->getUpdatePrevStateFunction());
// create a linear solver
auto linSolver = std::make_shared<GaussSeidel>();
auto linSolver = std::make_shared<ConjugateGradient>();
if (linSolver->getType() == imstk::LinearSolver<imstk::SparseMatrixd>::Type::GaussSeidel &&
dynaModel->isFixedBCImplemented())
{
LOG(WARNING) << "The GS solver may not be viable!";
}
// create a non-linear solver and add to the scene
auto nlSolver = std::make_shared<NewtonSolver>();
linSolver->setLinearProjectors(&projList);
nlSolver->setLinearSolver(linSolver);
nlSolver->setSystem(nlSystem);
scene->addNonlinearSolver(nlSolver);
......@@ -154,7 +152,8 @@ int main()
// Run the simulation
sdk->setActiveScene(scene);
sdk->startSimulation(SimulationStatus::PAUSED);
sdk->getViewer()->setBackgroundColors(Vec3d(0.3285, 0.3285, 0.6525), Vec3d(0.13836, 0.13836, 0.2748), true);
sdk->startSimulation();
return 0;
}
......@@ -32,7 +32,7 @@ using namespace imstk;
int main()
{
// SDK and Scene
auto sdk = std::make_shared<SimulationManager>(false, true);
auto sdk = std::make_shared<SimulationManager>();
auto scene = sdk->createNewScene("Rendering");
// Add IBL Probe
......
......@@ -31,7 +31,7 @@ using namespace imstk;
int main()
{
// SDK and Scene
auto sdk = std::make_shared<SimulationManager>(1);
auto sdk = std::make_shared<SimulationManager>();
auto sceneTest = sdk->createNewScene("HelloWorld");
// Add Plane
......
......@@ -20,6 +20,7 @@
=========================================================================*/
#include "imstkForceModelConfig.h"
#include <iostream>
namespace imstk
{
......@@ -70,10 +71,18 @@ ForceModelConfig::parseConfig(const std::string &configFileName)
vegaConfigFileOptions.printOptions();
}
// get the root directory
std::string rootDir;
const size_t last_slash_idx = configFileName.rfind('/');
if (std::string::npos != last_slash_idx)
{
rootDir = configFileName.substr(0, last_slash_idx);
}
// Store parsed string values
m_stringsOptionMap.emplace(optNameList.femMethodName, optList.femMethod);
m_stringsOptionMap.emplace(optNameList.invertibleMaterialName, optList.invertibleMaterial);
m_stringsOptionMap.emplace(optNameList.fixedDOFFilenameName, optList.fixedDOFFilename);
m_stringsOptionMap.emplace(optNameList.fixedDOFFilenameName, rootDir + std::string("/") + optList.fixedDOFFilename);
// Store parsed floating point values
m_floatsOptionMap.emplace(optNameList.dampingMassCoefficientName, optList.dampingMassCoefficient);
......
......@@ -514,20 +514,23 @@ FEMDeformableBodyModel::applyBoundaryConditions(SparseMatrixd &M, const bool wit
// Set column and row to zero.
for (auto & index : m_fixedNodeIds)
{
auto idx = static_cast<SparseMatrixd::Index>(index);
auto nodeIdx = static_cast<SparseMatrixd::Index>(index)*3;
for (int k = 0; k < M.outerSize(); ++k)
for (auto idx = nodeIdx; idx < nodeIdx + 3; idx++)
{
for (SparseMatrixd::InnerIterator i(M, k); i; ++i)
for (long long k = 0; k < M.outerSize(); ++k)
{
if (i.row() == idx || i.col() == idx)
for (SparseMatrixd::InnerIterator i(M, k); i; ++i)
{
i.valueRef() = 0.0;
}
if (i.row() == idx && i.col() == idx)
{
i.valueRef() = compliance;
if (i.row() == idx || i.col() == idx)
{
i.valueRef() = 0.0;
}
if (i.row() == idx && i.col() == idx)
{
i.valueRef() = compliance;
}
}
}
}
......@@ -539,7 +542,8 @@ FEMDeformableBodyModel::applyBoundaryConditions(Vectord &x) const
{
for (auto & index : m_fixedNodeIds)
{
x(index) = 0.0;
const auto _3Index = 3 * index;
x(_3Index) = x(_3Index + 1) = x(_3Index + 2) = 0.0;
}
}
......@@ -610,7 +614,10 @@ FEMDeformableBodyModel::getFunction()
(semiImplicit) ?
this->computeSemiImplicitSystemRHS(*m_previousState.get(), *m_currentState.get(), m_updateType) :
this->computeImplicitSystemRHS(*m_previousState.get(), *m_currentState.get(), m_updateType);
if (this->m_implementFixedBC)
{
applyBoundaryConditions(m_Feff);
}
return m_Feff;
};
}
......@@ -623,6 +630,10 @@ FEMDeformableBodyModel::getFunctionGradient()
{
this->computeImplicitSystemLHS(*m_previousState.get(), *m_currentState.get(), m_updateType);
if(this->m_implementFixedBC)
{
applyBoundaryConditions(m_Keff);
}
return m_Keff;
};
}
......
......@@ -248,6 +248,13 @@ public:
///
void setFixedSizeTimeStepping();
///
/// \brief Set the fixed BC implementation state
///
void enableFixedBC() { m_implementFixedBC = true; };
void disableFixedBC() { m_implementFixedBC = false; };
bool isFixedBCImplemented() const { return m_implementFixedBC; };
protected:
std::shared_ptr<InternalForceModel> m_internalForceModel; ///> Mathematical model for intenal forces
std::shared_ptr<TimeIntegrator> m_timeIntegrator; ///> Time integrator
......@@ -278,6 +285,10 @@ protected:
stateUpdateType m_updateType = stateUpdateType::deltaVelocity; ///> Update type of the model
bool m_damped = false; ///> Viscous or structurally damped system
// If this is true, the tangent stiffness and force vector will be modified to
// accommodate (the rows and columns will be nullified) the fixed boundary conditions
bool m_implementFixedBC = true;
};
} // imstk
......
......@@ -172,8 +172,13 @@ VTKRenderer::VTKRenderer(std::shared_ptr<Scene> scene, const bool enableVR)
#endif
///TODO : based on scene properties
// Customize background colors
m_vtkRenderer->SetBackground(0.66,0.66,0.66);
m_vtkRenderer->SetBackground2(0.406, 0.481, 0.497);
//m_vtkRenderer->SetBackground(0.66,0.66,0.66);
//m_vtkRenderer->SetBackground2(0.406, 0.481, 0.497);
m_vtkRenderer->SetBackground(0.3285, 0.3285, 0.6525);
m_vtkRenderer->SetBackground2(0.13836, 0.13836, 0.274);
m_vtkRenderer->GradientBackgroundOn();
///TODO : only on deformable objects with UseBounds (need vtk update)
......
......@@ -136,8 +136,8 @@ private:
///> Pointer to the Eigen's Conjugate gradient solver
Eigen::ConjugateGradient<SparseMatrixd> m_cgSolver;
std::vector<LinearProjectionConstraint> *m_FixedLinearProjConstraints;
std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints;
std::vector<LinearProjectionConstraint> *m_FixedLinearProjConstraints = nullptr;
std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints = nullptr;
};
} // imstk
......
......@@ -37,7 +37,7 @@ GaussSeidel::solve(Vectord& x)
return;
}
if (m_FixedLinearProjConstraints->size() == 0)
if (!m_FixedLinearProjConstraints)
{
this->gaussSeidelSolve(x);
}
......
......@@ -132,8 +132,8 @@ public:
private:
std::vector<LinearProjectionConstraint> *m_FixedLinearProjConstraints;
std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints;
std::vector<LinearProjectionConstraint> *m_FixedLinearProjConstraints = nullptr;
std::vector<LinearProjectionConstraint> *m_DynamicLinearProjConstraints = nullptr;
};
} // imstk
......
......@@ -106,6 +106,11 @@ public:
///
virtual bool isIterative() const = 0;
///
/// \brief Return the type of the solver
///
Type getType() { return m_type; };
protected:
Type m_type = Type::none; ///> Type of the scene object
double m_tolerance = 1.0e-4; ///> default tolerance
......
Markdown is supported
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