[SOLVED] “Could not read value for data field position” and segfault

Hello,

I have encountered an issue that might be a bug but I prefer to ask before opening an issue.
I have a scene in with bodies discredited with beam elements enter into contact.

When I increase the number of nodes used to discretise each body, I have a series of warnings that look like:
[WARNING] [MechanicalObject(DOFs)] Could not read value for data field position: 0. 0. 0. ..., 0. 0. 0.
After this, I also get:

[WARNING] [FixedConstraint(FixedConstraint)] Index 0 not valid, should be [0,0]. Constraint will be removed.

########## SIG 11 - SIGSEGV: segfault ##########
Segmentation fault: 11

The code I am running is the following (sorry for the length):

import sys
import os
import Sofa
import numpy as np
from pyquaternion import Quaternion
import pickle 

def save_obj(obj, name ):
    with open(name + '.pkl', 'wb') as f:
        pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)

def load_obj(name ):
    with open(name + '.pkl', 'rb') as f:
        return pickle.load(f)

def getOrientationQuaternion(a, b):
    assert np.allclose(np.linalg.norm(b), 1.)
    v = np.cross(a,b)
    s = np.linalg.norm(v) 
    if not np.allclose(s, 0):
        c = np.dot(a,b)
        vskew = np.array([[0, -v[2], v[1]], [v[2], 0 , -v[0]], [-v[1], v[0], 0]  ])
        # rotation matrix rotating a into b
        R = np.eye(3) + vskew + np.dot(vskew, vskew) * ((1-c)/(s**2))
    else:
        R = np.eye(3)
    assert np.allclose(R.dot(a), b)
    # components are not in the right order 
    qi = Quaternion(matrix=R).elements
    return np.array([qi[1], qi[2], qi[3], qi[0]])

def getCoordPointsAlongHellix(r, h, npt, tmax, c = 1):
    x = np.zeros((npt,3), dtype = float)
    # field of tangent vectors
    xp = np.zeros((npt,3), dtype = float)
    t = np.linspace(0,tmax, npt)
    x[:,0] = r * np.cos(t)
    x[:,1] = r * np.sin(t)
    x[:,2] = c * h * t

    xp[:,0] = -r * np.sin(t)
    xp[:,1] = r * np.cos(t)
    xp[:,2] = c * h 

    # https://math.stackexchange.com/a/476311/392320
    q = np.zeros((npt,4), dtype = float)
    a = np.array([1,0,0])
    # from this link, it looks like we need to provide the unit tangent vectors to the curve
    for i in range(npt):
        #tangent to the curve at the node
        b = xp[i]/np.linalg.norm(xp[i])
        q[i] = getOrientationQuaternion(a, b) 

        return x, q

class BeamFEMForceField (Sofa.PythonScriptController):

    def __init__(self, node, commandLineArguments) : 
        self.commandLineArguments = commandLineArguments
        print "Command line arguments for python : "+str(commandLineArguments)
        self.createGraph(node)
        return None;

    def createGraph(self,rootNode):
        UL = 1E-3
        # pitch length
        p = 115 * UL
        # heigth of the strand 
        h = 2 * p 
        # radius of the core wire 
        # found in the paper of Jiang 1999
        # careful because the diameter is given, not the radius
        rCore = 0.5 * 3.93 * UL
        rHelli = 0.5 * 3.73 * UL
        radiusSphereCollisionModelCore = str(0.95 * rCore )
        radiusSphereCollisionModelHeli = str(0.95 * rHelli)
        # Young's modulus
        E=188e9
        # Poisson's ratio
        nu = 0.3
        # number of nodes
        nn = int(self.commandLineArguments[1])
        assert nn > 1
        nbeam = nn - 1
        # target strand strain
        epsilon = 0.01 
        # take gravity into account
        gravity = 0
        BC = ['AxialForceAtTip', 'MoveNodeAtTip'][1]
        indexNodeBC = str(nn-1)
        DirectoryResults='./Tensile/'
        self.DirectoryResults = DirectoryResults
        if not os.path.exists(DirectoryResults):
            os.mkdir(DirectoryResults)

        self.rootNode = rootNode
        # rootNode
        rootNode.createObject('RequiredPlugin', name='SofaOpenglVisual')
        rootNode.createObject('RequiredPlugin', name='SofaPython')
        showCollisionModels = 1
        if showCollisionModels:
            rootNode.createObject('VisualStyle', displayFlags='showBehaviorModels showForceFields showCollisionModels')
        else:
            rootNode.createObject('VisualStyle', displayFlags='showBehaviorModels showForceFields')
        rootNode.createObject('DefaultPipeline', draw='0', depth='6', verbose='0')
        rootNode.createObject('BruteForceDetection', name='N2')
        rootNode.createObject('MinProximityIntersection', contactDistance= 0.005 * rCore , alarmDistance= 0.01 * rCore, name='Proximity')
        rootNode.createObject('DefaultContactManager', name='Response', response='default')
        if gravity:
            rootNode.gravity = [0,0,-9.81]
        else:
            rootNode.gravity = [0,0,0]

        rootNode.dt = 0.001
        rootNode.createObject('EulerImplicitSolver', rayleighStiffness='0', printLog='false', rayleighMass='0.1')
        rootNode.createObject('CGLinearSolver', threshold='1.0e-9', tolerance='1.0e-9', name='linear solver', iterations='25')

        # topology
        lines = np.zeros((nbeam, 2), dtype=int) 
        lines[:,0] = np.arange(0, nn -1)
        lines[:,1] = np.arange(1, nn)
        lines = str(lines.flatten()).replace('[', '').replace(']','')

        # rootNode/CentralBeam --> straight beam
        x, q = getCoordPointsAlongHellix(0, h, nn, tmax = 2 * np.pi, c = p)
        number_dofs_per_node = 7
        Coord = np.zeros((nn, number_dofs_per_node),dtype=float)
        Coord[:, :3] = x 
        Coord[:, 3:] = q
        lengthStrand = Coord[-1, 2]
        d = {"lengthStrand": lengthStrand}
        save_obj(d, DirectoryResults + 'infoSimu' )

        if BC == 'AxialForceAtTip':
            #index of the last node
            forceApplied = [0, 0, 10000, 0, 0, 0]
        elif BC == 'MoveNodeAtTip':
            index = str(nn-1)
            # easy to retrieve the displacement
            disp = epsilon * lengthStrand 
            keyTimes = np.zeros(3)
            keyTimes[0] = 0 
            keyTimes[1] = 10
            keyTimes[2] = 20
            movements = np.zeros((keyTimes.shape[0], 6), dtype=float)
            movements[1] = [0, 0, disp, 0, 0, 0]
            movements[2] = [0, 0, disp, 0, 0, 0]
            keyTimes = keyTimes.ravel().tolist()
            movements = movements.ravel().tolist()

        strCoord =  str(Coord.flatten()).replace('\n', '').replace('[', '').replace(']','')
        CentralBeam = rootNode.createChild('CentralBeam')
        self.CentralBeam = CentralBeam
        CentralBeam.createObject('MeshTopology', lines=lines, name='lines')
        CentralBeam.createObject('MechanicalObject', position=strCoord, name='DOFs', template='Rigid3d', showObject=0)
        CentralBeam.createObject('FixedConstraint', indices='0', name='FixedConstraint')
        CentralBeam.createObject('UniformMass', vertexMass='1 1 0.01 0 0 0 0.1 0 0 0 0.1', printLog='false', showAxisSizeFactor=0)
        CentralBeam.createObject('BeamFEMForceField', radius=str(rCore), name='FEM', poissonRatio=nu, youngModulus=E)
        Collision = CentralBeam.createChild('Collision')
        self.Collision = Collision
        CentralBeam.createObject('SphereModel', radius=radiusSphereCollisionModelCore, name='SphereCollision')
        if BC == 'AxialForceAtTip':
            CentralBeam.createObject('ConstantForceField', indices=indexNodeBC, showArrowSize='0.0005', printLog='1', forces=forceApplied)
        elif BC == 'MoveNodeAtTip':
            CentralBeam.createObject('LinearMovementConstraint', keyTimes=keyTimes, template='Rigid3d', movements=movements, indices=indexNodeBC)

        # rootNode/beamJ --> helicoidal beam
        rotAngle = 360/6
        x, q = getCoordPointsAlongHellix( (rCore+rHelli), h, nn, tmax = 2 * np.pi, c = p)
        Coord[:,:3] = x 
        Coord[:, 3:] = q
        strCoord =  str(Coord.flatten()).replace('\n', '').replace('[', '').replace(']','')
        for i in range(6):
            strbeam = 'HellicalBeam' + str(i) 
            HeliBeami = rootNode.createChild(strbeam)
            setattr(self, strbeam, HeliBeami)
            # HeliBeami.createObject('EulerImplicitSolver', rayleighStiffness='0', printLog='false', rayleighMass='0.1')
            # HeliBeami.createObject('BTDLinearSolver', printLog='false', template='BTDMatrix6d', verbose='false')
            rot = "0 0 " + str(i * rotAngle)
            HeliBeami.createObject('MechanicalObject', position=strCoord, name='DOFs', template='Rigid3d', rotation=rot)
            HeliBeami.createObject('MeshTopology', lines=lines, name='lines')
            HeliBeami.createObject('FixedConstraint', indices='0', name='FixedConstraint')
            HeliBeami.createObject('UniformMass', vertexMass='1 1 0.01 0 0 0 0.1 0 0 0 0.1', printLog='false',  showAxisSizeFactor=0)
            HeliBeami.createObject('BeamFEMForceField', radius=str(rHelli), name='FEM', poissonRatio=nu, youngModulus=E)
            HeliBeami.createObject('SphereModel', radius=radiusSphereCollisionModelHeli, name='SphereCollision')
            index = str(nn-1)
            if BC == 'AxialForceAtTip':
                HeliBeami.createObject('ConstantForceField', indices=indexNodeBC, showArrowSize='0.0005', printLog='1', forces=forceApplied)
            elif BC == 'MoveNodeAtTip':
                HeliBeami.createObject('LinearMovementConstraint', keyTimes=keyTimes, template='Rigid3d', movements=movements, indices=indexNodeBC)
            HeliBeami.createObject('Monitor', name="ReactionForce" + strbeam, indices=0,
            template="Rigid3d", showPositions=0, PositionsColor="1 0 1 1", ExportPositions=False,
            showVelocities=False, VelocitiesColor="0.5 0.5 1 1", ExportVelocities=False,
            showForces=0, ForcesColor="0.8 0.2 0.2 1", ExportForces=True, showTrajectories=0,
            TrajectoriesPrecision="0.1", TrajectoriesColor="0 1 1 1", sizeFactor="1",
            fileName=DirectoryResults + strbeam + 'ReactionForce')
            
        CentralBeam.createObject('Monitor', name="ReactionForceCentralBeam", indices='0',
        template="Rigid3d", showPositions=0, PositionsColor="1 0 1 1", ExportPositions=False,
        showVelocities=False, VelocitiesColor="0.5 0.5 1 1", ExportVelocities=False,
        showForces=0, ForcesColor="0.8 0.2 0.2 1", ExportForces=True, showTrajectories=0,
        TrajectoriesPrecision="0.1", TrajectoriesColor="0 1 1 1", sizeFactor="1",
        fileName=DirectoryResults + 'CentralBeamReactionForce')

        CentralBeam.createObject('Monitor', name="DisplacementEndNode", indices=indexNodeBC,
        template="Rigid3d", showPositions=0, PositionsColor="1 0 1 1", ExportPositions=1 ,
        showVelocities=False, VelocitiesColor="0.5 0.5 1 1", ExportVelocities=False,
        showForces=0, ForcesColor="0.8 0.2 0.2 1", ExportForces=0, showTrajectories=0,
        TrajectoriesPrecision="0.1", TrajectoriesColor="0 1 1 1", sizeFactor="1",
        fileName=DirectoryResults + 'CentralBeamDisplacementEnd')

        return 0;

    def onMouseButtonLeft(self, mouseX,mouseY,isPressed):
        ## usage e.g.
        #if isPressed : 
        #    print "Control+Left mouse button pressed at position "+str(mouseX)+", "+str(mouseY)
        return 0;

    def onKeyReleased(self, c):
        ## usage e.g.
        #if c=="A" :
        #    print "You released a"
        return 0;

    def initGraph(self, node):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onKeyPressed(self, c):
        ## usage e.g.
        #if c=="A" :
        #    print "You pressed control+a"
        return 0;

    def onMouseWheel(self, mouseX,mouseY,wheelDelta):
        ## usage e.g.
        #if isPressed : 
        #    print "Control button pressed+mouse wheel turned at position "+str(mouseX)+", "+str(mouseY)+", wheel delta"+str(wheelDelta)
        return 0;

    def storeResetState(self):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def cleanup(self):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onGUIEvent(self, strControlID,valueName,strValue):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onEndAnimationStep(self, deltaTime):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        t = self.rootNode.time
        # print "t="
        # print str(t)
        # print "\n"
        if t > 5:
            print "End of the loading reached. Post process starts"
            self.rootNode.animate = False
            print '\n \n '
            import os
            import matplotlib.pylab as plt
            plt.close('all')

            fileDir = self.DirectoryResults 
            fileExt = r".txt"

            L = [_ for _ in os.listdir(fileDir) if _.endswith(fileExt)]

            arr = np.zeros(np.loadtxt(fileDir+L[0]).shape)

            for i, file in enumerate(L):
                if 'ReactionForce' in file:
                    arr += np.loadtxt(fileDir+file)
                else:
                    posDirichlet = np.loadtxt(fileDir+file)

            t = np.loadtxt(fileDir+L[0])[:,0]
            RFz = arr[:,3]

            uDirichlet = posDirichlet - posDirichlet[0]

            d = load_obj(fileDir+'infoSimu')
            # length of the strand
            h = d['lengthStrand']
            # compute the strand axial strain
            epsilon = uDirichlet[:,3] / h

            filter = t <= 10

            plt.figure()
            plt.plot(t[filter], epsilon[filter])
            ax = plt.gca()
            ax.set_xlabel('time')
            ax.set_ylabel(r'$\epsilon$ strand $\frac{l}{L}$')
            plt.savefig(fileDir+ 'timeVSepsilon')

            plt.figure()
            plt.plot(epsilon[filter], RFz[filter], label='Sofa model with BFE')
            # currve from Costello theory. Seen in the paper of Jiang, 1999  
            # I just took 2 points on the curve of Fig. 5
            slopeCostello = (150E3)/0.011
            plt.plot(epsilon[filter], slopeCostello * epsilon[filter], label='Costello', linestyle='--')
            ax = plt.gca()
            ax.set_xlabel(r'$\epsilon$ strand $\frac{l}{L}$')
            ax.set_ylabel(r'$F_{z}$')
            ax.legend()
            plt.savefig(fileDir+ 'epsilonVSAxialLoad')
            quit()
        
        return 0;

    def onLoaded(self, node):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def reset(self):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onMouseButtonMiddle(self, mouseX,mouseY,isPressed):
        ## usage e.g.
        #if isPressed : 
        #    print "Control+Middle mouse button pressed at position "+str(mouseX)+", "+str(mouseY)
        return 0;

    def bwdInitGraph(self, node):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onScriptEvent(self, senderNode, eventName,data):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        return 0;

    def onMouseButtonRight(self, mouseX,mouseY,isPressed):
        ## usage e.g.
        #if isPressed : 
        #    print "Control+Right mouse button pressed at position "+str(mouseX)+", "+str(mouseY)
        return 0;

    def onBeginAnimationStep(self, deltaTime):
        ## Please feel free to add an example for a simple usage in /Users/marco.magliulo/mySofaCodes/myScriptsToGetStarted//Users/marco.magliulo/Software/sofa/src/applications/plugins/SofaPython/scn2python.py
        # nodeBlocked = self.rootNode.getChild('CentralBeam').getObject('FixedConstraint').findData('indices').value
        # print 'node blocked is ' + str(nodeBlocked)
        # import pdb; pdb.set_trace()
        return 0;

def createScene(rootNode):
    rootNode.findData('dt').value = '0.01'
    rootNode.findData('gravity').value = '0 0 -9.81'
    try : 
        sys.argv[0]
    except :
        commandLineArguments = []
    else :
        commandLineArguments = sys.argv
    myBeamFEMForceField = BeamFEMForceField(rootNode,commandLineArguments)
    return 0;

I launch the simulation as follows:
runSofa TensileTestOnSixPlusOneStrand.py --argv 150 -g batch --start -n 10000,
where 150 is the number of nodes used per body.

No problem when we have 100 nodes, but for 150 nodes and above, the warning messages and segfault error mentioned above occur.

Anyone has experienced something similar?

Cheers