import febio import sys, string import pickle import hex_cent import numpy as np from scipy.cluster.vq import kmeans2 def main(inpfile, pklfile): print "----- ----- ----- ----- ----- Start making cell model ----- ----- ----- ----- -----" mesh = febio.MeshDef(inpfile,'abq') model = febio.Model(modelfile=string.replace(inpfile,'.inp','lumped.feb'),steps=[{'Displace':'poro'}]) # element/nodes N = len(mesh.elements) elms = np.zeros((N,8),int) for j in xrange(N): elms[j,:] = mesh.elements[j][2:] N = len(mesh.nodes) nodes = np.zeros((N,3),float) for j in xrange(N): nodes[j,:] = mesh.nodes[j][1:] #element centroids N = elms.shape[0] el = np.asfortranarray(elms) n = np.asfortranarray(nodes) cent = np.asfortranarray(np.zeros((N,3),float)) hex_cent.hexcalc(el,n,cent) d_set = {} #dictionary to store new element sets in as lists eset = mesh.elsets['epcm'] r = [] ids = [] for eid in eset: r.append(np.linalg.norm(cent[eid-1,:])) ids.append(eid) r = np.array(r) k = np.linspace(np.amin(r),np.amax(r),12,endpoint=True) output = kmeans2(r,k,iter=1000,minit='matrix') inds = output[1] binsort = np.argsort(output[0]) newinds = [] for i in inds: newinds.append(binsort[i]) inds = newinds for i in xrange(len(eset)): try: d_set[inds[i]+1].append(ids[i]) except: d_set[inds[i]+1] = [] d_set[inds[i]+1].append(ids[i]) for i in d_set.keys(): mesh.addElementSet(setname=str(i),eids=d_set[i]) pcm_sets = np.sort(d_set.keys()) # material properties # .. E [MPa], perm [10^-15 m^4 / Ns] props = { 'ecell':{'E':'0.00159','v':'0.34','perm':'1.0e-3'}, 'emem':{'E':'0.1','v':'0.04','perm':'3.00e-6'}, '1':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '2':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '3':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '4':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '5':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '6':{'E':'0.1266767198046819','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '7':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '8':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '9':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '10':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '11':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}, '12':{'E':'0.18001242649289148','v':'0.04','perm':'4.00e-2','M':'4.638','alpha':'0.0848','e0':'4.00'}} # create materials mats = [] cnt = 0 for s in map(str,list(pcm_sets)): cnt += 1 mats.append(febio.MatDef(matid=cnt,mname=s,elsets=[s],mtype='biphasic',attributes={'phi0':str(1/(1+float(props[s]['e0'])))})) mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean', attributes={'E':props[s]['E'],'v':props[s]['v']}) mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-Holmes-Mow', attributes={'perm':props[s]['perm'],'M':props[s]['M'],'alpha':props[s]['alpha']}) model.addMaterial(mat=mats[-1]) mats.append(febio.MatDef(matid=cnt+1,mname='eecm',elsets=['eecm'],mtype='biphasic',attributes={'phi0':'0.2'})) mats[-1].addBlock(branch=1,btype='solid',mtype='solid mixture',attributes={'mat_axis':['local','1,2,4']}) mats[-1].addBlock(branch=2,btype='solid',mtype='neo-Hookean',attributes={'E':'0.27','v':'0.04'}) mats[-1].addBlock(branch=2,btype='solid',mtype='ellipsoidal fiber distribution',attributes={'ksi':'4.0,4.0.,4.0', 'beta':'2.0,2.0,2.0'}) mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-ref-ortho',attributes={'perm0':'0.0', 'perm1':'7.6e-3,7.6e-3,7.6e-3','perm2':'0,0,0','M0':'0','alpha0':'0', 'M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'}) model.addMaterial(mat=mats[-1]) mats.append(febio.MatDef(matid=cnt+2,mname='ecell',elsets=['ecell'],mtype='biphasic',attributes={'phi0':'0.2'})) mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':props['ecell']['E'],'v':props['ecell']['v']}) mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':props['ecell']['perm']}) model.addMaterial(mat=mats[-1]) mats.append(febio.MatDef(matid=cnt+3,mname='emem',elsets=['emem'],mtype='biphasic',attributes={'phi0':'0.2'})) mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':props['emem']['E'],'v':props['emem']['v']}) mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':props['emem']['perm']}) model.addMaterial(mat=mats[-1]) model.addGeometry(mesh=mesh,mats=mats) # load interpolated data from pickle file pkl_file = open(pklfile, 'rb') load_curves = pickle.load(pkl_file) pkl_file.close() nnodes = len(load_curves) # number of nodes print "... Nr of nodes retrieved from pickle:\t", nnodes ctrl = febio.Control() # prescribe load curves print "... ... ... Prescribing load curves" itr = 0 # itr : to name the load curves countn = 0 boundary = febio.Boundary(steps=1) for n in load_curves.keys(): # for every dictionary in load_curves ndata = load_curves[n] # .. which represents the information about one node node_id = ndata["nid"] # .. retrieve time, displacements, and fluid pressure ## node_id = n # ADJUSTED NOV 25 for rinterp.py time = ndata["t"] countn += 1 # count nodes nsteps = len(time) # number of time steps for idof in ["x", "y", "z", "p"]: dof_data = ndata[idof] load_curve = 2*['' for i in time] # make list : [t1, dof1, t2, dof2, ..] itr += 1 for s in xrange(nsteps): i = s*2 # index number for times t1, t2, .. j = s*2+1 # index number for dofs dof1, dof2, .. load_curve[i]=time[s] # .. add info to list "load_curve" load_curve[j]=dof_data[s] model.addLoadCurve(lc=str(itr), lctype="linear", points=load_curve) # write loadcurve boundary.addPrescribed(nodeid=node_id, dof=idof, lc=str(itr), scale=1.0) # create must points loadcurve print "... ... ... Prescribing must points curve" mp_lc_id = itr+1 # must points loadcurve ID tdiff = [0 for i in xrange(nsteps-1)] # list with time differences between two consecutive time pts for t in xrange(nsteps-1): tdiff[t] = time[t+1]-time[t] tdiff = [0]+tdiff mp = 2*['' for ts in time] # list with entries for must points for s in xrange(nsteps): i = 2*s j = 2*s + 1 mp[i] = time[s] mp[j] = float(tdiff[s]) print "... ... ... Length of mp curve:", len(mp)/2 model.addLoadCurve(lc=str(mp_lc_id),lctype='step',points=mp) ctrl.setAttributes({'time_steps':'3000','step_size':'0.1','plot_level':'PLOT_DEFAULT', 'time_stepper': {'aggressiveness': '0', 'dtmin': '0.001', 'dtmax': 'lc='+str(mp_lc_id), 'max_retries': '10', 'opt_iter': '10'}, 'max_ups':'0','max_refs':'10','dtol':'0.001','etol':'0','rtol':'0','lstol':'0.9','ptol':'0.1'}) model.addControl(ctrl, step=0) model.addBoundary(boundary) model.addOutput("nodal fluid flux") model.writeModel() print "----- ----- ----- ----- ----- Finished making cell model ----- ----- ----- ----- -----" if __name__ == "__main__": main(sys.argv[-2], sys.argv[-1])