forceChain.py

# Thanks to the useful discussion with:
# 1- Admin, 2- Stefan, 3-Vanessa
import os


def FlipString(stringToFlip):
    return stringToFlip[::-1]


def FindTimeFrame(direc, x=4):
    """
    :param x: is the number of integers padded in the file e.g. x = 3 if it is file_000, x= 4 if file_0000
    :param direc: is the name of the directory were the files are exists
    :return: list of integers representing the end of the padded files.
    """
    print(f'I am looking in side {direc} directory!')
    directory = os.listdir(direc)
    number = []
    print('and I am searching for files with ".data." extension')
    extensionString = '.data.'
    for fname in directory:
        if '.data.' in fname:

            # flip the name
            flippedName = FlipString(fname)
            num = flippedName[0:x]
            fixedName = FlipString(num)
            num = int(fixedName)
            number.append(num)

    print(f'number of the files with ".data." extension is {len(number)}!')
    name = ''
    for fname in directory:
        if '.data.' in fname:
            pathname, extension = os.path.splitext(fname)
            filename = pathname.split('/')

            name = filename[-1].split('.')
            name = name[0]
            continue
    return number, name


def Periodicity(periodicIn, dataFile, t):
    tStr = str(t)
    tt = tStr.zfill(4)
    # periodic in x then get xMin and xMax
    dataFile = dataFile + ".data." + tt
    dataObj = open(dataFile, 'r')
    header = dataObj.readlines()[0:1]
    dataObj.close()

    # extract the info and save it in [x_min, x_max]
    xCoord = []
    yCoord = []
    zCoord = []
    for items in header:
        item = items.lstrip()
        domain = item.split()[2:]
        xCoord = [float(domain[0]), float(domain[3])]
        yCoord = [float(domain[1]), float(domain[4])]
        zCoord = [float(domain[2]), float(domain[5])]

    # return type

    if periodicIn == 'X':
        return xCoord
    if periodicIn == 'Y':
        return yCoord
    if periodicIn == 'Z':
        return zCoord
    if periodicIn == 'XY' or periodicIn == 'YX':
        return xCoord, yCoord
    else:
        return 0


def ConditionCheck(position, rad, id1, id2, domain, periodicIn):
    """
    :param id2:
    :param position: List, contain Position of the particle [x,y,z].
    :param id1: integer, of the particle which is in contact with anther one.
    :param rad: float, for the particle with id1 represent the radius.
    :param domain: list of floats, contain the bounding box depending on the
    :param periodicIn: char shows which direction is periodic
    :return: boolean type [true or false]
    """
    idList = [id1, id2]
    for i in idList:
        if periodicIn == 'X':
            r = rad[str(i)]
            if position[str(i)][0] > (domain[0] + r) or position[str(i)][0] < (domain[1] - r):
                return True
        elif periodicIn == 'Y':
            if position[str(i)][1] > (domain[0] + rad[str(i)]) or position[str(i)][1] < (domain[1] + rad[str(i)]):
                return True
        elif periodicIn == 'XY':
            x = domain[0]
            xx = x[0]
            y = domain[1]
            r = 2.2 * rad[str(i)]
            p = position[str(i)]
            if (x[0] + r) < p[0] < (x[1] - r):
                if (y[0] + r) < p[1] < (y[1] - r):
                    return True
        else:
            return False


def ForceChains(file, time, periodicIn='None', comment="comment"):
    # Bounding domain
    domain = Periodicity(periodicIn, file, time)

    # time to string and padded
    tStr = str(time)
    tt = tStr.zfill(4)
    # output file
    fileName = file + "ForceChain_" + str(time) + '.vtp'
    # create a vtp file and write to it
    vtpObj = open(fileName, 'w')

    # Read particle data file ---------------------------------------------------------
    dataObj = open(file + ".data." + tt, 'r')
    # The range may need to be changed
    lines = dataObj.readlines()[1:]
    dataObj.close()

    # open the data file to extract the position and the radii of the particles -------
    nParticle = 0
    radius = dict()
    position = dict()
    for line in lines:
        l1 = line.lstrip()
        l1 = l1[:-2].split(' ')
        # print l1
        # b_id = int(l1[0])
        position[str(nParticle)] = [float(l1[0]), float(l1[1]), float(l1[2])]
        radius[str(nParticle)] = float(l1[6])
        nParticle += 1

    # Read fstat file to count number of interactions -------------------------------------------
    fstatObj = open(file + ".fstat." + tt, 'r')
    lines = fstatObj.readlines()[3:]
    fstatObj.close()
    # number of interactions excluding the one on the P.B.
    interaction = 0
    tmp = 0.0
    # get number of interactions
    for line in lines:
        line = line.lstrip()
        l1 = line[1:-6].split(' ')
        # print l1[:3],'sss',l1[3:]
        [id1, id2] = [int(i) for i in l1[1:3]]
        if id1 != tmp:
            tmp = id2
            # particle - Particle contact
            if id1 > 0 and id2 > 0:
                # if the system is periodic in x-direction
                if ConditionCheck(position, radius, id1, id2, domain, periodicIn):
                    # real number of interactions without the periodicity artifact
                    interaction += 1

    # No of interactions
    nIntrs = int(interaction)

    # Writing the vtp Header
    vtpObj.write("<?xml version='1.0'?>\n<!-- time %s-->\n" % time)
    vtpObj.write("<VTKFile type='PolyData' version='0.1' byte_order='LittleEndian'>\n<PolyData>\n")
    vtpObj.write(
        "<Piece NumberOfPoints='%s' NumberOfVerts='0' NumberOfLines='%s' NumberOfStrips='0' NumberOfPolys='0'>\n" % (
            str(2 * nIntrs), str(nIntrs)))

    # Write coords of the interacting particles (also taking into account possible periodicity
    vtpObj.write("<Points>\n<DataArray type='Float32' NumberOfComponents='3' format='ascii'>\n")
    tmp = 0.0
    for line in lines:
        line = line.lstrip()
        l1 = line[1:-6].split(' ')
        # print l1[:3],'sss',l1[3:]
        [id1, id2] = [int(i) for i in l1[1:3]]
        if id1 != tmp:
            tmp = id2
            # particle - Particle contact
            if id1 > 0 and id2 > 0:
                # if the system is periodic in x-direction
                if ConditionCheck(position, radius, id1, id2, domain, periodicIn):
                    pos = position[str(id1)]
                    vtpObj.write("%g %g %g\n" % (pos[0], pos[1], pos[2]))
                    pos = position[str(id2)]
                    vtpObj.write("%g %g %g\n" % (pos[0], pos[1], pos[2]))
    # end of this section
    vtpObj.write("</DataArray>\n</Points>\n<Lines>\n<DataArray type='Int32' Name='connectivity' format='ascii'>\n")

    # Mysterious Stuff down here
    ss = ''
    for con in range(2 * interaction):
        ss += ' ' + str(con)
    vtpObj.write(ss + '\n')
    vtpObj.write("</DataArray>\n<DataArray type='Int32' Name='offsets' format='ascii'>\n")
    ss = ''
    for con in range(interaction):
        ss += ' ' + str(con * 2 + 2)
    vtpObj.write(ss)
    vtpObj.write("\n</DataArray>\n</Lines>\n")

    # Normal force components
    name = "ForceChain"
    vtpObj.write("<PointData Scalars='%s'>\n" % name)
    name = 'fn'
    vtpObj.write("<DataArray type='Float32' Name='%s' format='ascii'>\n" % name)
    tmp = 0.0
    for line in lines:
        line = line.lstrip()
        l1 = line[:-2].split(' ')
        [id1, id2] = [int(i) for i in l1[1:3]]
        if id1 != tmp:
            tmp = id2
            if ConditionCheck(position, radius, id1, id2, domain, periodicIn):
                # forces to be shown
                [fn, ft] = [float(i) for i in l1[8:10]]
                # unit vectors
                [nx, ny, nz] = [float(i) for i in l1[10:13]]
                [fn_x, fn_y, fn_z] = [fn * nx, fn * ny, fn * nz]
                fn = (fn_x ** 2. + fn_y ** 2. + fn_z ** 2.) ** 0.5
                vtpObj.write("%g %g\n" % (fn, fn))

    vtpObj.write("</DataArray>\n")

    # Tangential Force components
    name = 'ft'
    vtpObj.write("<DataArray type='Float32' Name='%s' format='ascii'>\n" % name)
    tmp = 0.0
    for line in lines:
        line = line.lstrip()
        l1 = line.split(' ')
        [id1, id2] = [int(i) for i in l1[1:3]]
        if id1 != tmp:
            tmp = id2
            if ConditionCheck(position, radius, id1, id2, domain, periodicIn):
                # forces to be shown
                [fn, ft] = [float(i) for i in l1[8:10]]
                # unit vectors
                [tx, ty, tz] = [float(i.strip('\n')) for i in l1[13:16]]
                [ft_x, ft_y, ft_z] = [ft * tx, ft * ty, ft * tz]
                ft = (ft_x ** 2. + ft_y ** 2. + ft_z ** 2.) ** 0.5
                vtpObj.write("%g %g\n" % (ft, ft))

    vtpObj.write("</DataArray>\n</PointData>")

    vtpObj.write("\n</Piece>\n</PolyData>\n</VTKFile>")
    # close vtp file
    vtpObj.close()


def main(directory, x=4):
    # list of files numbers
    [timeFrame, name] = FindTimeFrame(directory, x)

    # Loop over the files
    for t in timeFrame:
        filePath = directory + '/' + name
        ForceChains(filePath, t, periodicIn='XY')


if __name__ == '__main__':
    main('data', x=4)