Commit 1c01c312 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Replace the way local matrices are

initialized from Vega matrices. The ode rhs
now considers damping only if the coefficients are
positive, non zero. It also accumulates contact forces
if any. The updateAttachedMesh in the VegaVolumetricMesh
class now takes an Eigen vector to do the update.
parent 08043796
Pipeline #4820 passed with stage
......@@ -38,19 +38,29 @@
VegaVolumetricMesh::VegaVolumetricMesh(bool generateMeshGraph) : generateGraph(generateMeshGraph)
{
}
//---------------------------------------------------------------------------
VegaVolumetricMesh::~VegaVolumetricMesh() {}
//---------------------------------------------------------------------------
std::shared_ptr<Graph> VegaVolumetricMesh::getMeshGraph()
{
return this->meshGraph;
}
//---------------------------------------------------------------------------
size_t VegaVolumetricMesh::getNumberOfVertices() const
{
return this->mesh->getNumVertices();
}
//---------------------------------------------------------------------------
size_t VegaVolumetricMesh::getNumberOfElements() const
{
return this->mesh->getNumElements();
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::attachSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, const double &radius)
{
// Keep copy of the mesh pointer
......@@ -63,6 +73,8 @@ void VegaVolumetricMesh::attachSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceM
}
this->generateWeigths(surfaceMesh,radius);
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::attachSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, const std::string& fileName, const double &radius)
{
// Keep copy of the mesh pointer
......@@ -75,10 +87,14 @@ void VegaVolumetricMesh::attachSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceM
}
this->readWeights(surfaceMesh,fileName,radius);
}
//---------------------------------------------------------------------------
std::shared_ptr<VolumetricMesh> VegaVolumetricMesh::getVegaMesh()
{
return this->mesh;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::setVegaMesh(std::shared_ptr<VolumetricMesh> newMesh)
{
this->mesh = newMesh;
......@@ -89,13 +105,16 @@ void VegaVolumetricMesh::setVegaMesh(std::shared_ptr<VolumetricMesh> newMesh)
meshGraph = std::make_shared<Graph>(*GenerateMeshGraph::Generate(this->mesh.get()));
}
}
void VegaVolumetricMesh::updateAttachedMeshes(double *q)
//---------------------------------------------------------------------------
void VegaVolumetricMesh::updateAttachedMeshes(const core::Vectord &q)
{
auto renderingMesh = this->getRenderingMesh();
auto *data = const_cast<double*>(q.data());
if(renderingMesh)
{
std::vector<core::Vec3d> displacements(renderingMesh->getNumberOfVertices(),core::Vec3d::Zero());
VolumetricMesh::interpolate(q,
VolumetricMesh::interpolate(data,
displacements.data()->data(),
renderingMesh->getNumberOfVertices(),
this->mesh->getNumElementVertices(),
......@@ -123,7 +142,7 @@ void VegaVolumetricMesh::updateAttachedMeshes(double *q)
auto collisionMesh = this->getCollisionMesh();
if(collisionMesh)
{
core::Matrix3MapType<double> displacementMap(q,3,this->mesh->getNumVertices());
core::Matrix3MapType<double> displacementMap(data,3,this->mesh->getNumVertices());
auto &restPositions = collisionMesh->getOrigVertices();
auto &vertices = collisionMesh->getVertices();
......@@ -139,45 +158,63 @@ void VegaVolumetricMesh::updateAttachedMeshes(double *q)
}
}
}
//---------------------------------------------------------------------------
const std::unordered_map<size_t,size_t>& VegaVolumetricMesh::getVertexMap() const
{
return this->vertexMap;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::setVertexMap(const std::unordered_map<size_t,size_t>& map)
{
this->vertexMap = map;
}
//---------------------------------------------------------------------------
const std::vector< size_t >& VegaVolumetricMesh::getFixedVertices() const
{
return this->fixedVertices;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::setFixedVertices(const std::vector< size_t >& dofs)
{
this->fixedVertices.clear();
this->fixedVertices = dofs;
}
//---------------------------------------------------------------------------
std::shared_ptr<SurfaceMesh> VegaVolumetricMesh::getAttachedMesh(const size_t i)
{
return i < this->attachedMeshes.size() ?
this->attachedMeshes.at(i) :
nullptr;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::setRenderDetail(int i, std::shared_ptr< RenderDetail > newRenderDetail)
{
this->attachedMeshes.at(i)->setRenderDetail(newRenderDetail);
}
//---------------------------------------------------------------------------
std::shared_ptr<SurfaceMesh> VegaVolumetricMesh::getRenderingMesh()
{
return attachedMeshes.size() > 1
? this->attachedMeshes.at(attachedMeshes.size() - 1)
: nullptr;
}
//---------------------------------------------------------------------------
std::shared_ptr<SurfaceMesh> VegaVolumetricMesh::getCollisionMesh()
{
return attachedMeshes.size() > 0
? this->attachedMeshes.at(0)
: nullptr;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::saveWeights(std::shared_ptr<SurfaceMesh> surfaceMesh, const std::string& filename) const
{
auto vertices = this->attachedVertices.at(surfaceMesh);
......@@ -199,6 +236,8 @@ void VegaVolumetricMesh::saveWeights(std::shared_ptr<SurfaceMesh> surfaceMesh, c
fileStream << std::endl;
fileStream.close();
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::readWeights(std::shared_ptr<SurfaceMesh> surfaceMesh, const std::string& filename, const double radius)
{
std::ifstream fileStream(filename.c_str());
......@@ -239,6 +278,8 @@ void VegaVolumetricMesh::readWeights(std::shared_ptr<SurfaceMesh> surfaceMesh, c
std::cout << "\tTotal # of weights read: " << weigths.size() / this->mesh->getNumElementVertices() << std::endl;
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::generateWeigths(std::shared_ptr<SurfaceMesh> surfaceMesh, double radius, const bool saveToDisk, const std::string& filename)
{
std::cerr << "Generating weights..." << std::endl;
......@@ -305,14 +346,20 @@ void VegaVolumetricMesh::generateWeigths(std::shared_ptr<SurfaceMesh> surfaceMes
this->saveWeights(surfaceMesh, filename);
}
}
//---------------------------------------------------------------------------
const std::vector< double >& VegaVolumetricMesh::getAttachedWeights(std::shared_ptr<SurfaceMesh> surfaceMesh) const
{
return this->attachedWeights.at(surfaceMesh);
}
//---------------------------------------------------------------------------
const std::vector< int >& VegaVolumetricMesh::getAttachedVertices(std::shared_ptr<SurfaceMesh> surfaceMesh) const
{
return this->attachedVertices.at(surfaceMesh);
}
//---------------------------------------------------------------------------
void VegaVolumetricMesh::translate(const Eigen::Translation3d& translation, bool setInitialPoints)
{
for(auto meshes : this->attachedMeshes)
......
......@@ -48,6 +48,9 @@ public:
///
/// \brief Constructor
///
/// \param generateMeshGraph True if you want to generate a mesh graph. The mesh
/// graphs is used by the time stepping method to apply Lagrangian damping.
///
VegaVolumetricMesh(bool generateMeshGraph = true);
///
......@@ -82,19 +85,19 @@ public:
void attachSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh, const std::string &fileName, const double &radius = 5.0);
///
/// \brief Return mesh
/// \brief Return underlying Vega mesh data structure.
///
std::shared_ptr<VolumetricMesh> getVegaMesh();
///
/// \brief Sets the vega mesh
/// \brief Sets the Vega underlying mesh.
///
void setVegaMesh(std::shared_ptr<VolumetricMesh> newMesh);
///
/// \brief Update nodes to local arrays
///
void updateAttachedMeshes(double *q);
void updateAttachedMeshes(const core::Vectord &q);
///
/// \brief Return the vertex map
......@@ -169,7 +172,8 @@ private:
// Vega mesh base object
std::shared_ptr<VolumetricMesh> mesh;
// Vega mesh graph
// Vega mesh graph. The mesh graphs is used by the time stepping method to
// apply Lagrangian damping.
std::shared_ptr<Graph> meshGraph;
// Generate graph for this mesh
......
......@@ -24,6 +24,10 @@
// SimMedTK includes
#include "SceneModels/DeformableSceneObject.h"
#include <cmath> // for std::isfinite()
DeformableSceneObject::DeformableSceneObject():
OdeSystem(),
integrationScheme(TimeIntegrator::ImplicitEuler)
......@@ -36,9 +40,9 @@ void DeformableSceneObject::applyContactForces()
for(const auto & cf : this->getContactForces())
{
auto i = cf.first;
f(i) += cf.second(0);
f(i + 1) += cf.second(1);
f(i + 2) += cf.second(2);
this->f(i) += cf.second(0);
this->f(i + 1) += cf.second(1);
this->f(i + 2) += cf.second(2);
}
}
......@@ -92,10 +96,18 @@ void DeformableSceneObject::update(const double dt)
this->odeSolver->solve(*this->currentState, *this->newState, dt);
// TODO: Check state validity
if(!std::isfinite(this->newState->getPositions().sum()) ||
!std::isfinite(this->newState->getVelocities().sum()) )
{
// TODO: log this and throw exception, this is a fatal error
return;
}
this->currentState.swap(this->previousState);
this->currentState.swap(this->newState);
// TODO: Check state validity
this->updateMesh();
}
//---------------------------------------------------------------------------
......
......@@ -68,6 +68,11 @@ public:
///
void update(const double dt);
///
/// \brief Update states
///
virtual void updateMesh() = 0;
///
/// \brief Reset the current state to the initial state
///
......
......@@ -106,7 +106,7 @@ public:
///
/// \brief Constructor
///
VegaConfiguration(const std::string &configurationFile, bool verbose = false);
VegaConfiguration(const std::string &configurationFile, bool verbose = true);
std::string vegaConfigFile; ///> Store configuration file.
......@@ -342,7 +342,7 @@ bool VegaFEMDeformableSceneObject::configure(const std::string &configFile)
this->setMassMatrix();
this->setTangentStiffnessMatrix();
this->setDampingMatrix();
this->setDampingMatrices();
this->setOdeRHS();
size_t numNodes = this->volumetricMesh->getNumberOfVertices();
......@@ -367,28 +367,25 @@ void VegaFEMDeformableSceneObject::initMassMatrix(bool saveToDisk)
&matrix,
true);
this->vegaMassMatrix.reset(matrix);
auto nnz = this->vegaMassMatrix->GetNumEntries();
this->massMatrixRowPointers.resize(this->vegaMassMatrix->GetNumRows()+1);
this->massMatrixColIndices.resize(nnz);
this->massMatrixValues.resize(nnz);
this->vegaMassMatrix->
GenerateCompressedRowMajorFormat(this->massMatrixValues.data(),
this->massMatrixRowPointers.data(),
this->massMatrixColIndices.data());
// Construct the Eigen mass matrix by mapping the arrays
this->M = Eigen::MappedSparseMatrix<double,Eigen::RowMajor>(
this->vegaMassMatrix->GetNumRows(),
this->vegaMassMatrix->GetNumColumns(),
this->vegaMassMatrix->GetNumEntries(),
this->massMatrixRowPointers.data(),
this->massMatrixColIndices.data(),
this->massMatrixValues.data());
auto rowLengths = matrix->GetRowLengths();
auto nonZeroValues = matrix->GetEntries();
auto columnIndices = matrix->GetColumnIndices();
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(matrix->GetNumEntries());
for(int i = 0, end = matrix->GetNumRows(); i < end; ++i)
{
for(int k = 0, k_end = rowLengths[i]; k < k_end; ++k)
{
triplets.emplace_back(i,columnIndices[i][k],nonZeroValues[i][k]);
}
}
this->M.resize(matrix->GetNumRows(),
matrix->GetNumColumns());
this->M.setFromTriplets(std::begin(triplets),std::end(triplets));
this->M.makeCompressed();
this->vegaMassMatrix.reset(matrix);
if(saveToDisk)
{
char name[] = "ComputedMassMatrix.mass";
......@@ -414,7 +411,6 @@ void VegaFEMDeformableSceneObject::initTangentStiffnessMatrix()
return;
}
this->vegaTangentStiffnessMatrix.reset(matrix);
if(!this->vegaMassMatrix)
{
......@@ -422,7 +418,7 @@ void VegaFEMDeformableSceneObject::initTangentStiffnessMatrix()
return;
}
this->vegaTangentStiffnessMatrix->BuildSubMatrixIndices(*this->vegaMassMatrix.get());
matrix->BuildSubMatrixIndices(*this->vegaMassMatrix.get());
if(!this->dampingMatrix)
{
......@@ -430,41 +426,49 @@ void VegaFEMDeformableSceneObject::initTangentStiffnessMatrix()
return;
}
this->vegaTangentStiffnessMatrix->BuildSubMatrixIndices(*this->dampingMatrix.get(), 1);
auto nnz = this->vegaTangentStiffnessMatrix->GetNumEntries();
this->tangentStiffnessMatrixRowPointers.resize(this->vegaTangentStiffnessMatrix->GetNumRows()+1);
this->tangentStiffnessMatrixColIndices.resize(nnz);
this->tangentStiffnessMatrixValues.resize(nnz);
matrix->BuildSubMatrixIndices(*this->dampingMatrix.get(), 1);
this->vegaTangentStiffnessMatrix->
GenerateCompressedRowMajorFormat(this->tangentStiffnessMatrixValues.data(),
this->tangentStiffnessMatrixRowPointers.data(),
this->tangentStiffnessMatrixColIndices.data());
// Construct the eigen stiffness matrix by mapping the arrays
this->K = Eigen::MappedSparseMatrix<double,Eigen::RowMajor>(
this->vegaTangentStiffnessMatrix->GetNumRows(),
this->vegaTangentStiffnessMatrix->GetNumColumns(),
nnz,
this->tangentStiffnessMatrixRowPointers.data(),
this->tangentStiffnessMatrixColIndices.data(),
this->tangentStiffnessMatrixValues.data());
auto rowLengths = matrix->GetRowLengths();
auto columnIndices = matrix->GetColumnIndices();
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(matrix->GetNumEntries());
for(int i = 0, end = matrix->GetNumRows(); i < end; ++i)
{
for(int k = 0, k_end = rowLengths[i]; k < k_end; ++k)
{
triplets.emplace_back(i,columnIndices[i][k],0.001);
}
}
this->K.resize(matrix->GetNumRows(),
matrix->GetNumColumns());
this->K.setFromTriplets(std::begin(triplets),std::end(triplets));
this->K.makeCompressed();
this->vegaTangentStiffnessMatrix.reset(matrix);
const auto &dampingStiffnessCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingStiffnessCoefficient");
const auto &dampingMassCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingMassCoefficient");
// Initialize the Raleigh damping matrix
this->C = this->M*dampingMassCoefficient + this->K*dampingStiffnessCoefficient;
}
//---------------------------------------------------------------------------
void VegaFEMDeformableSceneObject::initDampingMatrix()
{
auto dampingLaplacianCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingLaplacianCoefficient");
if(!(dampingLaplacianCoefficient > 0.0))
{
/// TODO: add to log
return;
}
auto meshGraph = this->volumetricMesh->getMeshGraph();
if(!meshGraph)
......@@ -482,30 +486,28 @@ void VegaFEMDeformableSceneObject::initDampingMatrix()
return;
}
this->dampingMatrix.reset(matrix);
matrix->ScalarMultiply(dampingLaplacianCoefficient);
auto rowLengths = matrix->GetRowLengths();
auto nonZeroValues = matrix->GetEntries();
auto columnIndices = matrix->GetColumnIndices();
std::vector<Eigen::Triplet<double>> triplets;
triplets.reserve(matrix->GetNumEntries());
for(int i = 0, end = matrix->GetNumRows(); i < end; ++i)
{
for(int k = 0, k_end = rowLengths[i]; k < k_end; ++k)
{
triplets.emplace_back(i,columnIndices[i][k],nonZeroValues[i][k]);
}
}
this->dampingMatrix->ScalarMultiply(
this->vegaFemConfig->floatsOptionMap.at("dampingLaplacianCoefficient"));
auto nnz = this->dampingMatrix->GetNumEntries();
this->dampingMatrixColPointers.resize(this->dampingMatrix->GetNumRows()+1);
this->dampingMatrixColIndices.resize(nnz);
this->dampingMatrixValues.resize(nnz);
this->dampingMatrix->
GenerateCompressedRowMajorFormat(this->dampingMatrixValues.data(),
this->dampingMatrixColPointers.data(),
this->dampingMatrixColIndices.data());
// Construct the Eigen damping matrix by mapping the arrays
this->D = Eigen::MappedSparseMatrix<double,Eigen::RowMajor>(
this->dampingMatrix->GetNumRows(),
this->dampingMatrix->GetNumColumns(),
this->dampingMatrix->GetNumEntries(),
this->dampingMatrixColPointers.data(),
this->dampingMatrixColIndices.data(),
this->dampingMatrixValues.data());
this->D.resize(matrix->GetNumRows(),
matrix->GetNumColumns());
this->D.setFromTriplets(std::begin(triplets),std::end(triplets));
this->D.makeCompressed();
this->dampingMatrix.reset(matrix);
}
//---------------------------------------------------------------------------
......@@ -689,9 +691,9 @@ void VegaFEMDeformableSceneObject::initForceModel()
stVKInternalForces.get(),
stVKStiffnessMatrix.get());
auto &uInitial = this->getInitialState()->getPositions();
auto &uCurrent = this->currentState->getPositions();
this->forceModel->GetInternalForce(uInitial.data(), uCurrent.data());
// auto &uInitial = this->getInitialState()->getPositions();
// auto &uCurrent = this->currentState->getPositions();
// this->forceModel->GetInternalForce(uInitial.data(), uCurrent.data());
break;
}
......@@ -770,7 +772,7 @@ std::vector< std::size_t > VegaFEMDeformableSceneObject::loadBoundaryConditions(
//---------------------------------------------------------------------------
void VegaFEMDeformableSceneObject::updateValuesFromMatrix(std::shared_ptr<SparseMatrix> matrix,
std::vector<double> &values)
double *values)
{
auto rowLengths = matrix->GetRowLengths();
auto nonZeroValues = matrix->GetEntries();
......@@ -794,9 +796,32 @@ void VegaFEMDeformableSceneObject::updateValuesFromMatrix(std::shared_ptr<Sparse
//---------------------------------------------------------------------------
void VegaFEMDeformableSceneObject::setOdeRHS()
{
auto odeRHS = [this](const OdeSystemState & s) -> const core::Vectord&
const auto &dampingStiffnessCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingStiffnessCoefficient");
const auto &dampingMassCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingMassCoefficient");
auto odeRHS = [&,this](const OdeSystemState & s) -> const core::Vectord&
{
this->f = -this->C * s.getVelocities() - this->K * s.getPositions();
double *data = const_cast<double*>(s.getPositions().data());
this->forceModel->GetInternalForce(data,this->f.data());
// Add the Raleigh damping force
if(dampingMassCoefficient > 0)
{
this->f += dampingMassCoefficient*this->M*s.getVelocities();
}
if(dampingStiffnessCoefficient > 0)
{
this->f += dampingStiffnessCoefficient*this->K*s.getVelocities();
}
// Apply contact forces
this->applyContactForces();
// Vega returns the negative of the force action on the material
this->f *= -1;
return this->f;
};
this->setFunction(odeRHS);
......@@ -813,8 +838,7 @@ void VegaFEMDeformableSceneObject::setTangentStiffnessMatrix()
GetTangentStiffnessMatrix(data,this->vegaTangentStiffnessMatrix.get());
this->updateValuesFromMatrix(this->vegaTangentStiffnessMatrix,
this->tangentStiffnessMatrixValues);
std::cout << this->K << std::endl;
this->K.valuePtr());
return this->K;
};
this->setJaconbianFx(tangentStiffness);
......@@ -831,19 +855,42 @@ void VegaFEMDeformableSceneObject::setMassMatrix()
}
//---------------------------------------------------------------------------
void VegaFEMDeformableSceneObject::setDampingMatrix()
void VegaFEMDeformableSceneObject::setDampingMatrices()
{
const auto &dampingStiffnessCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingStiffnessCoefficient");
this->vegaFemConfig->floatsOptionMap.at("dampingStiffnessCoefficient");
const auto &dampingMassCoefficient =
this->vegaFemConfig->floatsOptionMap.at("dampingMassCoefficient");
this->vegaFemConfig->floatsOptionMap.at("dampingMassCoefficient");
auto raleighDamping = [&,this](const OdeSystemState & /*s*/) -> const core::SparseMatrixd&
{
this->C = this->M*dampingMassCoefficient;
this->C += this->K*dampingStiffnessCoefficient;
if(dampingMassCoefficient > 0)
{
this->C = dampingMassCoefficient*this->M;
if(dampingStiffnessCoefficient > 0)
{
this->C += dampingStiffnessCoefficient*this->K;
}
}
else if(dampingStiffnessCoefficient > 0)
{
this->C = dampingStiffnessCoefficient*this->K;
}
return this->C;
};
this->setJaconbianFv(raleighDamping);
if(this->dampingMatrix)
{
auto lagrangianDamping = [this](const OdeSystemState & /*s*/) -> const core::SparseMatrixd&
{
return this->D;
};
this->setDamping(lagrangianDamping);
}
}
void VegaFEMDeformableSceneObject::updateMesh()
{
this->volumetricMesh->updateAttachedMeshes(this->currentState->getPositions());
}
......@@ -121,7 +121,7 @@ public:
/// \brief Helper function to interface Vega sparse matrices with Eigen matrices.
///
void updateValuesFromMatrix(std::shared_ptr<SparseMatrix> matrix,
std::vector<double> &values);
double *values);
///
/// \brief Set the ODE system right hand side function.
......@@ -141,7 +141,12 @@ public:
///
/// \brief Set the the Raleigh damping matrix function to be evaluated by the ODE solver.
///
void setDampingMatrix();
void setDampingMatrices();