import febio, string, sys import numpy as np import math def main(cfg): ccontent = parse_cfg(cfg) msh = ccontent['*MeshFile'][0] matf = ccontent['*MaterialDataFile'][0] fitvars = ccontent['*ConstantsToFit'] othervars = ccontent['*OtherConstants'] meshscale = float(ccontent['*Scale'][0]) mesh = febio.MeshDef(msh,'abq',scale=0.001) depths = [] esets = {} for i in xrange(len(mesh.elements)): corners = [] for n in mesh.elements[i][2:]: corners.append(mesh.nodes[n-1][1:]) z_depth = 0. for c in corners: z_depth += c[-1] z_depth /= float(len(corners)) z_depth = round(z_depth,3) if z_depth > meshscale*float(ccontent['*Thickness'][0]): continue try: esets[z_depth].append(i+1) except: esets[z_depth] = [] esets[z_depth].append(mesh.elements[i][1]) depths.append(z_depth) for i in esets.keys(): mesh.addElementSet(setname=str(i),eids=esets[i]) mesh.addElement(etype='quad4',corners=[0.,0.,0.,4.,4.*math.tan(math.radians(3.)),0.,4.,4.*math.tan(math.radians(3.)),1.,0.,0.,1.],name='plane') mesh.addElementData(elset='eplane',attributes={'thickness': '0.001,0.001,0.001,0.001'}) maxz = max(depths) minz = min(depths) mat_data = parse_mat(matf) fits = [] for v in fitvars: v = v.split(',') fits.append([v[1],fit_depth(mat_data[v[0]],mat_data[v[1]])]) other_constants = [] fid = open('fits.log','a') for f in fits: fid.write(f[0]+','+string.join(map(str,f[1]),',')+'\n') fid.close() for v in othervars: avg = 0 for i in mat_data[v]: avg += i other_constants.append([v, avg/len(mat_data[v])]) mat_defs = [] for i in xrange(len(depths)): #create a new febio.MatDef object for this depth mat_defs.append(febio.MatDef(matid=str(i+1),mname=str(depths[i]),elsets=[str(round(depths[i],3))],mtype='biphasic',attributes={'phi0': '0.2'})) ndepth = (depths[i]-minz)/(maxz-minz) ''' if ndepth > 0.9: w = [0.9,0.1,0.1] elif ndepth > 1./3.: w = [0.5,0.5,0.5] else: w = [0.1,0.1,0.9] ''' material_props = {} # weights for orientation Q = 0.59256083 B = 0.36698825 M = 0.65384271 nu = 0.04109246 #logistic function constants dmy = 0.8/(1+Q*np.exp(-B*(ndepth-M))**(1/nu)) #mu = 0.5 #sigma = 0.5 #dmy2 = 2*np.exp(-((ndepth-0.3)*2-mu)**2/(2*sigma**2))/(sigma*np.sqrt(2*np.pi)) #normal distribution to weight midzone orientations w = np.array([dmy,ndepth*(1-dmy)+(1-ndepth)*(dmy),1.0-dmy]) w = w/np.linalg.norm(w) for f in fits: v = f[0] p = f[1] material_props[v]= p[2]+p[1]*ndepth+p[0]*ndepth**2 for c in other_constants: material_props[c[0]] = c[1] ''' if depths[i] == minz: dmy = {'density': '1','E': '14000','v': '0.3'} mat_defs[i].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes=dmy) mat_defs[i].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':str(material_props['perm11']/3.)}) ''' mat_defs[i].addBlock(branch=1,btype='solid',mtype='solid mixture',attributes={'mat_axis': ['local','1,4,5']}) dmy = {'density': '1','E': str(material_props['E']),'v': str(material_props['v'])} mat_defs[i].addBlock(branch=2,btype='solid',mtype='neo-Hookean',attributes=dmy) dmy = {'ksi': str(material_props['ksi']*w[0])+','+str(material_props['ksi']*w[1])+','+str(material_props['ksi']*w[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': '0', 'perm1': str(material_props['perm11']*w[0])+','+str(material_props['perm11']*w[1])+','+str(material_props['perm11']*w[2]), 'M0': '0', 'alpha0': '2', 'M': '0,0,0', 'alpha': '2,2,2'} mat_defs[i].addBlock(branch=1,btype='permeability',mtype='perm-ref-ortho',attributes=dmy) #create a rigid body material for plane element mat_defs.append(febio.MatDef(matid=str(i+2),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+3),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'}) #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['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':'10000','auto_penalty':'1','two_pass':'0','search_tol':'0.01','symmetric_stiffness':'0','search_radius':'1','seg_up':'0','tension':'1','minaug':'0','maxaug':'10'}) boundary.addContact(step=0,ctype='sliding2',master=mesh.fsets['findbot'],slave=[mesh.fsets['findent'],mesh.fsets['fzp']],attributes={'penalty':'10','auto_penalty':'1','two_pass':'0','pressure_penalty':'1'}) 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() #write the model model.writeModel() def fit_depth(x,y): x = np.array(x) y = np.array(y) z = np.polyfit(x,y,2) return z def parse_mat(matf): fid = open(matf,'r') lines = map(string.rstrip,fid.readlines()) fid.close() content = {} headings = lines[0].split(',') for h in headings: content[h] = [] lines = lines[1:] for l in lines: l = l.split(',') cnt = 0 for h in headings: content[h].append(float(l[cnt])) cnt += 1 return content 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])