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=meshscale)
    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))
        if z_depth > meshscale*float(ccontent['*Thickness'][0]):
            continue
        z_depth = round(z_depth,3)
        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])
    
        
    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 = []       
    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)
        
        material_props = {}
        for f in fits:
            v = f[0]
            p = f[1]
            material_props[v]= p[0]+p[1]*ndepth+p[2]*ndepth**2
        for c in other_constants:
            material_props[c[0]] = c[1]
            
        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 = {'alpha': str(material_props['alpha']), 'beta': str(material_props['beta']), 'ksi': str(material_props['ksi']), 'theta': '0', 'phi': '90'}
        mat_defs[i].addBlock(branch=2,btype='solid',mtype='fiber-exp-pow',attributes=dmy)
        
        dmy = {'alpha': str(material_props['alpha']), 'beta': str(material_props['beta']), 'ksi': str(material_props['ksi']), 'theta': '90', 'phi': '90'}
        mat_defs[i].addBlock(branch=2,btype='solid',mtype='fiber-exp-pow',attributes=dmy)
        
        dmy = {'alpha': str(material_props['alpha']), 'beta': str(material_props['beta']), 'ksi': str(material_props['ksi']), 'theta': '0', 'phi': '0'}
        mat_defs[i].addBlock(branch=2,btype='solid',mtype='fiber-exp-pow',attributes=dmy)
        
        dmy = {'perm0': '0', 'perm1': str(material_props['perm11'])+','+str(material_props['perm12'])+','+str(material_props['perm13']), '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='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': '1', 'dtmin': '0.01', 'dtmax': 'lc=1', 'max_retries': '10', 'opt_iter': '10'}})
    
    #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['nzm'],dof='xyz')
    boundary.addFixed(nset=mesh.nsets['nxm'],dof='x')
    boundary.addFixed(nset=mesh.nsets['nym'],dof='y')
    boundary.addFixed(nset=mesh.nsets['nzp'],dof='p')
    boundary.addFixed(nset=mesh.nsets['nperi'],dof='p')
    boundary.addFixed(nset=mesh.nsets['nind'],dof='xy')
    boundary.addPrescribed(step=0,nset=mesh.nsets['nind'],dof='z',lc='2',scale='-.46')
    #boundary.addContact(step=0,ctype='facet-to-facet sliding',master=mesh.fsets['findbot'],slave=mesh.fsets['find'],attributes={'auto_penalty':'1','penalty':'100','two_pass':'1'})
    
    #model.addConstraint(step=0,matid=str(len(mat_defs)),dof={'trans_x':'fixed','trans_y':'fixed','trans_z':['prescribed','2','-1.0'],'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])