Merge branch 'master' of https://bitbucket.org/samop/trisurf-manager

Miha

2019-03-25 f039f4c3e44d77456521013030fb3e955e7d9a61

Merge branch 'master' of https://bitbucket.org/samop/trisurf-manager

 trisurf/VTKRendering.py

@@ -4,7 +4,7 @@
    from vtk import *
except:
    print("Vtk rendering works if you manually install vtk7 for python3")
    exit(1)
#    exit(1)

class MultiRender:
    def __init__(self,args,host):
@@ -182,3 +182,16 @@
        #self.render.RemoveActor(self.actor)
        
        return


class RenderVTUFile(Renderer):
    def __init__(self,args,filename=None,scalar_field='vertices_idx'):
        self.filename=filename
        Renderer.__init__(self,args,host=None, run=None, timestep=-1, scalar_field=scalar_field)

    def lastVTU(self):
        if(self.filename):
            return self.filename
        else:
            raise Exception('No file specified')


 trisurf/analyses.py

@@ -36,7 +36,7 @@
    from trisurf import VTKRendering as vtk
    import math
    from multiprocessing import Process
    table=trisurf.Statistics(run.Dir.fullpath(),filename='data_tspoststat.csv').getTable()
    table=trisurf.Statistics(run.Dir.fullpath(),filename='poststat.csv').getTable()
    def running_avg(col):
        import numpy as np
        avg=[]    
@@ -66,7 +66,7 @@
# these functions should be wrapped
# -------------------------------
@analysis('plotColumnFromPostProcess')
def plotColumnFromPostProcess(run, filename='data_tspoststat.csv', column='hbar', **kwargs):
def plotColumnFromPostProcess(run, filename='poststat.csv', column='hbar', **kwargs):
    import matplotlib.pyplot as plt

    def smooth(y, box_pts):

 trisurf/vtu-compress-in-development.py

New file
@@ -0,0 +1,93 @@
"""
Binary data in the file are compressed when the VTKFile element is of the form

 <VTKFile ... compressor="vtkZLibDataCompressor">

The data corresponding to a data array are stored in a set of blocks which are each compressed using the zlib library. The block structure allows semi-random access without decompressing all data. In uncompressed form all the blocks have the same size except possibly the last block which might be smaller. The data for one array begin with a header of the form

 [#blocks][#u-size][#p-size][#c-size-1][#c-size-2]...[#c-size-#blocks][DATA]

Each token is an integer value whose type is specified by "header_type" at the top of the file (UInt32 if no type specified). The token meanings are:

 [#blocks] = Number of blocks
 [#u-size] = Block size before compression
 [#p-size] = Size of last partial block (zero if it not needed)
 [#c-size-i] = Size in bytes of block i after compression

The [DATA] portion stores contiguously every block appended together. The offset from the beginning of the data section to the beginning of a block is computed by summing the compressed block sizes from preceding blocks according to the header. 


Usage example:
compressvtkfile('timestep_000099.vtu')
"""

import xml.etree.ElementTree as ET
import struct
import base64
import zlib

def element2binary(element, pack_type):
    t=element.text.split()
    if(pack_type=='q' or pack_type=='B'): #Int64 or UInt8
        s=[int(i) for i in t]
    elif(pack_type=='d'): #double
        s=[float(i) for i in t]	
    else:
        print("wrong type")
        exit(1)

    if(pack_type in 'qd'):
        mult=8
    else:
        mult=1

    #if(pack_type=='q'): #reduce integer to int32
    #    mult=4
    #    pack_type='L'
    #    element.set('type','UInt32')	
    b=struct.pack(pack_type*len(s),*s)
    #blen=struct.pack('I',mult*len(s))
    cmprs=zlib.compress(b,6)
    sizes=[1,len(b),len(b),len(cmprs)]
    header=struct.pack('I'*4,*sizes)
    element.text=base64.b64encode(header).decode('ascii')+base64.b64encode(cmprs).decode('ascii')
    element.set('format','binary')	
	
def compressvtkfile(filename):
    try:
        tree = ET.parse(filename)
    except:
        print("Error reading snapshot file")

    iterator=tree.iter()

    for element in iterator:
        #UInt8, Int64, Float64
        if element.get('type')=='UnstructuredGrid':
            element.set('compressor','vtkZLibDataCompressor')
        if element.get('type')=='Int64':
            element2binary(element,'q')

        if element.get('type')=='UInt8':
            element2binary(element,'B')
        if element.get('type')=='Float64':
            element2binary(element,'d')

    tree.write("test.vtu")

    #now compress the trisurf part:
    with open('test.vtu', 'r') as file:
        c=file.read()
    i=c.find('<trisurf ')
    j=c.find('>',i)
    k=c.find('</trisurf>',j)
    cs=c[j+1:k]
    fields=c[i:j+1].split()
    fields=[s.replace('false','true') for s in fields]
    starttag=' '.join(fields)
    cmprs=zlib.compress(cs.encode('ascii'),6)
    cmprs=base64.b64encode(cmprs)
    cc=c[0:i]+starttag+cmprs.decode('ascii')+c[k:]
    with open('test1.vtu', 'w') as file:
        file.write(cc)


 trisurf/wrapper.py

New file
@@ -0,0 +1,293 @@
from ctypes import *

TS_SUCCESS=0
TS_FAIL=1

TS_ID_FILAMENT=1

TS_COORD_CARTESIAN=0
TS_COORD_SPHERICAL=1
TS_COORD_CYLINDRICAL=2


class ts_coord(Structure):
    _fields_=[
        ("e1", c_double),
        ("e2", c_double),
        ("e3", c_double),
        ("coord_type", c_uint)
        ]
class ts_vertex(Structure):
    pass
class ts_bond(Structure):
    pass
class ts_triangle(Structure):
    pass
class ts_cell(Structure):
    pass
class ts_poly(Structure):
    pass
class ts_cluster(Structure):
    pass
ts_vertex._fields_=[
        ('idx',c_uint),
        ('x',c_double),
        ('y',c_double),
        ('z',c_double),
        ('neigh_no',c_uint),
        ('neigh', POINTER(POINTER(ts_vertex))),
        ('bond_length', POINTER(c_double)),
        ('bond_length_dual',POINTER(c_double)),
        ('curvature', c_double),
        ('energy', c_double),
        ('energy_h',c_double),
        ('tristar_no', c_uint),
        ('tristar', POINTER(POINTER(ts_triangle))),
        ('bond_no',c_uint),
        ('bond',POINTER(POINTER(ts_bond))),
        ('cell',POINTER(ts_cell)),
        ('xk',c_double),
        ('c',c_double),
        ('id', c_uint),
        ('projArea',c_double),
        ('relR', c_double),
        ('solAngle', c_double),
        ('grafted_poly', POINTER(ts_poly)),
        ('cluster',POINTER(ts_cluster)),
        ]
class ts_vertex_list(Structure):
    _fields_=[('n',c_uint), ('vtx',POINTER(POINTER(ts_vertex)))]

ts_bond._fields_=[('idx',c_uint),
        ('vtx1', POINTER(ts_vertex)),
        ('vtx2', POINTER(ts_vertex)),
        ('bond_length',c_double),
        ('bond_length_dual',c_double),
        ('tainted', c_char),
        ('energy',c_double),
        ('x',c_double),
        ('y',c_double),
        ('z',c_double),
    ]
class ts_bond_list(Structure):
    _fields_=[('n', c_uint),('bond',POINTER(POINTER(ts_bond)))]

ts_triangle._fields_=[
        ('idx',c_uint),
        ('vertex', POINTER(ts_vertex)*3),
        ('neigh_no',c_uint),
        ('neigh', POINTER(POINTER(ts_triangle))),
        ('xnorm', c_double),
        ('ynorm', c_double),
        ('znorm', c_double),
        ('area', c_double),
        ('volume', c_double),
        ('energy', c_double),
    ]

class ts_triangle_list(Structure):
    _fields_=[('n',c_uint),('a0',c_double),('tria', POINTER(POINTER(ts_triangle))),]


ts_cell._fields_=[
    ('idx', c_uint),
    ('vertex', POINTER(POINTER(ts_vertex))),
    ('nvertex', c_uint),
    ]		

class ts_cell_list(Structure):
    _fields_=[
        ('ncmax', c_uint*3),
        ('cellno', c_uint),
        ('cell',POINTER(POINTER(ts_cell))),
        ('dcell', c_double),
        ('shift', c_double),
        ('max_occupancy', c_double),
        ('dmin_interspecies', c_double)
    ]

class ts_spharm(Structure):
    _fields_=[
        ('l',c_uint),
        ('ulm', POINTER(POINTER(c_double))),
    #    ('ulmComplex', POINTER(POINTER(gsl_complex))), #poisci!!!!
        ('ulmComplex', POINTER(POINTER(c_double))), #temporary solution
        ('sumUlm2', POINTER(POINTER(c_double))),
        ('N', c_uint),
        ('co',POINTER(POINTER(c_double))),
        ('Ylmi', POINTER(POINTER(c_double))),
        ]

ts_poly._fields_=[
        ('vlist', POINTER(ts_vertex_list)),
        ('blist', POINTER(ts_bond_list)),
        ('grafted_vtx',POINTER(ts_vertex)),
        ('k', c_double),
    ]

class ts_poly_list(Structure):
    _fields_=[('n',c_uint),('poly',POINTER(POINTER(ts_poly)))]

class ts_tape(Structure):
    _fields_=[
        ('nshell',c_long),
        ('ncxmax',c_long),
        ('ncymax',c_long),
        ('nczmax',c_long),
        ('npoly',c_long),
        ('nmono',c_long),
        ('internal_poly',c_long),
        ('nfil',c_long),
        ('nfono', c_long),
        ('R_nucleus',c_long),
        ('R_nucleusX',c_double),
        ('R_nucleusY',c_double),
        ('R_nucleusZ',c_double),
        ('pswitch',c_long),
        ('constvolswitch',c_long),
        ('constareaswitch',c_long),
        ('stretchswitch',c_long),
        ('xkA0',c_double),
        ('constvolprecision',c_double),
        ('multiprocessing',c_char_p),
        ('brezveze0',c_long),
        ('brezveze1',c_long),
        ('brezveze2',c_long),
        ('xk0',c_double),
        ('dmax',c_double),
        ('dmin_interspecies',c_double),
        ('stepsize',c_double),
        ('kspring',c_double),
        ('xi',c_double),
        ('pressure',c_double),
        ('iterations',c_long),
        ('inititer',c_long),
        ('mcsweeps',c_long),
        ('quiet',c_long),
        ('shc',c_long),
        ('number_of_vertice_with_c0',c_long),
        ('c0', c_double),
        ('w', c_double),
        ('F', c_double),
    ]
		

class ts_vesicle(Structure):
    _fields_=[
        ('vlist', POINTER(ts_vertex_list)),
        ('blist', POINTER(ts_bond_list)),
        ('tlist', POINTER(ts_triangle_list)),
        ('clist', POINTER(ts_cell_list)),
        ('nshell', c_uint),
        ('bending_rigidity',c_double),
        ('dmax',c_double),
        ('stepsize',c_double),
        ('cm', c_double*3),
        ('volume', c_double),
        ('sphHarmonics',POINTER(ts_spharm)),
        ('poly_list', POINTER(ts_poly_list)),
        ('filament_list', POINTER(ts_poly_list)),
        ('spring_constant', c_double),
        ('pressure', c_double),
        ('pswitch', c_int),
        ('tape', POINTER(ts_tape)),
        ('R_nucleus', c_double),
        ('R_nucleusX', c_double),
        ('R_nucleusY', c_double),
        ('R_nucleusZ', c_double),
        ('nucleus_center', c_double *3 ),
        ('area', c_double),
    ]

ts_cluster._fields_=[('nvtx',c_uint),('idx',c_uint),('vtx', POINTER(POINTER(ts_vertex)))]

class ts_cluster_list(Structure):
    _fields_=[('n',c_uint),('cluster',POINTER(POINTER(ts_cluster)))]




ts=CDLL('libtrisurf.so')



#function call wrappers
def create_vesicle_from_tape(tape):
    """Using pointer for tape, it creates a vesicle, returning pointer to it."""
    ts.create_vesicle_from_tape.argtype=POINTER(ts_tape)
    ts.create_vesicle_from_tape.restype=POINTER(ts_vesicle)
    return ts.create_vesicle_from_tape(tape)

def parsetape(filename='tape'):
    """Loads tape with  filename (if not given it defaults to 'tape'). It returns a pointer to structure for tape"""
    ts.parsetape.restype=POINTER(ts_tape)
    ts.parsetape.argtype=[c_char_p]
    return ts.parsetape(filename.encode('ascii'))

def parseDump(filename):
    """Loads a vtu file with 'filename' and creates a vesicle returning pointer to it"""
    ts.parseDump.argtype=[c_char_p]
    ts.parseDump.restype=POINTER(ts_vesicle)
    vesicle=ts.parseDump(filename.encode('ascii'))
    return vesicle

def single_timestep(vesicle):
    """Makes a single timestep in simulations. Returns a tuple of vmsrt and bfrt (vertex move success rate and bond flip success rate)"""
    ts.single_timestep.argtype=[POINTER(ts_vesicle),POINTER(c_double),POINTER(c_double)]
    vmsrt=c_double(0.0)
    bfsrt=c_double(0.0)
    ts.single_timestep(vesicle,byref(vmsrt),byref(bfsrt))
    return (vmsrt.value, bfsrt.value)

def write_vertex_xml_file(vesicle,timestep_no=0):
    """Writes a vesicle into file with filename 'timestep_XXXXXX.vtu', where XXXXXX is a leading zeroed number given with timestep_no parameter (defaults to 0 if not given"""
    ts.write_vertex_xml_file.argtypes=[POINTER(ts_vesicle),c_int]
    ts.write_vertex_xml_file(vesicle,c_int(timestep_no))


def vesicle_free(vesicle):
    """Free memory of the whole vesicle"""
    ts.vesicle_free.argtype=[POINTER(ts_vesicle)]
    ts.vesicle_free(vesicle)

def vesicle_volume(vesicle):
    ts.vesicle_volume.argtype=[POINTER(ts_vesicle)]
    ts.vesicle_volume(vesicle)

def vesicle_area(vesicle):
    ts.vesicle_area.argtype=[POINTER(ts_vesicle)]
    ts.vesicle_area(vesicle)

def gyration_eigen(vesicle):
    ts.gyration_eigen.argtype=[POINTER(ts_vesicle), POINTER(c_double), POINTER(c_double), POINTER(c_double)]
    l1=c_double(0.0)
    l2=c_double(0.0)
    l3=c_double(0.0)
    ts.gyration_eigen(vesicle , byref(l1), byref(l2), byref(l3))
    return (l1.value, l2.value, l3.value)

def vesicle_meancurvature(vesicle):
    ts.vesicle_meancurvature.argtype=[POINTER(ts_vesicle)]
    ts.vesicle_meancurvature.restype=c_double
    return ts.vesicle_meancurvature(vesicle)

def init_cluster_list():
    ts.init_cluster_list.restype=POINTER(ts_cluster_list)
    ret=ts.init_cluster_list()
    return ret

def clusterize_vesicle(vesicle, cluster_list):
    ts.clusterize_vesicle.argtype=[POINTER(ts_vesicle), POINTER(ts_cluster_list)]
    ts.clusterize_vesicle(vesicle, cluster_list)

def cluster_list_free(cluster_list):
    """Free memory of cluster list"""
    ts.cluster_list_free.argtype=[POINTER(ts_cluster_list)]
    ts.cluster_list_free(cluster_list)

def stretchenergy(vesicle, triangle):
    ts.stretchenergy.argtype=[POINTER(ts_vesicle), POINTER(ts_triangle)]
    ts.stretchenergy(vesicle,triangle)