import febio
import sys, string
def main(inputf):
    mesh = febio.MeshDef(inputf,'abq')
    model = febio.Model(modelfile=string.replace(inputf,'.inp','.feb'),steps=[{'Displace':'poro'}])

    s_props = {
            'efscx':{'E':'0.27','v':'0.04','ksi':'10.0,1.0,1.0','beta':'2.0,2.0,2.0'},
            'efmcx':{'E':'0.27','v':'0.04','ksi':'4.0,4.0,4.0','beta':'2.0,2.0,2.0'},
            'efdcx':{'E':'0.27','v':'0.04','ksi':'1.0,1.0,10.0','beta':'2.0,2.0,2.0'},
            'etscx':{'E':'0.27','v':'0.04','ksi':'10.0,1.0,1.0','beta':'2.0,2.0,2.0'},
            'etmcx':{'E':'0.27','v':'0.04','ksi':'4.0,4.0,4.0','beta':'2.0,2.0,2.0'},
            'etdcx':{'E':'0.27','v':'0.04','ksi':'1.0,1.0,10.0','beta':'2.0,2.0,2.0'},
            'efscy':{'E':'0.27','v':'0.04','ksi':'1.0,10.0,1.0','beta':'2.0,2.0,2.0'},
            'efmcy':{'E':'0.27','v':'0.04','ksi':'4.0,4.0,4.0','beta':'2.0,2.0,2.0'},
            'efdcy':{'E':'0.27','v':'0.04','ksi':'1.0,1.0,10.0','beta':'2.0,2.0,2.0'},
            'etscy':{'E':'0.27','v':'0.04','ksi':'1.0,10.0,1.0','beta':'2.0,2.0,2.0'},
            'etmcy':{'E':'0.27','v':'0.04','ksi':'4.0,4.0,4.0','beta':'2.0,2.0,2.0'},
            'etdcy':{'E':'0.27','v':'0.04','ksi':'1.0,1.0,10.0','beta':'2.0,2.0,2.0'},
            'emenix':{'E':'0.55','v':'0.05','ksi':'1.0,200.0,1.0','beta':'2.0,2.0,2.0'},
            'emeniy':{'E':'0.55','v':'0.05','ksi':'200.0,1.0,1.0','beta':'2.0,2.0,2.0'}}

    f_props = {
            'efscx':{'perm0':'0','perm1':'7.6e-3,3.8e-3,3.8e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'efmcx':{'perm0':'0','perm1':'7.6e-3,7.6e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'efdcx':{'perm0':'0','perm1':'3.8e-3,3.8e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etscx':{'perm0':'0','perm1':'7.6e-3,3.8e-3,3.8e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etmcx':{'perm0':'0','perm1':'7.6e-3,7.6e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etdcx':{'perm0':'0','perm1':'3.8e-3,3.8e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'efscy':{'perm0':'0','perm1':'3.8e-3,7.6e-3,3.8e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'efmcy':{'perm0':'0','perm1':'7.6e-3,7.6e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'efdcy':{'perm0':'0','perm1':'3.8e-3,3.8e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etscy':{'perm0':'0','perm1':'3.8e-3,7.6e-3,3.8e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etmcy':{'perm0':'0','perm1':'7.6e-3,7.6e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'etdcy':{'perm0':'0','perm1':'3.8e-3,3.8e-3,7.6e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'emenix':{'perm0':'0','perm1':'0.63e-3,1.26e-3,0.63e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'},
            'emeniy':{'perm0':'0','perm1':'1.26e-3,0.63e-3,0.63e-3','perm2':'0.0','M0':'0','alpha0':'0','M':'4.638,4.638,4.638','alpha':'0.0848,0.0848,0.0848'}}
    mats = []
    cnt = 1
    for s in ['efscx','efmcx','efdcx','etscx','etmcx','etdcx','emenix','efscy','efmcy','efdcy','etscy','etmcy','etdcy','emeniy']:
        #solid phase
        mats.append(febio.MatDef(matid=cnt,mname=s,elsets=[s],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':s_props[s]['E'],'v':s_props[s]['v']})
        mats[-1].addBlock(branch=2,btype='solid',mtype='ellipsoidal fiber distribution',attributes={'ksi':s_props[s]['ksi'],
            'beta':s_props[s]['beta']})
        #fluid phase
        mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-ref-ortho',attributes={'perm0':f_props[s]['perm0'],
            'perm1':f_props[s]['perm1'],'perm2':f_props[s]['perm2'],'M0':f_props[s]['M0'],'alpha0':f_props[s]['alpha0'],
            'M':f_props[s]['M'],'alpha':f_props[s]['alpha']})
        model.addMaterial(mat=mats[-1])
        cnt += 1

    mats.append(febio.MatDef(matid=cnt,mname='femurx',elsets=['efbonex'],mtype='biphasic',attributes={'phi0':'0.2'}))
    mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':str(2000.0),'v':str(0.3)})
    mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':'3.8e-3'})
    model.addMaterial(mat=mats[-1])
    
    mats.append(febio.MatDef(matid=cnt+1,mname='tibiax',elsets=['etbonex'],mtype='biphasic',attributes={'phi0':'0.2'}))
    mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':str(2000.0),'v':str(0.3)})
    mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':'3.8e-3'})
    model.addMaterial(mat=mats[-1])

    mats.append(febio.MatDef(matid=cnt+2,mname='femurx',elsets=['efboney'],mtype='biphasic',attributes={'phi0':'0.2'}))
    mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':str(2000.0),'v':str(0.3)})
    mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':'3.8e-3'})
    model.addMaterial(mat=mats[-1])
    
    mats.append(febio.MatDef(matid=cnt+3,mname='tibiax',elsets=['etboney'],mtype='biphasic',attributes={'phi0':'0.2'}))
    mats[-1].addBlock(branch=1,btype='solid',mtype='neo-Hookean',attributes={'E':str(2000.0),'v':str(0.3)})
    mats[-1].addBlock(branch=1,btype='permeability',mtype='perm-const-iso',attributes={'perm':'3.8e-3'})
    model.addMaterial(mat=mats[-1])
    model.addGeometry(mesh=mesh,mats=mats)
    ctrl = febio.Control()
    ctrl.setAttributes({'time_steps':'3000','step_size':'0.1','plot_level':'PLOT_DEFAULT','analysis':'static',
        'time_stepper': {'aggressiveness': '0', 'dtmin': '0.001', 'dtmax': 'lc=1', 'max_retries': '10', 'opt_iter': '10'},
        'max_ups':'0','max_refs':'10','dtol':'0.01','etol':'0','ptol':'0.1','lstol':'0.9'})

    model.addControl(ctrl,step=0)

    #add load curve definitions
    model.addLoadCurve(lc='1',lctype='step',points=[0,0,0.5,0.1,1,0.1,5,1,50,5,200,20,300,50])
    model.addLoadCurve(lc='2',lctype='linear',points=[0,0,0.5,1,300,1])

    #define boundary conditions
    boundary = febio.Boundary(steps=1)
    
    boundary.addFixed(nset=mesh.nsets['nfxm'],dof='x')
    boundary.addFixed(nset=mesh.nsets['nfym'],dof='y')
    boundary.addFixed(nset=mesh.nsets['ntxm'],dof='x')
    boundary.addFixed(nset=mesh.nsets['ntym'],dof='y')
    boundary.addFixed(nset=mesh.nsets['nmxm'],dof='x')
    boundary.addFixed(nset=mesh.nsets['nmym'],dof='y')
    boundary.addFixed(nset=mesh.nsets['nfzp'],dof='xy')
    boundary.addFixed(nset=mesh.nsets['ntzm'],dof='xyz')
    boundary.addFixed(nset=mesh.nsets['nmfix'],dof='xy')
    boundary.addFixed(nset=mesh.nsets['nffix'],dof='xy')
    boundary.addFixed(nset=mesh.nsets['ntfix'],dof='xy')
    boundary.addFixed(nset=mesh.nsets['nfzp'],dof='p')
    boundary.addFixed(nset=mesh.nsets['ntzm'],dof='p')
    boundary.addFixed(nset=mesh.nsets['nmout'],dof='p')
    boundary.addFixed(nset=mesh.nsets['nfout'],dof='p')
    boundary.addFixed(nset=mesh.nsets['ntout'],dof='p')

    boundary.addPrescribed(step=0,nset=mesh.nsets['nfzp'],dof='z',lc='2',scale='-0.7')
    #contact
    boundary.addContact(step=0,ctype='sliding2',master=[mesh.fsets['ffzm']],slave=[mesh.fsets['fmzp']],attributes={'two_pass':'1','auto_penalty':'1','penalty':'2','pressure_penalty':'1','symmetric_stiffness':'0'})
    boundary.addContact(step=0,ctype='sliding2',master=[mesh.fsets['ftzp']],slave=[mesh.fsets['fmzm']],attributes={'two_pass':'1','auto_penalty':'1','penalty':'2','pressure_penalty':'1','symmetric_stiffness':'0'})

    #fill boundary conditions into model
    model.addBoundary(boundary)
    model.addOutput(output='contact pressure')


    model.writeModel()

            
if __name__ == "__main__":
    main(sys.argv[-1])