# -*- coding: mbcs -*-
"""
cuboid AF Model
main model of cuboid lamellar annulus fibrosus, loaded in displacement (parrallel or perpendicular to lamellar plane)
"""
## classical python modules
import math
## my abaqus module
import abaqusTools

## default abaqus modules
from abaqus import *
backwardCompatibility.setValues(reportDeprecated=False)
from abaqusConstants import *
#some caeModules seem not to be imported by abaqus cae in noGUI,
#let's import all of them here!!
from caeModules import *
#-----------------------------------------------------
def getParameters(_p={}):
    param = {}
    #GEOMETRY
    param['nbLam'] = 12
    param['height'] = 12.#radial dimension
    param['width'] = 5.#disc height
    param['length'] = 25.#circumferential direction
    #MATERIAL
    param['matType'] = 'neoHooke'#'Holzapfel' or 'neoHooke'
    param['holzapfelParameters'] = (0.022, 9.09e-4, 23.92, 1045.7, 0.)#ovine anterior - neoHokean model fit on those parameters
    param['fiberDirection'] = [math.pi/6.,-math.pi/6.,math.pi/6.,-math.pi/6.,math.pi/6.,-math.pi/6.,math.pi/6.,-math.pi/6.,math.pi/6.,-math.pi/6.,math.pi/6.,-math.pi/6.]#list of fiber angles in the (theta,z) plane
    #INTERACTIONS
    param['interfaceType'] = 'Tie'                  #'Frictionless', 'Tie', 'Friction', 'Cohesive', CohesiveFriction, 'Rough'
    param['contactStiffness'] = 1.                  #irrelevant if param['interfaceType']=='Tie'
    param['frictionCoef'] = 0.1                     #irrelevant if param['interfaceType']!='Friction'
    param['cohesivePenalties'] = (5.0, 5.0, 5.0)    #irrelevant if param['interfaceType'] does not contain 'Cohesive'
    #STEP
    param['timePeriod'] = 1.
    param['initialInc'] = 1e-5# needed for contact detection, can be larger for Tie
    param['maxInc'] = 0.08
    try:param['minInc'] = min(1e-6,_p['initialInc']/1000)
    except(KeyError):param['minInc'] = min(1e-6,param['initialInc']/1000)
    #LOAD
    param['load'] = 'horizDispl'#vertical displacement = in the plane of fibres;horizDipl = perpendicular to the fibre plane
    try:param['displ'] = 0.2*_p['length']
    except(KeyError):param['displ'] = 0.2*param['length']
    #JOB
    param['modelName'] = 'defaultName'
    param['scratchDir'] = '.'
    param['numCpus'] = 8
    param['saveCaeFile'] = True
    #
    param.update(_p)
    return param
#-----------------------------------------------------
def setElementType(eleTypeString):
    from abaqusConstants import C3D8R,C3D8RH,ENHANCED 
    if eleTypeString=='C3D8RH': eleType = mesh.ElemType(C3D8RH, hourglassControl=ENHANCED)
    else:eleType = mesh.ElemType(C3D8R, hourglassControl=ENHANCED)
    return eleType
#-----------------------------------------------------
def createCuboidPart(width,length,height,name,model):
    myRect = model.ConstrainedSketch(name+'_sketch',sheetSize=250.)
    myRect.rectangle((0,0),(width,length))
    myCuboid = model.Part(name+'_part',dimensionality=THREE_D,type=DEFORMABLE_BODY)
    myCuboid.BaseSolidExtrude(sketch=myRect, depth=height)
    return myCuboid
#-----------------------------------------------------  
def getEandNuPerp(holzapfelParam):
    c10 = holzapfelParam[0]
    d = holzapfelParam[1]
    #elastic equivalent for the full material, no fibre
    if d==0.: nu = 0.499
    else: nu = (3-2*c10*d)/(6+2*d*c10)
    E = 4*(1+nu)*c10
    return E,nu
#-----------------------------------------------------  
def getEandNu(holzapfelParam):
    c10 = holzapfelParam[0]
    d = holzapfelParam[1]
    k = holzapfelParam[2]
    #elastic equivalent for the full material, initial fibre stiffness
    if d==0.: nu = 0.499
    else: nu = (6-2*(4./3.*k+2*c10)*d)/(12+4*d*c10)
    E = 8./3.*k+4*(1+nu)*c10
    return E,nu
##################################################################################    
def setAnalysis(p):
    # check parameter consistency
    assert (len(p['fiberDirection']) == p['nbLam']), "number of fibre directions as to be equal to the number of parts!!"
	# MODEL
    myModel = mdb.Model(p['modelName'])
    myAssembly = myModel.rootAssembly
    abaqusTools.deleteDefaultModel()

    instances = list()
    directions = list()
    lowerFace = list()
    upperFace = list()
    innerFace = list()
    outerFace = list()
    
    for rec in range(p['nbLam']):
        # create cuboid
        rectangleName = 'cylinder%d'%(rec)
        myPart = createCuboidPart(p['width'],p['height'],p['length']/p['nbLam'],rectangleName,myModel)
        # coordinate system
        csysCyl = myPart.DatumCsysByThreePoints(coordSysType=CYLINDRICAL,origin=(0.,0.,0.),point1=(1.,0.,0.),point2=(0.,0.,-1.))
        # create material
        myMat = myModel.Material(name='AFLam%i'%rec])

        if isinstance(p['holzapfelParameters'],list) and len(p['holzapfelParameters']) == p['nbLam']:#there is one set of parameters per lamella
            matParam = p['holzapfelParameters'][rec]
        elif len(p['holzapfelParameters']) == 5:#there is one set of parameters
            matParam = p['holzapfelParameters']
        else: raise("parameter 'holzapfelParameters' of unknown type or wrong length")

        if p['matType'] == 'Holzapfel':
            if matParam[1]==0.:# incompressible
                myMat.Hyperelastic(table=(matParam,),materialType=ANISOTROPIC,anisotropicType=HOLZAPFEL, behaviorType=INCOMPRESSIBLE,localDirections=1)
            else:
                myMat.Hyperelastic(table=(matParam,),materialType=ANISOTROPIC,anisotropicType=HOLZAPFEL, behaviorType=COMPRESSIBLE,localDirections=1)
        elif p['matType'] == 'neoHooke': 
            if p['load'] == 'vertDispl': E,nu = getEandNu(matParam)
            elif p['load'] == 'horizDispl': E,nu = getEandNuPerp(matParam)
            if matParam[1]==0.:# incompressible
                myMat.Hyperelastic(testData=OFF,table=((E/(4*(1.+nu)),0),),materialType=ISOTROPIC,type=NEO_HOOKE, behaviorType=INCOMPRESSIBLE)
            else:
                myMat.Hyperelastic(testData=OFF,table=((E/(4*(1.+nu)),6*(1-2.*nu)/E),),materialType=ISOTROPIC,type=NEO_HOOKE, behaviorType=COMPRESSIBLE)

        # assign material
        abaqusTools.assignMaterialToPart(p['myMaterialName'][rec],myPart,myModel,orientation=None)
        instances.append(abaqusTools.createInstanceAndAddtoAssembly(myPart,myAssembly,translate=(0,0,p['length']/p['nbLam']*rec)))
        if p['matType'] == 'Holzapfel': directions.append((math.cos(p['fiberDirection'][rec]),math.sin(p['fiberDirection'][rec]),0.))
        # Find the inner/outer faces using coordinates.
        innerFacePoint = (p['width']/2.,p['height']/2.,p['length']/p['nbLam']*rec)
        innerFace.append(instances[rec].faces.findAt((innerFacePoint,) ))
        outerFacePoint = (p['width']/2.,p['height']/2.,p['length']/p['nbLam']*(rec+1))
        outerFace.append(instances[rec].faces.findAt((outerFacePoint,)))
        lowerFacePoint = (p['width']/2.,0.,p['length']/p['nbLam']*(rec+0.5))
        lowerFace.append(instances[rec].faces.findAt((lowerFacePoint,) ))
        upperFacePoint = (p['width']/2.,p['height'],p['length']/p['nbLam']*(rec+0.5))
        upperFace.append(instances[rec].faces.findAt((upperFacePoint,)))
        myAssembly.seedPartInstance(regions=(instances[rec],), size=.3*p['length']/p['nbLam'])
        myAssembly.generateMesh(regions=(instances[rec],))
        myAssembly.setElementType((instances[rec].cells,), (setElementType('C3D8RH'),))
    if p['nbLam']>1:
        ## CONSTRAINTS - same for all interfaces!!
        for nb in range(1,p['nbLam']):
            masterRegion = myAssembly.Surface(name='master%d'%(nb),side1Faces=outerFace[nb-1])
            slaveRegion = myAssembly.Surface(name='slave%d'%(nb),side1Faces=innerFace[nb])
            from interactions import Interactions
            inter = Interactions(myModel)
            inter.setMasterSlave(masterRegion,slaveRegion)
            inter.setName('interface%d'%(nb))
            if p['interfaceType'] == 'Tie': inter.setInteractionToTie()
            elif p['interfaceType'] == 'Friction': inter.setFrictionBehaviour('Friction')
            elif p['interfaceType'] == 'Cohesive': inter.setCohesiveBehaviour(useDefaultBehaviour=False,penalties=p['cohesivePenalties'])
            elif p['interfaceType'] == 'CohesiveFriction':
                inter.setCohesiveBehaviour(useDefaultBehaviour=False,penalties=p['cohesivePenalties'])
                inter.setFrictionBehaviour('Friction')
            elif p['interfaceType'] == 'Rough': inter.setFrictionBehaviour('Rough')
            inter.createInteraction()
	## STEP
    myModel.StaticStep(name='Load',previous='Initial',timePeriod=p['timePeriod'],initialInc=p['initialInc'],nlgeom=ON,
    maxInc=p['maxInc'],minInc=p['minInc'],maxNumInc=10000)

	## LOAD/BC - after step as the step names are used!!!
    if p['load'] == 'vertDispl': myModel.PinnedBC(name='Fixed',createStepName='Load',region=tuple(upperFace))
    else: myModel.PinnedBC(name='Fixed',createStepName='Load',region=tuple(innerFace[0]))

    if p['load'] == 'horizDispl': myModel.DisplacementBC(name='Displ',createStepName='Load',region=tuple(outerFace[-1]),u1=0.,u2=0.,u3=p['displ'])#positive displ is tension
    elif p['load'] == 'vertDispl': myModel.DisplacementBC(name='Displ',createStepName='Load',region=tuple(lowerFace),u1=0.,u2=-1.*p['displ'],u3=0.)#positive displ is tension
    else: raise Exception("NOT IMPLEMENTED")
   
    ## SETS FOR OUTPUT ANALYSIS
    myAssembly.Set(faces=innerFace[0], name='innerFace')
    myAssembly.Set(faces=outerFace[-1], name='outerFace')
    myAssembly.Set(faces=tuple(lowerFace), name='lowerFace')
    myAssembly.Set(faces=tuple(upperFace), name='upperFace')
    ## JOB
    from jobCreation import JobDefinition
    myJobDef = JobDefinition(p['modelName'])
    myJobDef.setScratchDir(p['scratchDir'])
    if p['matType'] == 'Holzapfel':
        myJobDef.setToFibrous()        
        myJobDef.fibreDirections(directions)
    myNewJob = myJobDef.create()    
    if p['numCpus']>1: myNewJob.setValues(numCpus=p['numCpus'],numDomains=p['numCpus'],multiprocessingMode=THREADS)
    if p['saveCaeFile']:mdb.saveAs(myNewJob.name)
    #-------------------------------------------------------
    return myNewJob,mdb
##################################################################################