Commit b97c0778 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

ENH: Finish interface to vega force models.

parent 84f3390c
Pipeline #4779 passed with stage
......@@ -83,7 +83,7 @@ void VegaVolumetricMesh::setVegaMesh(std::shared_ptr<VolumetricMesh> newMesh)
{
this->mesh = newMesh;
if(nullptr != this->mesh && generateGraph)
if(nullptr != this->mesh && this->generateGraph)
{
meshGraph.reset();
meshGraph = std::make_shared<Graph>(*GenerateMeshGraph::Generate(this->mesh.get()));
......
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
// Rensselaer Polytechnic Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
// SimMedTK includes
#include "Simulators/DeformableSceneObject.h"
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
// Rensselaer Polytechnic Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
#ifndef SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
#define SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
// SimMedTK includes
#include "Core/SceneObject.h"
#include "TimeIntegrators/TimeIntegrator.h"
#include "TimeIntegrators/BackwarEuler.h"
#include "TimeIntegrators/ForwardEuler.h"
///
/// \brief Base class for all deformable scene objects.
///
class DeformableSceneObject : public SceneObject, public OdeSystem
{
public:
///
/// \brief Constructor
///
DeformableSceneObject() :
OdeSystem(),
integrationScheme(TimeIntegrator::ImplicitEuler)
{
}
///
/// \brief Destructor
///
virtual ~DeformableSceneObject() = default;
///
/// \brief Append the contact forces (if any) to external forces
///
void applyContactForces()
{
for(const auto &cf : this->getContactForces())
{
int i = cf.first;
f(i) += cf.second(0);
f(i+1) += cf.second(1);
f(i+2) += cf.second(2);
}
}
///
/// \brief Set the integration scheme used to solve the ODE system.
///
void setTimeIntegrator(TimeIntegrator::IntegratorType integrator)
{
this->integrationScheme = integrator;
}
///
/// \brief Initialize the ode solver.
///
bool init()
{
auto thisPointer = this->safeDownCast<DeformableSceneObject>();
switch(integrationScheme)
{
case TimeIntegrator::ImplicitEuler:
{
this->odeSolver = std::make_shared<BackwardEuler>(thisPointer);
}
case TimeIntegrator::ExplicitEuler:
{
this->odeSolver = std::make_shared<ForwardEuler>(thisPointer);
}
default:
{
std::cerr << "Invalid time integration scheme." << std::endl;
}
return false;
}
}
///
/// \brief Update states
///
void update(double dt)
{
if(!this->odeSolver)
{
std::cerr << "Ode solver needs to be set." << std::endl;
return;
}
this->odeSolver->solve(this->currentState,this->newState,dt);
this->currentState.swap(this->previousState);
this->currentState.swap(this->newState);
// TODO: Check state validity
}
///
/// \brief Reset the current state to the initial state
///
virtual void resetToInitialState()
{
*this->currentState = *this->initialState;
*this->previousState = *this->initialState;
}
///
/// \brief Return the current state.
///
std::shared_ptr<OdeSystemState> getCurrentState()
{
this->currentState;
}
///
/// \brief Return the previous state.
///
std::shared_ptr<OdeSystemState> getPreviousState()
{
this->previousState;
}
protected:
std::shared_ptr<TimeIntegrator> odeSolver; ///> Integration scheme
// Consecutive system states
std::shared_ptr<OdeSystemState> currentState; ///> Current model state
std::shared_ptr<OdeSystemState> previousState;///> Previous model state
std::shared_ptr<OdeSystemState> newState; ///> Storage for the next state
core::SparseMatrixd M; ///> Mass matrix
core::SparseMatrixd C; ///> Damping matrix
core::SparseMatrixd K; ///> Stiffness matrix
core::Vectord f; ///> Force loads
TimeIntegrator::IntegratorType integrationScheme; ///> Integration scheme used.
};
#endif // SMVEGAFEMSCENEOBJECT_DEFORMABLE_H
......@@ -41,10 +41,6 @@ SceneObjectDeformable::SceneObjectDeformable() :
RenderDelegate::RendererType::VTK));
}
SceneObjectDeformable::~SceneObjectDeformable()
{
}
void SceneObjectDeformable::applyContactForces()
{
for(const auto &cf : this->getContactForces())
......
......@@ -41,15 +41,13 @@ class SceneObjectDeformable : public SceneObject, public OdeSystem
{
public:
///
/// \brief Constructor
///
SceneObjectDeformable();
/// \brief Destructor
~SceneObjectDeformable();
/// \brief rest the object to inital configuration
/// and reset initial states
virtual void resetToInitialState() = 0;
virtual ~SceneObjectDeformable() = default;
/// \brief Load specified meshes
virtual void loadVolumeMesh() = 0;
......@@ -168,6 +166,9 @@ public:
}
}
///
/// \brief Update states
///
void update(double dt)
{
if(!this->odeSolver)
......@@ -176,14 +177,39 @@ public:
return;
}
this->odeSolver->solve(currentState,newState,dt);
this->odeSolver->solve(this->currentState,this->newState,dt);
currentState.swap(previousState);
currentState.swap(newState);
this->currentState.swap(this->previousState);
this->currentState.swap(this->newState);
// TODO: Check state validity
}
///
/// \brief Reset the current state to the initial state
///
virtual void resetToInitialState()
{
*this->currentState = *this->initialState;
*this->previousState = *this->initialState;
}
///
/// \brief Return the current state.
///
std::shared_ptr<OdeSystemState> getCurrentState()
{
this->currentState;
}
///
/// \brief Return the previous state.
///
std::shared_ptr<OdeSystemState> getPreviousState()
{
this->previousState;
}
protected:
friend class SceneObjectDeformableRenderDelegate;
......@@ -224,12 +250,10 @@ protected:
std::shared_ptr<OdeSystemState> previousState;
std::shared_ptr<OdeSystemState> newState;
std::shared_ptr<OdeSystemState> initialState;
core::Matrixd M; ///> Mass matrix
core::Matrixd C; ///> Damping matrix
core::Matrixd K; ///> Stiffness matrix
core::Vectord f; ///> Force loads
core::SparseMatrixd M; ///> Mass matrix
core::SparseMatrixd C; ///> Damping matrix
core::SparseMatrixd K; ///> Stiffness matrix
core::Vectord f; ///> Force loads
TimeIntegrator::IntegratorType integrationScheme;
......
// This file is part of the SimMedTK project.
// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
// Rensselaer Polytechnic Institute
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
//---------------------------------------------------------------------------
//
// Authors:
//
// Contact:
//---------------------------------------------------------------------------
#include "VegaFEMDeformableSceneObject.h"
This diff is collapsed.
......@@ -126,9 +126,9 @@ public:
///
/// \brief Get the initial velocities and positions of the system.
///
const std::shared_ptr<OdeSystemState> getState() const
const std::shared_ptr<OdeSystemState> getInitialState() const
{
return this->state;
return this->initialState;
}
///
......@@ -136,9 +136,9 @@ public:
///
/// \param newState Current positions and velocities.
///
void setState(std::shared_ptr<OdeSystemState> newState) const
void setInitialState(std::shared_ptr<OdeSystemState> newState) const
{
this->state = newState;
this->initialState = newState;
}
private:
......@@ -147,7 +147,7 @@ private:
MatrixFunctionType Mass; ///> Function to evaluate the mass matrix.
FunctionType F; ///> Right hand side function
std::shared_ptr<OdeSystemState> state; ///> Initial state of the system.
std::shared_ptr<OdeSystemState> initialState; ///> Initial state of the system.
};
#endif // ODESYSTEM_H
......@@ -25,6 +25,7 @@
#define ODESYSTEMSTATE_H
#include "Core/Vector.h"
#include "Core/Matrix.h"
///
/// \brief Stores the state of the differential equation.
......@@ -32,11 +33,18 @@
class OdeSystemState
{
public:
/// Default constructor/destructor
///
/// \brief Default constructor/destructor.
///
OdeSystemState() = default;
~OdeSystemState() = default;
/// Constructor with system size.
OdeSystemState(const size_t size)
{
this->resize(size);
}
///
/// \brief Set the derivative with respect to v of the right hand side.
///
......@@ -82,14 +90,78 @@ public:
///
/// \return Reference to velocities.
///
void resize(size_t size)
void resize(const size_t size)
{
positions.resize(size);
velocities.resize(size);
positions.setZero();
velocities.setZero();
}
///
/// \brief Return vector containing the indices of fixed dofs
///
/// \return Reference to vector indices.
///
void setBoundaryConditions(const std::vector<size_t> &boundaryConditions)
{
this->fixedVertices = boundaryConditions;
}
///
/// \brief Apply boundary conditions to sparse matrix.
///
/// \param M Sparse matrix container.
/// \param withCompliance True if the fixed vertices should have complieance.
///
void applyBoundaryConditions(core::SparseMatrixd &M, bool withCompliance = true)
{
double compliance = withCompliance ? 1.0 : 0.0;
for(auto &index : this->fixedVertices)
{
M.middleRows(index,1).setZero();
M.middleCols(index,1).setZero();
M.coeffRef(index,index) = compliance;
}
}
///
/// \brief Apply boundary conditions to dense matrix.
///
/// \param M Dense matrix container.
/// \param withCompliance True if the fixed vertices should have complieance.
///
void applyBoundaryConditions(core::Matrixd &M, bool withCompliance = true)
{
double compliance = withCompliance ? 1.0 : 0.0;
for(auto &index : this->fixedVertices)
{
M.middleRows(index,1).setZero();
M.middleCols(index,1).setZero();
M(index,index) = compliance;
}
}
///
/// \brief Apply boundary conditions to a vector.
///
/// \param x vector container.
///
void applyBoundaryConditions(core::Vectord &x)
{
for(auto &index : this->fixedVertices)
{
x(index) = 0.0;
}
}
private:
core::Vectord positions; ///> State position.
core::Vectord velocities; ///> State velocity.
std::vector<size_t> fixedVertices; ///> Fixed dof ids
};
#endif // ODESYSTEMSTATE_H
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