Multimaterial hyperelasticity using sub-topologies

[alxbilger]

To define regions of an object where materials are different, I thought about using sub-topologies, defined using BoxROI.

<?xml version="1.0" ?>
<Node name="root" dt="0.005" showBoundingTree="0" gravity="0 -9 0">

    <Node name="plugins">
        <RequiredPlugin name="Sofa.Component.AnimationLoop"/> <!-- Needed to use components [FreeMotionAnimationLoop] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Correction"/> <!-- Needed to use components [LinearSolverConstraintCorrection] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Solver"/> <!-- Needed to use components [GenericConstraintSolver] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
        <RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
        <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
        <RequiredPlugin name="Sofa.Component.LinearSystem"/> <!-- Needed to use components [ConstantSparsityPatternSystem] -->
        <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.HyperElastic"/> <!-- Needed to use components [TetrahedronHyperelasticityFEMForceField] -->
        <RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
        <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier] -->
        <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
        <RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2TetraTopologicalMapping] -->
        <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
        <RequiredPlugin name="Sofa.GUI.Component"/> <!-- Needed to use components [ConstraintAttachButtonSetting] -->
    </Node>

    <VisualStyle displayFlags="showForceFields showBehaviorModels" />
    <ConstraintAttachButtonSetting/> <!-- The presence of this component sets the mouse interaction to Lagrangian-based constraints at the GUI launch -->
    <FreeMotionAnimationLoop />
    <GenericConstraintSolver maxIterations="200" tolerance="1.0e-8"/>

    <Node name="hyperelasticity">
        <EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />

        <ConstantSparsityPatternSystem template="CompressedRowSparseMatrix"/>
        <SparseLDLSolver template="CompressedRowSparseMatrix"/>

        <RegularGridTopology name="hexaGrid" min="0 0 0" max="1 1 2.7" n="4 4 11" p0="0 0 0"/>

        <MechanicalObject name="mechObj"/>
        <MeshMatrixMass totalMass="10"/>
        <LinearSolverConstraintCorrection/>

        <Node name="tetras">
            <TetrahedronSetTopologyContainer name="Container" position="@../mechObj.position"/>
            <TetrahedronSetTopologyModifier name="Modifier" />
            <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />
            <Hexa2TetraTopologicalMapping name="mapping" input="@../" output="@Container" />

            <Node name="material1">
                <BoxROI name="materialBox" box="-0.1 -0.1 -0.1  1.1 1.1 1.36" position="@../Container.position" tetrahedra="@../Container.tetrahedra"/>

                <TetrahedronSetTopologyContainer name="Container" position="@../Container.position"
                                                 tetrahedra="@materialBox.tetrahedraInROI"/>
                <TetrahedronSetTopologyModifier name="Modifier" />
                <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />

                <TetrahedronHyperelasticityFEMForceField name="FEM" topology="@Container"
                    ParameterSet="3448.2759 31034.483" materialName="StableNeoHookean"/>

            </Node>

            <Node name="material2">
                <BoxROI name="materialBox" box="-0.1 -0.1 1.34  1.1 1.1 2.8" position="@../Container.position" tetrahedra="@../Container.tetrahedra"/>

                <TetrahedronSetTopologyContainer name="Container" position="@../Container.position"
                                                 tetrahedra="@materialBox.tetrahedraInROI"/>
                <TetrahedronSetTopologyModifier name="Modifier" />
                <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />

                <TetrahedronHyperelasticityFEMForceField name="FEM"
                    ParameterSet="5000 7000 10" materialName="MooneyRivlin"/>
            </Node>

        </Node>

        <BoxROI name="box" drawBoxes="true" box="0 0 0  1 1 0.05"/>
        <FixedProjectiveConstraint indices="@box.indices"/>
    </Node>

</Node>

Visually, the force field looks what is expected, i.e two different regions are displayed:

However it crashes at some point during the simulation.

I also believe that the fact that the force field works on edges is problematic to handle the interface between both materials. I suspect that the edges at the interface is computed twice (once in each material).

[epernod]

Hello Mr @alxbilger ,
there is an engine called : SubsetTopology which is similar to BoxRoi but can handle the renumbering of indices vs the indices in your output edges, triangles, etc..
which might be the cause of your crash.
There are several examples in: https://github.com/sofa-framework/sofa/tree/master/examples/Component/Engine/Select

One is doing something closed to what you want but there are springs in between. I suppose you can avoid that using the option localIndices. I don't remember exactly... I did that in my previous SOFA life.

[alxbilger]

Thanks @epernod. I actually tried to use SubsetTopology, but in the end I did not understand why it was needed so I removed it.

I understand from your comment that you consider both objects independently, and connect them using springs. It's something that I wanted to avoid. Even if we use geometric constraints, I am not sure that the modeling is accurate.

[epernod]

I was thinking again at this. It should definitively work...
Your BBox are overlapping such as the points on the interface are in both topology of the FEM ?
Are you sure TetrahedronHyperelasticityFEMForceField is using the edges? If yes could you try with a more commun TetrahedronFEMForceField with just diffrent your modulus.

[epernod]

Here you are @alxbilger :

<?xml version="1.0" ?>
<Node name="root" dt="0.005" showBoundingTree="0" gravity="0 -9 0">

    <Node name="plugins">
        <RequiredPlugin name="Sofa.Component.AnimationLoop"/> <!-- Needed to use components [FreeMotionAnimationLoop] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Correction"/> <!-- Needed to use components [LinearSolverConstraintCorrection] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Lagrangian.Solver"/> <!-- Needed to use components [GenericConstraintSolver] -->
        <RequiredPlugin name="Sofa.Component.Constraint.Projective"/> <!-- Needed to use components [FixedProjectiveConstraint] -->
        <RequiredPlugin name="Sofa.Component.Engine.Select"/> <!-- Needed to use components [BoxROI] -->
        <RequiredPlugin name="Sofa.Component.LinearSolver.Direct"/> <!-- Needed to use components [SparseLDLSolver] -->
        <RequiredPlugin name="Sofa.Component.LinearSystem"/> <!-- Needed to use components [ConstantSparsityPatternSystem] -->
        <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.HyperElastic"/> <!-- Needed to use components [TetrahedronHyperelasticityFEMForceField] -->
        <RequiredPlugin name="Sofa.Component.StateContainer"/> <!-- Needed to use components [MechanicalObject] -->
        <RequiredPlugin name="Sofa.Component.Topology.Container.Dynamic"/> <!-- Needed to use components [TetrahedronSetGeometryAlgorithms,TetrahedronSetTopologyContainer,TetrahedronSetTopologyModifier] -->
        <RequiredPlugin name="Sofa.Component.Topology.Container.Grid"/> <!-- Needed to use components [RegularGridTopology] -->
        <RequiredPlugin name="Sofa.Component.Topology.Mapping"/> <!-- Needed to use components [Hexa2TetraTopologicalMapping] -->
        <RequiredPlugin name="Sofa.Component.Visual"/> <!-- Needed to use components [VisualStyle] -->
        <RequiredPlugin name="Sofa.GUI.Component"/> <!-- Needed to use components [ConstraintAttachButtonSetting] -->
    </Node>

    <VisualStyle displayFlags="showForceFields showBehaviorModels" />
    <ConstraintAttachButtonSetting/> <!-- The presence of this component sets the mouse interaction to Lagrangian-based constraints at the GUI launch -->
    <FreeMotionAnimationLoop />
    <GenericConstraintSolver maxIterations="200" tolerance="1.0e-8"/>

    <Node name="grid">
        <RegularGridTopology name="hexaGrid" min="0 0 0" max="1 1 2.7" n="4 4 11" p0="0 0 0"/>
        <MechanicalObject name="mechObj"/>
        <Node name="tetras">
            <TetrahedronSetTopologyContainer name="Container" position="@../mechObj.position"/>
            <TetrahedronSetTopologyModifier name="Modifier" />
            <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />
            <Hexa2TetraTopologicalMapping name="mapping" input="@../hexaGrid" output="@Container" />
        </Node>
    </Node>

    <Node name="hyperelasticity">
        <EulerImplicitSolver name="odesolver" rayleighStiffness="0.1" rayleighMass="0.1" />

        <ConstantSparsityPatternSystem template="CompressedRowSparseMatrix"/>
        <SparseLDLSolver template="CompressedRowSparseMatrix"/>

        <TetrahedronSetTopologyContainer name="Container" src="@../grid/tetras/Container"/>
        <TetrahedronSetTopologyModifier name="Modifier" />
        <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />

        <MechanicalObject name="mechObj"/>
        <MeshMatrixMass totalMass="10"/>
        <LinearSolverConstraintCorrection/>

        <Node name="material1">
            <BoxROI name="materialBox" box="-0.1 -0.1 -0.1  1.1 1.1 1.36" position="@../Container.position" tetrahedra="@../Container.tetrahedra"/>

            <TetrahedronSetTopologyContainer name="Container" position="@../Container.position"
                                             tetrahedra="@materialBox.tetrahedraInROI"/>
            <TetrahedronSetTopologyModifier name="Modifier" />
            <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />

            <TetrahedronHyperelasticityFEMForceField name="FEM" topology="@Container"
                ParameterSet="3448.2759 31034.483" materialName="StableNeoHookean"/>

        </Node>

        <Node name="material2">
            <BoxROI name="materialBox" box="-0.1 -0.1 1.34  1.1 1.1 2.8" position="@../Container.position" tetrahedra="@../Container.tetrahedra"/>

            <TetrahedronSetTopologyContainer name="Container" position="@../Container.position"
                                             tetrahedra="@materialBox.tetrahedraInROI"/>
            <TetrahedronSetTopologyModifier name="Modifier" />
            <TetrahedronSetGeometryAlgorithms template="Vec3" name="GeomAlgo" />

            <TetrahedronHyperelasticityFEMForceField name="FEM"
                ParameterSet="5000 7000 10" materialName="MooneyRivlin"/>
        </Node>

        <BoxROI name="box" drawBoxes="true" box="0 0 0  1 1 0.05"/>
        <FixedProjectiveConstraint indices="@box.indices"/>
    </Node>

</Node>

I'm not 100% sure what was the problem in your version but I suspect a conflict between the grid and the topology containers.
In a general way I suggest to always init the grid and topolgy in a another node (as you will do for a loader) from the simulated model.