Skip to content

Commit

Permalink
Add all built-in solvers and several ordering methods (#1)
Browse files Browse the repository at this point in the history
* Introduce LLT, LU and QR solvers

* Choice of the ordering method

* Add metis support

* Optimal plugins

* Link with Metis

* find package metis

* File encoding

* Another try

* Find metis for v21.12

* Fix CMakeLists.txt

* unused variable

* Remove find_package of metis

* Missing template keyword

* find_package SofaFramework

* Make sure the linear solver module is found

* Replace hard-coded value

* Condition on finding the metis include file

* Fix missing guard

* Remove warning directive
  • Loading branch information
alxbilger authored Apr 14, 2022
1 parent f807b77 commit 703c2ca
Show file tree
Hide file tree
Showing 28 changed files with 710 additions and 121 deletions.
42 changes: 34 additions & 8 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,29 +1,55 @@
cmake_minimum_required(VERSION 3.12)
project(EigenLinearSolvers VERSION 1.0 LANGUAGES CXX)

find_package(Sofa.Component.LinearSolver QUIET)
if (NOT Sofa.Component.LinearSolver_FOUND)
find_package(SofaBaseLinearSolver QUIET)
find_package(SofaFramework REQUIRED)
find_package(Sofa.Component.LinearSolver.Direct QUIET)
if (NOT Sofa.Component.LinearSolver.Direct_FOUND)
find_package(SofaSparseSolver QUIET)
endif()

if (NOT Sofa.Component.LinearSolver.Direct_FOUND AND NOT SofaSparseSolver_FOUND)
message(FATAL "Cannot find the linear solver module in SOFA")
endif()

# List all files
set(EIGENLINEARSOLVERS_SRC_DIR src/${PROJECT_NAME})
set(HEADER_FILES
${EIGENLINEARSOLVERS_SRC_DIR}/config.h.in
${EIGENLINEARSOLVERS_SRC_DIR}/FindMetis.h

${EIGENLINEARSOLVERS_SRC_DIR}/SimplicialLDLTTraits.h
${EIGENLINEARSOLVERS_SRC_DIR}/SimplicialLLTTraits.h
${EIGENLINEARSOLVERS_SRC_DIR}/SparseLUTraits.h
${EIGENLINEARSOLVERS_SRC_DIR}/SparseQRTraits.h

${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].inl

${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver[CRS].h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenDirectSparseSolver[CRS].inl

${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].inl

${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT[CRS].h

${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU[CRS].h

${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR.h
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR[CRS].h
)
set(SOURCE_FILES
${EIGENLINEARSOLVERS_SRC_DIR}/init.cpp

${EIGENLINEARSOLVERS_SRC_DIR}/EigenConjugateGradient[CRS].cpp
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLDLT[CRS].cpp
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSimplicialLLT[CRS].cpp
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseLU[CRS].cpp
${EIGENLINEARSOLVERS_SRC_DIR}/EigenSparseQR[CRS].cpp
)
set(README_FILES
README.md
Expand All @@ -33,10 +59,10 @@ set(README_FILES
add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${README_FILES})

# Link the plugin library to its dependency(ies).
if (Sofa.Component.LinearSolver_FOUND)
target_link_libraries(${PROJECT_NAME} Sofa.Component.LinearSolver)
elseif(SofaBaseLinearSolver_FOUND)
target_link_libraries(${PROJECT_NAME} SofaBaseLinearSolver)
if (Sofa.Component.LinearSolver.Direct_FOUND)
target_link_libraries(${PROJECT_NAME} Sofa.Component.LinearSolver.Direct)
elseif(SofaSparseSolver_FOUND)
target_link_libraries(${PROJECT_NAME} SofaSparseSolver)
endif()

# Create package Config, Version & Target files.
Expand Down
17 changes: 9 additions & 8 deletions scenes/Cylinder_EigenSimplicialLDLT.scn
Original file line number Diff line number Diff line change
@@ -1,17 +1,18 @@
<Node name="root" gravity="-1.8 0 100" dt="0.001">
<RequiredPlugin name="SofaBoundaryCondition"/>
<RequiredPlugin name="SofaImplicitOdeSolver"/>
<RequiredPlugin name="SofaLoader"/>
<RequiredPlugin name="SofaMiscForceField"/>
<RequiredPlugin name="SofaOpenglVisual"/>
<RequiredPlugin name="SofaSimpleFem"/>
<RequiredPlugin name="EigenLinearSolvers"/>
<RequiredPlugin name="EigenLinearSolvers"/> <!-- Needed to use components [EigenSimplicialLDLT] -->
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Mapping.Linear"/> <!-- Needed to use components [BarycentricMapping] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [MeshMatrixMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <!-- Needed to use components [TetrahedronFEMForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.GL.Component.Rendering3D"/> <!-- Needed to use components [OglModel] -->

<Node name="DeformableObject">

<EulerImplicitSolver name="odeImplicitSolver" />

<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3d"/>
<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3d" ordering="Metis"/>

<MeshGmshLoader name="loader" filename="mesh/truthcylinder1.msh" />
<TetrahedronSetTopologyContainer src="@loader" name="topologyContainer"/>
Expand Down
17 changes: 10 additions & 7 deletions scenes/FEMBAR_EigenConjugateGradient.scn
Original file line number Diff line number Diff line change
@@ -1,13 +1,16 @@
<Node name="root" dt="0.02" gravity="0 -10 0">

<Node name="plugins">
<RequiredPlugin name="SofaBoundaryCondition"/>
<RequiredPlugin name="SofaImplicitOdeSolver"/>
<RequiredPlugin name="SofaEngine"/>
<RequiredPlugin name="EigenLinearSolvers"/> <!-- Needed to use components [EigenConjugateGradient] -->
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <!-- Needed to use components [HexahedronFEMForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
<RequiredPlugin name="SofaEngine"/> <!-- Needed to use components [BoxROI] -->

<RequiredPlugin name="SofaSimpleFem"/>
<RequiredPlugin name="SofaPreconditioner"/>
<RequiredPlugin name="EigenLinearSolvers"/>
</Node>

<VisualStyle displayFlags="showBehaviorModels showForceFields" />
Expand All @@ -16,7 +19,7 @@
<DefaultVisualManagerLoop name="visualLoop"/>

<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
<EigenConjugateGradient template="EigenSparseMatrixd" iterations="1000" tolerance="1e-9" preconditioner="diagonal"/>
<EigenConjugateGradient template="CompressedRowSparseMatrixMat3x3d" iterations="1000" tolerance="1e-9" preconditioner="diagonal"/>

<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
<MechanicalObject name="DoFs" />
Expand Down
17 changes: 10 additions & 7 deletions scenes/FEMBAR_EigenSimplicialLDLT.scn
Original file line number Diff line number Diff line change
@@ -1,13 +1,16 @@
<Node name="root" dt="0.02" gravity="0 -10 0">

<Node name="plugins">
<RequiredPlugin name="SofaBoundaryCondition"/>
<RequiredPlugin name="SofaImplicitOdeSolver"/>
<RequiredPlugin name="SofaEngine"/>
<RequiredPlugin name="EigenLinearSolvers"/> <!-- Needed to use components [EigenSimplicialLDLT] -->
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <!-- Needed to use components [HexahedronFEMForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
<RequiredPlugin name="SofaEngine"/> <!-- Needed to use components [BoxROI] -->

<RequiredPlugin name="SofaSimpleFem"/>
<RequiredPlugin name="SofaPreconditioner"/>
<RequiredPlugin name="EigenLinearSolvers"/>
</Node>

<VisualStyle displayFlags="showBehaviorModels showForceFields" />
Expand All @@ -16,7 +19,7 @@
<DefaultVisualManagerLoop name="visualLoop"/>

<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3d"/>
<EigenSimplicialLDLT template="CompressedRowSparseMatrixMat3x3d" ordering="Natural"/>

<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
<MechanicalObject name="DoFs" />
Expand Down
33 changes: 33 additions & 0 deletions scenes/FEMBAR_EigenSimplicialLLT.scn
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
<Node name="root" dt="0.02" gravity="0 -10 0">

<Node name="plugins">
<RequiredPlugin name="EigenLinearSolvers"/> <!-- Needed to use components [EigenSimplicialLLT] -->
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <!-- Needed to use components [HexahedronFEMForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
<RequiredPlugin name="SofaEngine"/> <!-- Needed to use components [BoxROI] -->

</Node>

<VisualStyle displayFlags="showBehaviorModels showForceFields" />

<DefaultAnimationLoop name="animationLoop"/>
<DefaultVisualManagerLoop name="visualLoop"/>

<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
<EigenSimplicialLLT template="CompressedRowSparseMatrixMat3x3d" ordering="AMD"/>

<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
<MechanicalObject name="DoFs" />

<UniformMass name="mass" totalMass="320" />
<HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.45" method="large" />

<BoxROI name="box" box="-10 -1 -0.0001 -5 4 0.0001"/>
<FixedConstraint indices="@box.indices" />

</Node>
33 changes: 33 additions & 0 deletions scenes/FEMBAR_EigenSparseLU.scn
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
<Node name="root" dt="0.02" gravity="0 -10 0">

<Node name="plugins">
<RequiredPlugin name="EigenLinearSolvers"/> <!-- Needed to use components [EigenSparseLU] -->
<RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedConstraint] -->
<RequiredPlugin name="Sofa.Component.Mass"/> <!-- Needed to use components [UniformMass] -->
<RequiredPlugin name="Sofa.Component.ODESolver.Backward"/> <!-- Needed to use components [EulerImplicitSolver] -->
<RequiredPlugin name="Sofa.Component.SolidMechanics.FEM.Elastic"/> <!-- Needed to use components [HexahedronFEMForceField] -->
<RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
<RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
<RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
<RequiredPlugin name="SofaEngine"/> <!-- Needed to use components [BoxROI] -->

</Node>

<VisualStyle displayFlags="showBehaviorModels showForceFields" />

<DefaultAnimationLoop name="animationLoop"/>
<DefaultVisualManagerLoop name="visualLoop"/>

<EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />
<EigenSparseLU template="CompressedRowSparseMatrixMat3x3d" ordering="Natural"/>

<RegularGridTopology name="grid" nx="4" ny="4" nz="20" xmin="-9" xmax="-6" ymin="0" ymax="3" zmin="0" zmax="19" />
<MechanicalObject name="DoFs" />

<UniformMass name="mass" totalMass="320" />
<HexahedronFEMForceField name="FEM" youngModulus="4000" poissonRatio="0.45" method="large" />

<BoxROI name="box" box="-10 -1 -0.0001 -5 4 0.0001"/>
<FixedConstraint indices="@box.indices" />

</Node>
2 changes: 1 addition & 1 deletion src/EigenLinearSolvers/EigenConjugateGradient.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ namespace EigenLinearSolvers
{

/**
* EigenConjugateGradient is a SOFA component allowing to solve a of linear equations. This is an essential component
* EigenConjugateGradient is a SOFA component allowing to solve a set of linear equations. This is an essential component
* in a SOFA simulation.
*
* The class is empty because it is partially specialized in other files.
Expand Down
14 changes: 14 additions & 0 deletions src/EigenLinearSolvers/EigenDirectSparseSolver.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
#pragma once
#include <EigenLinearSolvers/config.h>

namespace EigenLinearSolvers
{

template<class TMatrix, class TVector, class EigenSolver>
class EigenDirectSparseSolver
{
public:
virtual ~EigenDirectSparseSolver() = default;
};

}
83 changes: 83 additions & 0 deletions src/EigenLinearSolvers/EigenDirectSparseSolver[CRS].h
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
#pragma once
#include <EigenLinearSolvers/config.h>

#include <EigenLinearSolvers/EigenDirectSparseSolver.h>

#if __has_include(<sofa/component/linearsolver/iterative/MatrixLinearSolver.h>)
#include <sofa/component/linearsolver/iterative/MatrixLinearSolver.h>
#else
#include <SofaBaseLinearSolver/MatrixLinearSolver.h>
#endif

#include <variant>
#include <Eigen/SparseCore>

#include <sofa/helper/OptionsGroup.h>
#include <EigenLinearSolvers/FindMetis.h>

namespace EigenLinearSolvers
{

/**
* Partial template specialization of EigenDirectSparseSolver for a matrix of type CompressedRowSparseMatrix
*/
template<class TBlockType, class EigenSolver>
class EigenDirectSparseSolver<
sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>,
sofa::linearalgebra::FullVector<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real>,
EigenSolver >
: public sofa::component::linearsolver::MatrixLinearSolver<
sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>,
sofa::linearalgebra::FullVector<typename sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType>::Real> >
{
public:
typedef sofa::linearalgebra::CompressedRowSparseMatrix<TBlockType> Matrix;
using Real = typename Matrix::Real;
typedef sofa::linearalgebra::FullVector<Real> Vector;

SOFA_ABSTRACT_CLASS(SOFA_TEMPLATE3(EigenDirectSparseSolver, Matrix, Vector, EigenSolver),
SOFA_TEMPLATE2(sofa::component::linearsolver::MatrixLinearSolver, Matrix, Vector));

using NaturalOrderSolver = typename EigenSolver::NaturalOrderSolver;
using AMDOrderSolver = typename EigenSolver::AMDOrderSolver;
using COLAMDOrderSolver = typename EigenSolver::COLAMDOrderSolver;
#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1
using MetisOrderSolver = typename EigenSolver::MetisOrderSolver;
#endif

~EigenDirectSparseSolver() override = default;

void init() override;
void reinit() override;

using EigenSparseMatrix = Eigen::SparseMatrix<Real, Eigen::RowMajor>;
using EigenSparseMatrixMap = Eigen::Map<EigenSparseMatrix>;
using EigenVectorXdMap = Eigen::Map<Eigen::Matrix<Real, Eigen::Dynamic, 1> >;

void solve (Matrix& A, Vector& x, Vector& b) override;
void invert(Matrix& A) override;

protected:

sofa::core::objectmodel::Data<sofa::helper::OptionsGroup> d_orderingMethod;
unsigned int m_selectedOrderingMethod { std::numeric_limits<unsigned int>::max() };

std::variant<NaturalOrderSolver, AMDOrderSolver, COLAMDOrderSolver
#if EIGENLINEARSOLVERS_HAS_METIS_INCLUDE == 1
, MetisOrderSolver
#endif
> m_solver;

Eigen::ComputationInfo getSolverInfo() const;
void updateSolverOderingMethod();

sofa::linearalgebra::CompressedRowSparseMatrix<Real> Mfiltered;
std::unique_ptr<EigenSparseMatrixMap> m_map;

typename sofa::linearalgebra::CompressedRowSparseMatrix<Real>::VecIndex MfilteredrowBegin;
typename sofa::linearalgebra::CompressedRowSparseMatrix<Real>::VecIndex MfilteredcolsIndex;

EigenDirectSparseSolver();
};

}
Loading

0 comments on commit 703c2ca

Please sign in to comment.