Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[all] Sofa add mechanical matrix mapper #721

Merged
Merged
Show file tree
Hide file tree
Changes from 15 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 0 additions & 3 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -127,9 +127,6 @@ option(SOFA_WITH_DEVTOOLS "Compile with developement extra features." ON)
### Threading
option(SOFA_WITH_THREADING "Compile sofa with thread-safetiness support (PARTIAL/EXPERIMENTAL)" ON)

### Experimental
option(SOFA_WITH_EXPERIMENTAL_FEATURES "Compile sofa with exprimental features regarding sparse matrices treatments under mappings" OFF)

### Testing
option(SOFA_BUILD_TESTS "Compile the automatic tests for Sofa, along with the gtest library." ON)
if(SOFA_BUILD_TESTS)
Expand Down
7 changes: 0 additions & 7 deletions SofaKernel/SofaFramework/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -42,12 +42,6 @@ else()
set(SOFA_WITH_THREADING_ "0")
endif()

if(SOFA_WITH_EXPERIMENTAL_FEATURES)
set(SOFA_WITH_EXPERIMENTAL_FEATURES_ "1")
else()
set(SOFA_WITH_EXPERIMENTAL_FEATURES_ "0")
endif()

if(SOFA_WITH_DEPRECATED_COMPONENTS)
set(SOFA_WITH_DEPRECATED_COMPONENTS_ "1")
else()
Expand Down Expand Up @@ -175,7 +169,6 @@ set(SOFA_BUILD_OPTIONS_SRC
sharedlibrary_defines.h
build_option_deprecated_components.h
build_option_dump_visitor.h
build_option_experimental_features.h
build_option_opengl.h
build_option_threading.h
)
Expand Down
6 changes: 3 additions & 3 deletions SofaKernel/framework/sofa/core/behavior/BaseMechanicalState.h
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@
#ifndef SOFA_CORE_BEHAVIOR_BASEMECHANICALSTATE_H
#define SOFA_CORE_BEHAVIOR_BASEMECHANICALSTATE_H

#include <sofa/config/build_option_experimental_features.h>

#include <sofa/core/BaseState.h>
#include <sofa/core/MultiVecId.h>
Expand Down Expand Up @@ -222,11 +221,12 @@ class SOFA_CORE_API BaseMechanicalState : public virtual BaseState

/// build the jacobian of the constraint in a baseMatrix
virtual void getConstraintJacobian(const ConstraintParams* params, sofa::defaulttype::BaseMatrix* J,unsigned int & off) = 0;
#if SOFA_WITH_EXPERIMENTAL_FEATURES() == 1

/// fill the jacobian matrix (of the constraints) with identity blocks on the provided list of nodes(dofs)
virtual void buildIdentityBlocksInJacobian(const sofa::helper::vector<unsigned int>& list_n, core::MatrixDerivId &mID) = 0;
#endif

/// Renumber the constraint ids with the given permutation vector
virtual void renumberConstraintId(const sofa::helper::vector<unsigned>& renumbering) = 0;

class ConstraintBlock
{
Expand Down
40 changes: 40 additions & 0 deletions SofaKernel/framework/sofa/simulation/MechanicalVisitor.h
Original file line number Diff line number Diff line change
Expand Up @@ -2077,6 +2077,46 @@ class SOFA_SIMULATION_CORE_API MechanicalAccumulateMatrixDeriv : public BaseMech
bool reverseOrder;
};

class SOFA_SIMULATION_CORE_API MechanicalRenumberConstraint : public MechanicalVisitor
{
public:
MechanicalRenumberConstraint(const sofa::core::MechanicalParams* mparams , const sofa::helper::vector<unsigned> &renumbering)
: MechanicalVisitor(mparams) , renumbering(renumbering)
{
#ifdef SOFA_DUMP_VISITOR_INFO
setReadWriteVectors();
#endif
}
virtual Result fwdMechanicalState(simulation::Node* /*node*/, core::behavior::BaseMechanicalState* mm)
{
mm->renumberConstraintId(renumbering);
return RESULT_PRUNE;
}
// This visitor must go through all mechanical mappings, even if isMechanical flag is disabled
virtual bool stopAtMechanicalMapping(simulation::Node* /*node*/, core::BaseMapping* /*map*/)
{
return false; // !map->isMechanical();
}


/// Return a class name for this visitor
/// Only used for debugging / profiling purposes
virtual const char* getClassName() const { return "MechanicalRenumberConstraint"; }

virtual bool isThreadSafe() const
{
return false;
}
#ifdef SOFA_DUMP_VISITOR_INFO
void setReadWriteVectors()
{
}
#endif

protected:
const sofa::helper::vector<unsigned> &renumbering;
};

/** Apply the constraints as filters to the given vector.
This works for simple independent constraints, like maintaining a fixed point.
*/
Expand Down
9 changes: 3 additions & 6 deletions SofaKernel/modules/SofaBaseMechanics/MechanicalObject.h
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,6 @@
#define SOFA_COMPONENT_MECHANICALOBJECT_H
#include "config.h"

#include <sofa/config/build_option_experimental_features.h>

#include <sofa/core/behavior/MechanicalState.h>
#include <sofa/core/topology/BaseMeshTopology.h>

Expand Down Expand Up @@ -147,9 +145,7 @@ class MechanicalObject : public sofa::core::behavior::MechanicalState<DataTypes>
Data< VecDeriv > vfree; ///< free velocity coordinates of the degrees of freedom
Data< VecCoord > x0; ///< rest position coordinates of the degrees of freedom
Data< MatrixDeriv > c; ///< constraints applied to the degrees of freedom
#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)
Data< MatrixDeriv > m; ///< mappingJacobian applied to the degrees of freedom
#endif
Data< VecCoord > reset_position; ///< reset position coordinates of the degrees of freedom
Data< VecDeriv > reset_velocity; ///< reset velocity coordinates of the degrees of freedom
#endif
Expand Down Expand Up @@ -308,6 +304,8 @@ class MechanicalObject : public sofa::core::behavior::MechanicalState<DataTypes>

/// @}

/// Renumber the constraint ids with the given permutation vector
void renumberConstraintId(const sofa::helper::vector< unsigned >& renumbering) override;


/// @name Integration related methods
Expand Down Expand Up @@ -376,9 +374,8 @@ class MechanicalObject : public sofa::core::behavior::MechanicalState<DataTypes>
virtual void resetConstraint(const core::ConstraintParams* cparams) override;

virtual void getConstraintJacobian(const core::ConstraintParams* cparams, sofa::defaulttype::BaseMatrix* J,unsigned int & off) override;
#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)

virtual void buildIdentityBlocksInJacobian(const sofa::helper::vector<unsigned int>& list_n, core::MatrixDerivId &mID) override;
#endif
/// @}

/// @name Debug
Expand Down
16 changes: 5 additions & 11 deletions SofaKernel/modules/SofaBaseMechanics/MechanicalObject.inl
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@
#ifndef SOFA_COMPONENT_MECHANICALOBJECT_INL
#define SOFA_COMPONENT_MECHANICALOBJECT_INL

#include <sofa/config/build_option_experimental_features.h>

#include <SofaBaseMechanics/MechanicalObject.h>
#include <sofa/core/visual/VisualParams.h>
#include <SofaBaseLinearSolver/SparseMatrix.h>
Expand Down Expand Up @@ -87,9 +85,7 @@ MechanicalObject<DataTypes>::MechanicalObject()
, vfree(initData(&vfree, "free_velocity", "free velocity coordinates of the degrees of freedom"))
, x0(initData(&x0, "rest_position", "rest position coordinates of the degrees of freedom"))
, c(initData(&c, "constraint", "constraints applied to the degrees of freedom"))
#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)
, m(initData(&m, "mappingJacobian", "mappingJacobian applied to the degrees of freedom"))
#endif
, reset_position(initData(&reset_position, "reset_position", "reset position coordinates of the degrees of freedom"))
, reset_velocity(initData(&reset_velocity, "reset_velocity", "reset velocity coordinates of the degrees of freedom"))
, restScale(initData(&restScale, (SReal)1.0, "restScale", "optional scaling of rest position coordinates (to simulated pre-existing internal tension).(default = 1.0)"))
Expand Down Expand Up @@ -149,9 +145,7 @@ MechanicalObject<DataTypes>::MechanicalObject()
setVecDeriv(core::VecDerivId::freeVelocity().index, &vfree);
setVecDeriv(core::VecDerivId::resetVelocity().index, &reset_velocity);
setVecMatrixDeriv(core::MatrixDerivId::constraintJacobian().index, &c);
#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)
setVecMatrixDeriv(core::MatrixDerivId::mappingJacobian().index, &m);
#endif

// These vectors are set as modified as they are mandatory in the MechanicalObject.
x .forceSet();
Expand Down Expand Up @@ -2544,12 +2538,10 @@ void MechanicalObject<DataTypes>::resetConstraint(const core::ConstraintParams*
MatrixDeriv *c = c_data.beginEdit(cParams);
c->clear();
c_data.endEdit(cParams);
#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)
Data<MatrixDeriv>& m_data = *this->write(core::MatrixDerivId::mappingJacobian());
MatrixDeriv *m = m_data.beginEdit(cParams);
m->clear();
m_data.endEdit(cParams);
#endif
}

template <class DataTypes>
Expand Down Expand Up @@ -2580,7 +2572,6 @@ void MechanicalObject<DataTypes>::getConstraintJacobian(const core::ConstraintPa
off += this->getSize() * N;
}

#if(SOFA_WITH_EXPERIMENTAL_FEATURES()==1)
template <class DataTypes>
void MechanicalObject<DataTypes>::buildIdentityBlocksInJacobian(const sofa::helper::vector<unsigned int>& list_n, core::MatrixDerivId &mID)
{
Expand All @@ -2607,9 +2598,12 @@ void MechanicalObject<DataTypes>::buildIdentityBlocksInJacobian(const sofa::help
cMatrix->endEdit();

}
#endif


template <class DataTypes>
void MechanicalObject<DataTypes>::renumberConstraintId(const sofa::helper::vector<unsigned>& /*renumbering*/)
{
this->serr << "MechanicalObject<DataTypes>::renumberConstraintId not implemented in the MatrixDeriv constraint API" << this->sendl;
Copy link
Contributor

@damienmarchal damienmarchal Jul 25, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please don't use serr/sout anymore...
but more fundamentally... a pure virtual method is ...pure specifically to force the object inheriting from the parent to not forget to implement them. Providing an empty implementation like this one sounds like skipping over the safeguard that was put in place so it compile and let the user handle the problem (which I doubt it can).

Don't we have a better way to handle those cases ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the override is mandatory, a pure virtual function is the way to go. This compiler-bypassing behavior should not be allowed.

If the override is not mandatory, then non-virtual + msg_warning("You didn't override") seems nice.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with your Guillaume :)

}

template <class DataTypes>
std::list< core::behavior::BaseMechanicalState::ConstraintBlock > MechanicalObject<DataTypes>::constraintBlocks( const std::list<unsigned int> &indices) const
Expand Down
26 changes: 26 additions & 0 deletions examples/Components/animationloop/MechanicalMatrixMapper.html
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
<html>


<body bgcolor="#F0F8FF">

<div id="contenu">

<h1>Mechanical matrix mapper</h1>

<p>In this scene we demonstrate how to use the <FONT COLOR="#0B2161">MechanicalMatrixMapper</FONT> component. This component allows to map mechanical matrices (Stiffness, Mass) through a mapping.
<p>This is needed in SOFA scenes having these two following particularities:
<p> -- There are using a direct solver (e.g. SparseLDLSolver) that, unlike
iterative solvers, need to build the mechanical matrices.
<p> -- They involves ForceFields that implement addKToMatrix (i.e. that compute internal forces such as e.g. TetrahedronFEMForceField,
TetrahedronHyperElasticityFEMForceField, but not ConstantForceField which only contributes to the right-hand side) and that
ARE USED UNDER mappings. In this scene, the <FONT COLOR="#0B2161">HexaHedronFEMForceField</FONT> is used under a <FONT COLOR="#0B2161">SubsetMultiMapping</FONT> and a <FONT COLOR="#0B2161">RigidMapping</FONT>), and the MechanicalMatrixMapper is needed to map that ForceField to the degrees of freedom closer to the root in the scene graph.

<p>Without this component, such a scene either crashes or gives unlogical behaviour.

<br>

</div>

</body>

</html>
48 changes: 48 additions & 0 deletions examples/Components/animationloop/MechanicalMatrixMapper.pyscn
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
import Sofa

def createScene(rootNode):

rootNode.createObject('RequiredPlugin', name='ModelOrderReduction')
rootNode.createObject('VisualStyle', displayFlags='showCollisionModels hideMappings showForceFields')

meshDivision = rootNode.createChild('meshDivision')
meshDivision.createObject('RegularGridTopology', name='topology', n=[2, 2, 4] , min=[-1, -1, -2], max=[1, 1, 2])
meshOfStructure = meshDivision.createObject('Mesh',name='beamMesh', src="@topology")
meshDivision.createObject('MechanicalObject', template='Vec3d')
boxFixed= meshDivision.createObject('BoxROI',template="Vec3d", name='box_roi_fix',box=[-1, -1, -2.05, 1, 1, -1.5] , position="@beamMesh.position" )
boxDeformable= meshDivision.createObject('BoxROI',template="Vec3d", name='box_roi_deformable',box=[-1, -1, -2.05, 1, 1, 1.5], position="@beamMesh.position" )
boxRigid= meshDivision.createObject('BoxROI',template="Vec3d", name='box_roi_rigid',box=[-1, -1, 1.500001, 1, 1, 2.05], position="@beamMesh.position")

SolverNode= rootNode.createChild('SolverNode')
SolverNode.createObject('EulerImplicit',verbose='false')
SolverNode.createObject('SparseLDLSolver', name="ldlsolveur")
SolverNode.createObject('MechanicalMatrixMapper', template='Vec3d,Rigid', object1='@./deformablePartNode/beamPart1Mech', object2='@./RigidNode/rigid1', nodeToParse='@./deformablePartNode/FEMNode')


##### Deformable Part of the BEAM (Main Body)
deformablePartNode= SolverNode.createChild('deformablePartNode')
deformablePartNode.createObject('MechanicalObject', template='Vec3d',name='beamPart1Mech', position="@"+boxDeformable.getPathName()+".pointsInROI")

##### Rigid Part of the BEAM (Top)
RigidNode= SolverNode.createChild('RigidNode')
RigidNode.createObject("MechanicalObject",template="Rigid",name="rigid1", position=[0, 0, 2, 0, 0, 0, 1], showObject=True, showObjectScale=0.5)
RigidifiedNode= RigidNode.createChild('RigidifiedNode')

RigidifiedNode.createObject("MechanicalObject",name="rigidMecha",template="Vec3d", position="@"+boxRigid.getPathName()+".pointsInROI")
RigidifiedNode.createObject("RigidMapping",globalToLocalCoords="true")

##### COMBINED BEAM
FEMNode= deformablePartNode.createChild('FEMNode')
RigidifiedNode.addChild(FEMNode)

FEMNode.createObject('Mesh',name='meshInput',src="@"+meshOfStructure.getPathName())
FEMNode.createObject('MechanicalObject', template='Vec3d',name='beamMecha')
FEMNode.createObject('HexahedronFEMForceField', name='HexaFF', src="@meshInput", poissonRatio=0.49, youngModulus=2000)
FEMNode.createObject('RestShapeSpringsForceField', name='restShapeFF', points="@"+boxFixed.getPathName()+'.indices', stiffness=10000, angularStiffness=10000)
FEMNode.createObject('UniformMass', totalmass=0.1)
FEMNode.createObject('ConstantForceField', name='xMoins', indices=15, forces=[1000, 0, 0])
FEMNode.createObject("SubsetMultiMapping",name="subsetMapping",template="Vec3d,Vec3d", input="@../beamPart1Mech @../../RigidNode/RigidifiedNode/rigidMecha",output="@./beamMecha", indexPairs="0 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 10 0 11 1 0 1 1 1 2 1 3")

return rootNode


9 changes: 6 additions & 3 deletions modules/SofaGeneralAnimationLoop/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,17 @@ set(SOURCE_FILES

list(APPEND HEADER_FILES
MultiStepAnimationLoop.h
MultiTagAnimationLoop.h )
MultiTagAnimationLoop.h
MechanicalMatrixMapper.h
MechanicalMatrixMapper.inl)

list(APPEND SOURCE_FILES
MultiStepAnimationLoop.cpp
MultiTagAnimationLoop.cpp )
MultiTagAnimationLoop.cpp
MechanicalMatrixMapper.cpp)

add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES})
target_link_libraries(${PROJECT_NAME} PUBLIC SofaSimulationCommon)
target_link_libraries(${PROJECT_NAME} PUBLIC SofaSimulationCommon SofaBaseLinearSolver)
set_target_properties(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-DSOFA_BUILD_GENERAL_ANIMATION_LOOP")
set_target_properties(${PROJECT_NAME} PROPERTIES PUBLIC_HEADER "${HEADER_FILES}")

Expand Down
77 changes: 77 additions & 0 deletions modules/SofaGeneralAnimationLoop/MechanicalMatrixMapper.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@
/******************************************************************************
* SOFA, Simulation Open-Framework Architecture, version 1.0 RC 1 *
* (c) 2006-2018 MGH, INRIA, USTL, UJF, CNRS *
* *
* This library is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by *
* the Free Software Foundation; either version 2.1 of the License, or (at *
* your option) any later version. *
* *
* This library is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this library; if not, write to the Free Software Foundation, *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
*******************************************************************************
* SOFA :: Modules *
* *
* Authors: The SOFA Team and external contributors (see Authors.txt) *
* *
* Contact information: contact@sofa-framework.org *
******************************************************************************/
#include "MechanicalMatrixMapper.inl"
#include <sofa/defaulttype/Vec3Types.h>
#include <sofa/defaulttype/RigidTypes.h>
#include <sofa/core/ObjectFactory.h>

namespace sofa
{

namespace component
{

namespace interactionforcefield
{

using namespace sofa::defaulttype;

//////////////////////////////////////////// FACTORY //////////////////////////////////////////////
SOFA_DECL_CLASS(MechanicalMatrixMapper)

int MechanicalMatrixMapperClass = core::RegisterObject("This component allows to map the stiffness (and mass) matrix through a mapping.")
#ifdef SOFA_WITH_FLOAT
.add< MechanicalMatrixMapper<Vec3fTypes, Rigid3fTypes> >()
.add< MechanicalMatrixMapper<Vec3fTypes, Vec3fTypes> >()
.add< MechanicalMatrixMapper<Vec1fTypes, Rigid3fTypes> >()
.add< MechanicalMatrixMapper<Vec1fTypes, Vec1fTypes> >()
#endif
#ifdef SOFA_WITH_DOUBLE
.add< MechanicalMatrixMapper<Vec3dTypes, Rigid3dTypes> >(true)
.add< MechanicalMatrixMapper<Vec3dTypes, Vec3dTypes> >(true)
.add< MechanicalMatrixMapper<Vec1dTypes, Rigid3dTypes> >(true)
.add< MechanicalMatrixMapper<Vec1dTypes, Vec1dTypes> >(true)
#endif
;
////////////////////////////////////////////////////////////////////////////////////////////////////////

#ifdef SOFA_WITH_DOUBLE
template class MechanicalMatrixMapper<Vec3dTypes, Rigid3dTypes>;
template class MechanicalMatrixMapper<Vec3dTypes, Vec3dTypes>;
template class MechanicalMatrixMapper<Vec1dTypes, Rigid3dTypes>;
template class MechanicalMatrixMapper<Vec1dTypes, Vec1dTypes>;
#endif
#ifdef SOFA_WITH_FLOAT
template class MechanicalMatrixMapper<Vec3fTypes, Rigid3fTypes>;
template class MechanicalMatrixMapper<Vec3fTypes, Vec3fTypes>;
template class MechanicalMatrixMapper<Vec1fTypes, Rigid3fTypes>;
template class MechanicalMatrixMapper<Vec1fTypes, Vec1fTypes>;
#endif

} // namespace forcefield

} // namespace component

} // namespace sofa
Loading