import febio, string, sys import math def main(cfg): ccontent = parse_cfg(cfg) msh = ccontent['*MeshFile'][0] meshscale = float(ccontent['*Scale'][0]) material_props = {'E':[0.669374921544209]*6,'v':[0.197]*6,'ksi':10.0,'perm11':2.019e-3,'beta':2.0,'M':[0.0]*6,'alpha':'2,2,2'} mesh = febio.MeshDef(msh,'abq',scale=meshscale) mesh.addElement(etype='quad4',corners=[0.,0.,-0.6,3.,3.*math.tan(math.radians(1.)),-0.6,3.,3.*math.tan(math.radians(1.)),1.,0.,0.,1.],name='plane') mesh.addElementData(elset='eplane',attributes={'thickness': '1e-6,1e-6,1e-6,1e-6'}) mat_defs = [] esets = ['ecs','ecsp','ecm','ecmp','ecd','ecdp'] ori_w = [[1.0,0.1,0.1],[.9,0.1,0.1], [1.0]*3,[1.0]*3, [0.1,0.1,1.0],[0.1,0.1,1.0]] kw = [[1.25,0.8,0.8],[1.25,0.8,0.8], [1.0,1.0,1.0],[1.0,1.0,1.0], [0.75,0.75,1.1],[0.75,0.75,1.1]] for i in xrange(len(esets)): mat_defs.append(febio.MatDef(matid=str(i+1),mname=esets[i],elsets=esets[i],mtype='biphasic', attributes={'phi0': '0.3'})) if i%2 == 0: mat_defs[i].addBlock(branch=1,btype='solid',mtype='solid mixture',attributes={'mat_axis':['local','1,2,4']}) else: mat_defs[i].addBlock(branch=1,btype='solid',mtype='solid mixture',attributes={'mat_axis':['local','3,2,1']}) dmy = {'density': '1','E': str(material_props['E'][i]),'v': str(material_props['v'][i])} mat_defs[i].addBlock(branch=2,btype='solid',mtype='neo-Hookean',attributes=dmy) dmy = {'ksi': str(material_props['ksi']*ori_w[i][0])+','+str(material_props['ksi']*ori_w[i][1])+','+str(material_props['ksi']*ori_w[i][2]), 'beta': str(material_props['beta'])+','+str(material_props['beta'])+','+str(material_props['beta'])} mat_defs[i].addBlock(branch=2,btype='solid',mtype='ellipsoidal fiber distribution',attributes=dmy) dmy = {'perm0': str(material_props['perm11']/10.0), 'perm1': str(material_props['perm11']*kw[i][0])+','+str(material_props['perm11']*kw[i][1])+','+str(material_props['perm11']*kw[i][2]),'perm2':'0,0,0', 'M0': '1', 'alpha0': '0', 'M': string.join([str(material_props['M'][i])]*3,','), 'alpha': material_props['alpha']} mat_defs[i].addBlock(branch=1,btype='permeability',mtype='perm-ref-ortho',attributes=dmy) #create bone material mat_defs.append(febio.MatDef(matid=str(i+2),mname='bone',elsets=['bone'],mtype='biphasic')) mat_defs[i+1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':str(2000.0),'v':str(0.3)}) mat_defs[i+1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':str(350)}) #create a rigid body material for plane element mat_defs.append(febio.MatDef(matid=str(i+3),mname='plane',elsets=['eplane'],mtype='rigid body',attributes={'density': '1.0', 'center_of_mass': '0,0,0'})) #create a rigid body material for the indenter mat_defs.append(febio.MatDef(matid=str(i+4),mname='indenter',elsets=['eindent'],mtype='rigid body',attributes={'density':'1.0','center_of_mass': '0,0,0'})) #create a model object modelfilename = string.replace(msh,'.inp','.feb') model = febio.Model(modelfile=modelfilename,steps=[{'Displace': 'poro'}]) #fill in material definitions of model for m in mat_defs: model.addMaterial(m) #fill in the geometry block model.addGeometry(mesh=mesh,mats=mat_defs) #create a control block object ctrl = febio.Control() ctrl.setAttributes({'time_steps': '2743','step_size': '1','plot_level':'PLOT_MUST_POINTS','time_stepper': {'aggressiveness': '0', 'dtmin': '0.01', 'dtmax': 'lc=1', 'max_retries': '10', 'opt_iter': '10'},'max_ups':'0','max_refs':'25','dtol':'0.01','etol':'0','lstol':'0.9'}) #add control object to the first and only step model.addControl(ctrl,step=0) #add load curve definitions model.addLoadCurve(lc='1',lctype='linear',points=[0,0,14.3,14.3,15.3,1,16.3,1,17.3,1,18.3,1,19.3,1,20.3,1,25,4.7,30,5,35,5,40,5,45,5,50,5,75,25,100,25,125,25,150,25,175,25,200,25,225,25,250,25,350,100,450,100,550,100,650,100,750,100,850,100,914.3,64.3,928.6,14.3,929.6,1,930.6,1,931.6,1,932.6,1,933.6,1,934.6,1,939.3,4.7,944.3,5,949.3,5,954.3,5,959.3,5,964.3,5,989.3,25,1014.3,25,1039.3,25,1064.3,25,1089.3,25,1114.3,25,1139.3,25,1164.3,25,1264.3,100,1364.3,100,1464.3,100,1564.3,100,1664.3,100,1764.3,100,1828.6,64.3,1842.9,14.3,1843.9,1,1844.9,1,1845.9,1,1846.9,1,1847.9,1,1848.9,1,1853.6,4.7,1858.6,5,1863.6,5,1868.6,5,1873.6,5,1878.6,5,1903.6,25,1928.6,25,1953.6,25,1978.6,25,2003.6,25,2028.6,25,2053.6,25,2078.6,25,2178.6,100,2278.6,100,2378.6,100,2478.6,100,2578.6,100,2678.6,100,2742.9,64.3]) model.addLoadCurve(lc='2',lctype='linear',points=[0.,0.,14.3,.1,914.3,.1,928.6,.2,1828.6,.2,1842.9,.3,2742.9,.3]) #define boundary conditions boundary = febio.Boundary(steps=1) boundary.addFixed(nset=mesh.nsets['nxm'],dof='xy') boundary.addFixed(nset=mesh.nsets['nxp'],dof='xy') boundary.addFixed(nset=mesh.nsets['nym'],dof='y') boundary.addFixed(nset=mesh.nsets['nzm'],dof='z') boundary.addFixed(nset=mesh.nsets['nzm'],dof='p') boundary.addContact(step=0,ctype='sliding-tension-compression',master=mesh.fsets['fplane'],slave=mesh.fsets['fyp'],attributes={'laugon':'0','tolerance':'0.2','gaptol':'0','penalty':'100','auto_penalty':'1','two_pass':'0','search_tol':'0.01','symmetric_stiffness':'0','search_radius':'0.1','seg_up':'0','tension':'1','minaug':'0','maxaug':'10'}) boundary.addContact(step=0,ctype='sliding2',master=mesh.fsets['findbot'],slave=[mesh.fsets['fzp']],attributes={'penalty':'5','auto_penalty':'1','two_pass':'0','pressure_penalty':'1','symmetric_stiffness':'0'}) model.addConstraint(step=0,matid=str(len(mat_defs)-1),dof={'trans_x':'fixed','trans_y':'fixed','trans_z':'fixed','rot_x':'fixed','rot_y':'fixed','rot_z':'fixed'}) model.addConstraint(step=0,matid=str(len(mat_defs)),dof={'trans_x':'fixed','trans_y':'fixed','trans_z':['prescribed','2','-.46'],'rot_x':'fixed','rot_y':'fixed','rot_z':'fixed'}) #fill boundary conditions into model model.addBoundary(boundary) #add output request to model model.addOutput(['nodal fluid flux']) #write the model model.writeModel() def parse_cfg(cfg): fid = open(cfg,'r') lines = map(string.rstrip,fid.readlines()) fid.close() content = {} for line in lines: if '**' in line: continue elif '*' in line: keywrd = line content[keywrd] = [] else: content[keywrd].append(line) return content if __name__ == '__main__': main(sys.argv[-1])