Error using Triangle topo with RegularGrid

[hugtalbot]

Using a RegularGrid in 2D, a TriangleSetTopologyContainer returns the following error:

Cannot find edge 1 [12, 0] in triangle 0

and the structure edges in triangles is not created. I guess this would be due to the re-ordering of elements (hexa/quads/tri). Any idea @epernod ?

Example scene:

from stlib.scene import MainHeader
import Sofa
import math
from time import sleep
from splib.animation import animate, AnimationManager

def createScene(rootNode):
    rootNode.gravity = [0,0,0]
    rootNode.dt = 0.001
    rootNode.createObject("VisualStyle",  displayFlags="hideVisual showBehaviorModels showForceFields")
    rootNode.createObject("RequiredPlugin", name="SofaPython")

    returnMembrane = Membrane(rootNode, [])
    return 0;

class Membrane(Sofa.PythonScriptController):
    def createGraph(self, rootNode):
        return;

    def initGraph(self, rootNode):  
        self.root = rootNode
        #-----------------------------------------------------------------------------------------------------------------------------------------
        # Properties of the membrane
        self.MembraneLength = 1.0 # in mm
        self.YoungModSI = 6000.0 # in Pascal (SI)
        self.StiffnessSI = 300.0 # Force per unit length
        self.scale = 1 # in mm

        # Properties of the laser
        self.EXModulationFreq = 1.0 # in Hz
        self.EXOmega  = 2*math.pi*self.EXModulationFreq
        self.forceAmpl = 10.0
        self.LaserspotDiameter = 0.05

        # Simulation options
        self.root.findData("gravity").value = [0,0,0]
        self.positionCenter = [0,0,0] 
        self.DrivingForce =  [0,0,0]
        self.IterationIndices = [0,0]

        # Measurement options
        self.nbrOfPoints = 10.0
        self.LaserSpacing = self.MembraneLength/self.nbrOfPoints
        self.TotalPoints = self.nbrOfPoints*self.nbrOfPoints
        #------------------------------------------------------------------------------------------------------------------------------------------
        #------------------------------------------------------------------------------------------------------------------------------------------
        # Define square node and all corresponding properties
        square = rootNode.createChild("SquareGravity")

        square.createObject("EulerImplicit", name="eulerImplcit_Scheme" , printLog="false", rayleighStiffness="0.0", rayleighMass="0.0")
        square.createObject("CGLinearSolver", name="cgLinearSolver", iterations="200", tolerance="1.0e-9", threshold="1.0e-9")

        square.createObject("RegularGrid", name="RegularGrid", nx="11", ny="11", nz="1", xmin="0", xmax=str(self.MembraneLength), ymin="0", ymax=str(self.MembraneLength), zmin="0", zmax="0", drawTriangles="0")

        #square.createObject("MeshGmshLoader", name="GmshLoader", filename="mesh/square3.msh", scale="1", createSubelements="true")

        #square.createObject("RegularGridTopology", name="grid", n="6 6 1", min="0 0 0", max="1 1 1", p0="0 1 0", computeHexaList="0")

        square.createObject("TriangleSetTopologyContainer", name="TopologyContainer", src = "@RegularGrid")
        square.createObject("TriangleSetTopologyAlgorithms", name="TopologyAlgo", src = "@RegularGrid")
        # square.createObject("TriangleSetTopologyContainer", name="TopologyContainer", src = "@GmshLoader")
        square.createObject("TriangleSetTopologyModifier", name="TopologyModifier")
        square.createObject("TriangleSetGeometryAlgorithms", name="TopologyAlgorithm", template="Vec3d")

        self.MO = square.createObject("MechanicalObject")

        square.createObject("DiagonalMass", massDensity="0.15")

        square.createObject("BoxROI", name="boxROI1", box="-0.01 -0.01 -0.01 1.01 0.01 0.01", drawBoxes="1")
        square.createObject("FixedConstraint", indices="@boxROI1.indices" )
        square.createObject("BoxROI", name="boxROI2", box="-0.01 -0.01 -0.01 0.01 1.01 0.01", drawBoxes="1")
        square.createObject("FixedConstraint", indices="@boxROI2.indices")
        square.createObject("BoxROI", name="boxROI3", box="-0.01 0.99 -0.01 1.01 1.01 0.01", drawBoxes="1")
        square.createObject("FixedConstraint", indices="@boxROI3.indices")
        square.createObject("BoxROI", name="boxROI4", box="0.99 -0.01 -0.01 1.01 1.01 0.01", drawBoxes="1")
        square.createObject("FixedConstraint", indices="@boxROI4.indices")

        square.createObject("TriangularFEMForceField", name="FEM", youngModulus=str(self.YoungModSI/self.scale**2), poissonRatio="0.3", method="large")
        # square.createObject("TriangularBendingSprings", name="FEM-Bend", stiffness=str(self.StiffnessSI/self.scale), damping="0.0")
        square.createObject("DiagonalVelocityDampingForceField", name="FEM-Damp", dampingCoefficient="10")
        self.membrane = square

        self.CFF = square.createObject("ConstantForceField", name="FEM-CFF", force="0 0 0")

        # Define not modulated plate
        rest = rootNode.createChild("RestSquare")

        rest.createObject("RegularGrid", name="RG", nx="11", ny="11", nz="1", xmin="0", xmax=str(self.MembraneLength), ymin="0", ymax=str(self.MembraneLength), zmin="0", zmax=str(self.MembraneLength))
        #rest.createObject("MeshGmshLoader", name="GmshLoader", filename="mesh/square3.msh", scale="1", createSubelements="true" )

        #rest.createObject("TriangleSetTopologyContainer", name="TopologyContainer", src = "@GmshLoader")
        rest.createObject("TriangleSetTopologyContainer", name="TopologyContainer", src = "@RG")
        rest.createObject("TriangleSetGeometryAlgorithms", name="TopologyAlgorithm", template="Vec3d")

        rest.createObject("MechanicalObject", name="REST")      
# 
        # Define sphere which represents the laser
        self.ROI = rest.createObject("SphereROI", name="boxROI1", centers="0 0 0", radii=str(self.LaserspotDiameter), listening="1", drawSphere="1")

        return(rootNode)

---

Suggested labels:

[epernod]

Back to this issue. I dig into the code. The problem comes from the fact that the RegularGridTopology is computing a Quad topology and then creating the triangles by splitting the quads.

For example for a grid of nx="3" ny="3" nz="1" There will be :

• 9 positions
• 12 edges
• 8 triangles
• 4 quads

Instead of:

• 9 positions
• 16 edges
• 8 triangles

So using <TriangleSetTopologyContainer src="@grid" /> will force the topology to init itself on non-coherent buffers. If you use: <TriangleSetTopologyContainer triangles="@grid.triangles" position="@grid.position" /> it will recompute the good number of edges and the scene will works.

Thus I see several ways to solve the issue:

1. Easy way but not fixing the issue: just change your scene.
2. Change the code to always recompute the edges when using a TriangleSetTopologyContainer (or Quad, Tetra, Hexa). this add some computation but it ensures the topology will be correct (I'm still thinking if there is some cases where you wouldn't that computation)
3. Add an option in the TriangleSetTopologyContainer to force computation of edges. Similar to point 2. but only on demands ... Not very fan of this fix.

[epernod]

Solution chose: change Regular grid to compute: triangles + consistent edges An an option to compute only quads + edges (on demand)