# --------------------------------------------------------------
# Python script for the stability analysis of cylindrical shells with 
# Reduced Stiffness Method
#
#
# Dr.-Ing. Ronald Wagner
# Siemens Mobility
# 18.04.2019
# --------------------------------------------------------------

#cmd_line ='abaqus cae script=Cone_LRSM_clamped_FC_005.py'
#os.system(cmd_line)
#print(cmd_line)

#---------------------------------------------------------------

# Define the navigation path

#---------------------------------------------------------------

#import os
#myScriptPath = 'C:/Users/Heinz/Desktop/Paper_LRSM_Cone/Abaqus Temp'
#myErgPath = '/res/'
#myTempPath = '/tmp/'

#os.chdir(myScriptPath+myTempPath)

#---------------------------------------------------------------

# Import the required libraries

#---------------------------------------------------------------


from abaqus import *
from abaqusConstants import *
import regionToolset
import __main__
import section
import regionToolset
import part
import material
import assembly
import step
import interaction
import load
import mesh
import job
import sketch
import visualization
import xyPlot
import connectorBehavior
import odbAccess
from operator import add

import numpy as np
#------------------------------------------------------------------------------ 
#
# section for functions - START - DO Not TOUCH
#
#------------------------------------------------------------------------------ 
#------------------------------------------------------------------------------ 

def createGeometrySet_Face(x,y,z,setname):
    """
    
    This function creates a Geometry Set (Typ - face)
    
    Input:
    Coordinate of "face" (x,y,z - Float)
    Name of "face" (setname - String)
    Ouput:
    -    

    """
    face = ()    
    p = myModel.parts['Part-1']
    f = p.faces
    myface = f.findAt((x,y,z),)
    face = face + (f[myface.index:myface.index+1], )
    p.Set(faces=face, name=setname)   

#------------------------------------------------------------------------------ 
    
def create_GeneralStiffness(Stiff_Name,A11, A12, A22, A13, A23, A33, B11, B21, B31, D11, B12, B22, B32, D12, D22, B13, B23, B33, D13, D23, D33):
    """
    
    This function creates a General Shell Stiffness Section (ABD - Matrix)
    
    Input:
    Nummer der Section (Stiff_Name - String)
    Components of ABD Matrix (A11...D33 - Float)
    
    Ouput:
    -   

    """
    p = myModel.parts['Part-1']
    myModel.GeneralStiffnessSection(name=Stiff_Name, referenceTemperature=None, stiffnessMatrix=(A11, A12, A22, A13, A23, A33, B11, B21, B31, D11, B12, B22, B32, D12, D22, B13, B23, B33, D13, D23, D33), applyThermalStress=0, poissonDefinition=DEFAULT, useDensity=OFF)
      
#------------------------------------------------------------------------------ 

def create_LRSM_Surface(LRSM_Radius):
    """
    
    This function creates the LRSM surface
    
    Input:
    Width of LRSM surface (LRSM_Radius - Float)
    
    Ouput:
    -   

    """
    y1 = 0
    p = myModel.parts['Part-1']
    v = p.vertices
    myPoint1 = p.DatumPointByCoordinate(coords=(0.0, (myRadius+my_Small_Radius)/2.0,  myHeight/2.0))
    myPoint2= p.DatumPointByCoordinate(coords=(10.0, (myRadius+my_Small_Radius)/2.0,  myHeight/2.0))
    myPoint3= p.DatumPointByCoordinate(coords=(0.0, myRadius,  0.0 ))
    myID_P1 = myPoint1.id
    myID_P2 = myPoint2.id
    myID_P3 = myPoint3.id
    d = p.datums
    myDatum_P3 = p.DatumPlaneByThreePoints(point1=d[myID_P1], point2=d[myID_P2], point3=d[myID_P3])
    myID = myDatum_P3.id
    f = p.faces[:]
    pickedRegions = f
    p.deleteMesh(regions=pickedRegions)
    f, e, d = p.faces, p.edges, p.datums
    ursprung_an = mySemi_Vertex_Angle*(-1)
    L = myHeight/(np.cos(np.pi*ursprung_an/180))
    L_half = (myHeight/2.0)/(np.cos(np.pi*ursprung_an/180))
    R_half = L_half*np.sin(np.pi*ursprung_an/180)
    t = p.MakeSketchTransform(sketchPlane=d[myID], sketchUpEdge=e[1], sketchPlaneSide=SIDE1, origin=(0.0, myRadius+R_half, myHeight/2.0))
    s = myModel.ConstrainedSketch(name='__profile__', sheetSize=882.92, gridSpacing=22.07, transform=t)
    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=SUPERIMPOSE)
    p.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES)
    #s.EllipseByCenterPerimeter(center=(0.0, 0.0), axisPoint1=(0.0, LRSM_Radius*1.0), axisPoint2=(-LRSM_Radius*3.0, 0))
    s.CircleByCenterPerimeter(center=(0.0, y1), point1=(LRSM_Radius, y1))
#    s.CircleByCenterPerimeter(center=(0.0, y1), point1=(LRSM_Radius-8, y1))
    #s.rectangle(point1=(-1.0,-LRSM_Radius), point2=(1.0, LRSM_Radius ))
    #s.rectangle(point1=(-LRSM_Radius, -2.0), point2=(LRSM_Radius, 2.0 ))
    #s.unsetPrimaryObject()
    #s.CircleByCenterPerimeter(center=(-LRSM_Radius,0.0), point1=(-LRSM_Radius-LRSM_Radius*0.25,0.0))
    #s.CircleByCenterPerimeter(center=(0.0,LRSM_Radius), point1=(0.0,LRSM_Radius+LRSM_Radius*0.25))
    f = p.faces[:]
    pickedFaces = f
    f1, e1, d2 = p.faces, p.edges, p.datums
    p.PartitionFaceBySketchDistance(sketchPlane=d2[myID], sketchUpEdge=e1[1], faces=pickedFaces, sketchPlaneSide=SIDE1, sketch=s, distance=my_Small_Radius)
    s.unsetPrimaryObject()
    del myModel.sketches['__profile__']

#------------------------------------------------------------------------------ 
    
def createConeGeometry(radiusR,height,angle):
    """
    
    This function creates the conical shells or cylinder geometry
    
    Input:
    Radius of the cone, cylinder (radiusR - Float)
    Height, Length of the cone, cylinder (height - Float)
    Semi-vertex angle of the cone, cylinder (angle - Float)
    Ouput:
    -   

    """
    s = myModel.ConstrainedSketch(name='__profile__', sheetSize=500.0)
    g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
    s.setPrimaryObject(option=STANDALONE)
    s.ConstructionLine(point1=(0.0, -250.0), point2=(0.0, 250.0))
    s.FixedConstraint(entity=g[2])
    
    s.Line(point1=(-radiusR, 0.0), point2=(-radiusR+height*tan(angle*np.pi/180), height))
    p = myModel.Part(name='Part-1', dimensionality=THREE_D, type=DEFORMABLE_BODY)
    p = myModel.parts['Part-1']
    p.BaseShellRevolve(sketch=s, angle=360.0, flipRevolveDirection=OFF)
    s.unsetPrimaryObject()
    a = myModel.rootAssembly
    a.DatumCsysByDefault(CARTESIAN)
    p = myModel.parts['Part-1']
    a.Instance(name='Part-1-1', part=p, dependent=ON)
    a = myModel.rootAssembly
    a.rotate(instanceList=('Part-1-1', ), axisPoint=(0.0, 0.0, 0.0), axisDirection=(1.0, 0.0, 0.0), angle=90.0)
    p = myModel.parts['Part-1']
    a.Instance(name='Part-1-2', part=p, dependent=ON)
    a.rotate(instanceList=('Part-1-2', ), axisPoint=(0.0, 0.0, 0.0), axisDirection=(1.0, 0.0, 0.0), angle=90.0)
    a.InstanceFromBooleanMerge(name='Part-2', instances=(a.instances['Part-1-1'], 
    a.instances['Part-1-2'], ), originalInstances=SUPPRESS, domain=GEOMETRY)
    del myModel.parts['Part-1']
    myModel.parts.changeKey(fromName='Part-2', toName='Part-1')
    p = myModel.parts['Part-1']
    myDatum = p.DatumCsysByThreePoints(name='Zylinder_KOS', coordSysType=CYLINDRICAL, origin=(0.0, 0.0, 0.0), line1=(1.0, 0.0, 0.0), line2=(0.0, 1.0, 0.0))
        
    myID = myDatum.id
    return myID
    
#------------------------------------------------------------------------------ 

def createGeometrySet_Edge(x,y,z,setname):
    """
    
    This function creates a Geometry Set (Typ - Edge)
    
    Input:
    Coordinates of the "Edge" (x,y,z - Float)
    Name "Edge" (setname - String)
    Ouput:
    -    

    """
    edge = ()    
    p = myModel.parts['Part-1']
    e = p.edges
    myedge = e.findAt((x,y,z),)
    edge = edge + (e[myedge.index:myedge.index+1], )
    p.Set(edges=edge, name=setname)   

#------------------------------------------------------------------------------
    
def createMaterial(e11,e22,nu12,g12,g13,g23,materialname):
    """
    
    This function creates a Material modell
    
    Input:
    Material parameters (e11,e22,nu12,g12,g13,g23 - Float)
    Name "Material" (materialname - String)
    Ouput:
    -    

    """
    myMaterial = myModel.Material(name=materialname)
    myMaterial.Density(table=((myLaminateDensity, ), ))
    myMaterial.Elastic(type=LAMINA, table=((e11,e22,nu12,g12,g13,g23), ))
    return

#------------------------------------------------------------------------------

def createReferencePoint(x,y,z):
    """
    
    This function creates a Reference point (Typ - reference point)
    
    Input:
    Coordinates of the Reference point (x,y,z - Float)
    Ouput:
    Name und Position of the Reference point    

    """
    refpointname = a.ReferencePoint(point=(x, y, z))
    refpointposition = r.findAt((x, y, z), )
    return refpointname,refpointposition

#------------------------------------------------------------------------------

def createMesh(partname,ele_size):
    """
    
    This function creates a FE mesh with S4R shell elements
    
    Input:
    Element size (ele_size - Float)
    Ouput:
    -   

    """
    p = myModel.parts[partname]
    p.seedPart(size=ele_size, deviationFactor=0.1, minSizeFactor=0.1)
    p = myModel.parts[partname]
    elemType1 = mesh.ElemType(elemCode=S4R, elemLibrary=STANDARD)
    elemType2 = mesh.ElemType(elemCode=S3R, elemLibrary=STANDARD)
    p = myModel.parts[partname]
    f = p.faces[:]
    faces = f
    pickedRegions =(faces, )
    p.setMeshControls(regions=f, elemShape=QUAD, algorithm=MEDIAL_AXIS)
    p.setElementType(regions=pickedRegions, elemTypes=(elemType1, elemType2))
    p.generateMesh()
    a = myModel.rootAssembly
    a.regenerate()
    p = myModel.parts[partname]
    regions = p.sets['Set-Cone']
    p.flipNormal(regions=regions)
    
#------------------------------------------------------------------------------

def createRigidBodyInteraction(RefPoint,RefEdge,RefName):
    """
    
    This function creates a "Rigid Body Interaction" between
    a reference point and an edge
    
    Input:
    Reference point and Edge for the Interaction (RefPoint,RefEdge - Float)
    Name of the Interaction (RefName - String)
    Ouput:
    -   

    """
    myRefPoint = (RefPoint, )
    myRegion = regionToolset.Region(referencePoints=myRefPoint)
    myIndex = RefEdge.index
    myRefRegion = regionToolset.Region(edges=e[myIndex:myIndex+1])
    myRigidbody = myModel.RigidBody(name=RefName, refPointRegion=myRegion, tieRegion=myRefRegion)
    
#------------------------------------------------------------------------------

def createStaticStep(Stepname,Pre_Stepname,numInc,Initial_inc,Max_inc,Min_inc):
    """
    
    This function creates a "Static Simulation" Step
    
    Input:
    Number of Increments, Start-Increment, max. Increment and min. Increment (numInc,Initial_inc,Max_inc,Min_inc - Float)
    Name of "actuell" und des "prior" Steps (Stepname,Pre_Stepname - String)
    Ouput:
    -   

    """
    myModel.StaticStep(name=Stepname, previous=Pre_Stepname, maxNumInc=numInc, initialInc=Initial_inc, maxInc=Max_inc,minInc=Min_inc, extrapolation=NONE, nlgeom=ON)
    
    # comment / uncomment for nonlinear analysis with artifical damping
    
    #myModel.steps[Stepname].setValues(maxNumInc=numInc, stabilizationMagnitude=1e-09, stabilizationMethod=DAMPING_FACTOR, continueDampingFactors=False, adaptiveDampingRatio=None, minInc=1e-15)

#------------------------------------------------------------------------------

def createDispBC_Edge_RP(RP,RPname,BCname,Step_name,trans1,trans2,trans3,rot1,rot2,rot3):
    """
    
    This function creates a Displacement boundary condition at a Reference point
    
    Input:
    Name of the point (RP - interne Variable)
    Name of the reference point (RPname - string)
    Name of the boundary condition (BCname - string)
    Name of the Steps (Step_name - string)
    Displacement and Rotation degree of freedom (trans1,trans2,trans3,rot1,rot2,rot3 - Float - oder UNSET wenn Freiheitsgrad uneingeschränkt ist)
    
    Ouput:
    -   

    """
    a = myModel.rootAssembly
    p = a.instances['Part-1-1']
    r = a.referencePoints
    refPoints1=(RP, )
    a.Set(referencePoints=refPoints1, name=RPname)
    myRefPoint = (RP, )
    myRegion = regionToolset.Region(referencePoints=myRefPoint)
    myBC = myModel.DisplacementBC(name=BCname, createStepName=Step_name, region=myRegion, u1=trans1, u2=trans2, u3=trans3, ur1=rot1, ur2=rot2, ur3=rot3)

#------------------------------------------------------------------------------

def createPlane_by_Prinicpal(Type,Offset):
    """
    
    This function creates a Principal-Plane
    
    Input:
    Typ of the Prinicpal-Plane (XY,XZ,YZ - interne Variable)
    Offset of the Plane with respect to the global coordinate system (Offset - Float)
    
    Ouput:
    Number/ID of the Principal Plane

    """
    p = myModel.parts['Part-1']
    myPlane = p.DatumPlaneByPrincipalPlane(principalPlane=Type, offset=Offset)
    myID = myPlane.id
    return myID

#------------------------------------------------------------------------------

def createVertice_Set(x,y,z,VerticeSet_name):
    """
    
    This function creates a Vertice Set (Typ - Vertice)
    
    Input:
    Coordinates of "Vertice" (x,y,z - Float)
    Name "Vertice" (setname - String)
    Ouput:
    -    

    """
    a = myModel.rootAssembly
    v = a.instances['Part-1-1'].vertices
    myVertice = v.findAt((x,y,z),)
    myVertice = (v[myVertice.index:myVertice.index+1],)
    a.Set(vertices=myVertice, name=VerticeSet_name)

#------------------------------------------------------------------------------

def createDispBC_Edge_Vertice(VerticeSet,BCname,Step_name,trans1,trans2,trans3,rot1,rot2,rot3):
    """
    
    This function creates a Displacement boundary condition at a Vertice
    
    Input:
    Name of the Vertice (VerticeSet - interne Variable)
    Name of the Boundary condition (BCname - string)
    Name of the Steps (Step_name - string)
    Displacement and Rotation degree of freedom (trans1,trans2,trans3,rot1,rot2,rot3 - Float - oder UNSET wenn Freiheitsgrad uneingeschränkt ist)
    
    Ouput:
    -   

    """
    a = myModel.rootAssembly
    myModel.DisplacementBC(name=BCname, createStepName=Step_name, region=region, u1=UNSET, u2=-myPerturbation, u3=UNSET, ur1=UNSET, ur2=UNSET, ur3=UNSET, amplitude=UNSET, fixed=OFF, distributionType=UNIFORM, fieldName='', localCsys=None)

#------------------------------------------------------------------------------

def createEdge_Set(x,y,z,EdgeSet_name):
    """
    
    This function creates a Edge Set (Typ - Edge)
    
    Input:
    Coordinates of the "Edge" (x,y,z - Float)
    Name of the "Edge" (EdgeSet_name - String)
    Ouput:
    -    

    """
    e = p.edges
    myEdge = ()
    Edge = e.findAt((x,y,z),)
    myEdge = myEdge + (e[Edge.index:Edge.index+1], )
    p.Set(edges=myEdge, name=EdgeSet_name)    


def createFace_Partition_byPlane(ID_Plane):
    """
    
    This function creates a Geometry Partition with a plane   
    Input:
    Number/ID of the plane (ID_Plane - integer)
    
    Ouput:
    -   

    """
    p = myModel.parts['Part-1']
    f = p.faces[:]
    pickedFaces = f
    d = p.datums
    p.PartitionFaceByDatumPlane(datumPlane=d[ID_Plane], faces=pickedFaces)    


def CLT(myE11,myE22,myG12,myNu12,myLaminate,myLaminateThickness,myPlyNumber):
    #--------------------------------------------------------------------------
    #
    # This function calculates the ABD Composite Stiffness Matrix Components
    #
    #--------------------------------------------------------------------------

    i=0    
    myThickness = myPlyNumber*myLaminateThickness
    Z = []
    Z.append(-myThickness/2.0)
    
    for j in range(1,myPlyNumber+1,1):
        Z.append(Z[j-1]+myLaminateThickness)
    
    # Calcuate Reduced Stiffnesses
    
    myQ11 = []
    myQ12 = []
    myQ16 = []
    myQ22 = []
    myQ26 = []
    myQ66 = []
    
    myQ11_S = []
    myQ12_S = []
    myQ16_S = []
    myQ22_S = []
    myQ26_S = []
    myQ66_S = []
    
    for i in range(0,1,1):
        myQ11.append((myE11**2)/(myE11-myNu12**2*myE22))
        myQ12.append((myNu12*myE11*myE22)/(myE11-myNu12**2*myE22))
        myQ16.append(0)
        myQ22.append((myE11*myE22)/(myE11-myNu12**2*myE22))
        myQ26.append(0)
        myQ66.append(myG12)
    
    i=0    
    for j in range(0,myPlyNumber,1):
        myQ11_S.append(myQ11[i]*np.cos(myLaminate[j]*np.pi/180.0)**4 + 2*(myQ12[i]+2*myQ66[i])*(np.cos(myLaminate[j]*np.pi/180.0))**2*(np.sin(myLaminate[j]*np.pi/180.0))**2 + myQ22[i]*(np.sin(myLaminate[j]*np.pi/180.0))**4)
        myQ12_S.append(myQ12[i]*((np.cos(myLaminate[j]*np.pi/180.0)**4)+(np.sin(myLaminate[j]*np.pi/180.0)**4))+ (myQ11[i]+myQ22[i]-4*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)**2*np.sin(myLaminate[j]*np.pi/180.0)**2)
        myQ16_S.append((myQ11[i]-myQ12[i]-2*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)**3*np.sin(myLaminate[j]*np.pi/180.0) - (myQ22[i] - myQ12[i] - 2*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)*np.sin(myLaminate[j]*np.pi/180.0)**3)
        myQ22_S.append(myQ11[i]*np.sin(myLaminate[j]*np.pi/180.0)**4 + 2*(myQ12[i]+2*myQ66[i])*(np.cos(myLaminate[j]*np.pi/180.0))**2*(np.sin(myLaminate[j]*np.pi/180.0))**2 + myQ22[i]*(np.cos(myLaminate[j]*np.pi/180.0))**4)
        myQ26_S.append((myQ11[i]-myQ12[i]-2*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)*np.sin(myLaminate[j]*np.pi/180.0)**3 - (myQ22[i] - myQ12[i] - 2*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)**3*np.sin(myLaminate[j]*np.pi/180.0))
        myQ66_S.append((myQ11[i] + myQ22[i] - 2*myQ12[i] - 2*myQ66[i])*np.cos(myLaminate[j]*np.pi/180.0)**2*np.sin(myLaminate[j]*np.pi/180.0)**2 + myQ66[i]*(np.cos(myLaminate[j]*np.pi/180.0)**4+np.sin(myLaminate[j]*np.pi/180.0)**4))
    
    
    
    # Calcualte A Matrix
    #-------------
            
    A11_v = []
    A12_v = []
    A16_v = []
    A22_v = []
    A26_v = []
    A66_v = []
    i = 1
    for j in range(0,myPlyNumber,1):
        A11_v.append(myQ11_S[j]*(Z[i]-Z[i-1]))
        A12_v.append(myQ12_S[j]*(Z[i]-Z[i-1]))
        A16_v.append(myQ16_S[j]*(Z[i]-Z[i-1]))
        A22_v.append(myQ22_S[j]*(Z[i]-Z[i-1]))
        A26_v.append(myQ26_S[j]*(Z[i]-Z[i-1]))
        A66_v.append(myQ66_S[j]*(Z[i]-Z[i-1]))
        i = i+1
    
    A11 = sum(A11_v)
    A12 = sum(A12_v)
    A16 = sum(A16_v)
    A22 = sum(A22_v)
    A26 = sum(A26_v)
    A66 = sum(A66_v)
    
    # Calcualte B Matrix
    #-------------
            
    B11_v = []
    B12_v = []
    B16_v = []
    B22_v = []
    B26_v = []
    B66_v = []
    i = 1
    for j in range(0,myPlyNumber,1):
        B11_v.append(0.5*myQ11_S[j]*(Z[i]**2-Z[i-1]**2))
        B12_v.append(0.5*myQ12_S[j]*(Z[i]**2-Z[i-1]**2))
        B16_v.append(0.5*myQ16_S[j]*(Z[i]**2-Z[i-1]**2))
        B22_v.append(0.5*myQ22_S[j]*(Z[i]**2-Z[i-1]**2))
        B26_v.append(0.5*myQ26_S[j]*(Z[i]**2-Z[i-1]**2))
        B66_v.append(0.5*myQ66_S[j]*(Z[i]**2-Z[i-1]**2))
        i = i+1
    
    
    
    
    B11 = sum(B11_v)
    B12 = sum(B12_v)
    B16 = sum(B16_v)
    B22 = sum(B22_v)
    B26 = sum(B26_v)
    B66 = sum(B66_v)
    
    # Calcualte D Matrix
    #-------------
            
    D11_v = []
    D12_v = []
    D16_v = []
    D22_v = []
    D26_v = []
    D66_v = []
    i = 1
    for j in range(0,myPlyNumber,1):
        D11_v.append((1/3.0)*myQ11_S[j]*(Z[i]**3-Z[i-1]**3))
        D12_v.append((1/3.0)*myQ12_S[j]*(Z[i]**3-Z[i-1]**3))
        D16_v.append((1/3.0)*myQ16_S[j]*(Z[i]**3-Z[i-1]**3))
        D22_v.append((1/3.0)*myQ22_S[j]*(Z[i]**3-Z[i-1]**3))
        D26_v.append((1/3.0)*myQ26_S[j]*(Z[i]**3-Z[i-1]**3))
        D66_v.append((1/3.0)*myQ66_S[j]*(Z[i]**3-Z[i-1]**3))
        i = i+1
    
    
    D11 = sum(D11_v)
    D12 = sum(D12_v)
    D16 = sum(D16_v)
    D22 = sum(D22_v)
    D26 = sum(D26_v)
    D66 = sum(D66_v)
    return A11,A12,A16,A22,A26,A66,B11,B12,B16,B22,B26,B66,D11,D12,D16,D22,D26,D66

#------------------------------------------------------------------------------ 
#------------------------------------------------------------------------------ 
#
# section for functions - END
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# Name of the Finite Element Analysis Modell    
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------    

# name of the numerical model

myName = ['S0_Yield']

# Limit for the outer loop (number of shells investigated)

Limit = 1

#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# Material    
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Layup and number of plies

myLaminate1 = [45,-45,0,90,90,0,-45,45]
myLayerNumber_v = [len(myLaminate1)]
myShell = [myLaminate1]

# Material parameter

myE1_v = [208000]
myE2_v = [208000]
myG12_v = [80000]
myNu12_v = [0.3]
myG23_v = [80000]

# ply thickness

myLaminateThickness_v = [0.254/8*15]
myCore_v = [12.7]

#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# Gemetry   
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# Cylinder length and radius

myHeight_v = [126.48*4]                  
myRadius_v = [126.48]
mySemi_Vertex_Angle = 0.0


#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# Finite Element Analysis  
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# element length
# can be estimated with = 0.5*np.sqrt(myRadius*myLayerNumber*myLaminateThickness)

myMesh_Size_v = [0.5*np.sqrt(myRadius_v[0]*myLayerNumber_v[0]*myLaminateThickness_v[0])]
myMesh_Size_v = [4.5]
# axial shortening / rotation of the cylinder (simulation is displacement controlled)

my_disp_v = [0.5]

# number of cores for the simulation

myCpu = 4

#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------
#
# LRSM Analysis 
#
#------------------------------------------------------------------------------   
#------------------------------------------------------------------------------
#------------------------------------------------------------------------------

# components of ABD Stiffness matrix

A11,A12,A13,A22,A23,A33,B11,B12,B13,B22,B23,B33,D11,D12,D13,D22,D23,D33 = CLT(myE1_v[0],myE2_v[0],myG12_v[0],myNu12_v[0],myLaminate1,myLaminateThickness_v[0],len(myLaminate1))

# membrane stiffness reduction factor

LRSM_Factor = 1000

# start and end of iterations for LRSM (20 increments)
    
my_START = [1]
my_END = [21]

# outer loop for different shells with vectors for result extraction
 
for jc in range(0,1,1):
  somefloats_f = []
  somefloats_f2 = []
  somefloats_f3 = []#
  somefloats_u = []
  pert_L_v = []
  
  # inner loop for LRSM iteration
  
  for ic in range(my_START[0],my_END[0],1):
    # Define the material variables   
    myHeight = myHeight_v[jc]    # Cylinder length [mm]                                                                
    myRadius = myRadius_v[0]    # Cylinder radius [mm]                                                              
    myLaminateThickness = myLaminateThickness_v[0]   # Laminate Thickness [mm]
    myCore = myCore_v[0]
    myE1 = myE1_v[0]   # Laminate E1 [MPa]                                                                     
    myE2 = myE2_v[0] # Laminate E2 [MPa]                                                                           
    myG12 = myG12_v[0]    # Laminate G12 [MPa]                                                                    
    myNu12 = myNu12_v[0]    # Laminate nu12 [-]                                                                         
    myG13 = myG12   # Outplane shear modulus [MPa]                                                          
    myG23 = myG12    # Outplane shear modulus[MPa]                                                           
    myNu13 = myNu12   # Poissons coefficient [-]                                                            
    myNu23 = myNu13   # Poissons coefficient [-] 
    myLaminateDensity = 7.8e-9    # Laminate density [t/mm³]
    myLayerNumber = myLayerNumber_v[0]   # Number of lamina layers                                                          
    myLayerOrientation = myShell[0]    # Orientation lamina layers 
    disp = my_disp_v[jc]  # axial shortening
     # perturbation displacement for SPDA
    Mesh_Size = myMesh_Size_v[0] # mesh
    SlantLength = myHeight/(np.cos(mySemi_Vertex_Angle*np.pi/180))
    my_Small_Radius = myRadius-np.sin(mySemi_Vertex_Angle*np.pi/180)*SlantLength
    Average_Radius = (my_Small_Radius+myRadius)/(2*np.cos(mySemi_Vertex_Angle*np.pi/180))
    myLoadKDF_v = [Average_Radius*0.4/my_END[0]]
    myPerturbation = myLoadKDF_v[0]*ic
    L_R = int(SlantLength/Average_Radius)
    R_t = int(Average_Radius/(myLayerNumber*myLaminateThickness))
    Zb = int((SlantLength**2*(np.sqrt(1-0.3**2)))/(Average_Radius*myLayerNumber*myLaminateThickness))
    An = int(mySemi_Vertex_Angle)
    myString = str(myName[0])+'_Z_'+str(Zb)+'_ß_'+str(An)+'_LRSM_FC_Loop-'+str(ic)
    myString_result = str(myName[0])+'_Z_'+str(Zb)+'_ß_'+str(An)
    myModel = mdb.Model(name=myString)
    
    
    myID = createConeGeometry(myRadius,myHeight,mySemi_Vertex_Angle)
    
    # Create Cylinder Material
         
    createMaterial(myE1, myE2, myNu12, myG12, myG13, myG23,'CFRP_Cone')
    myModel.materials['CFRP_Cone'].Plastic(table=((235.0, 0.0), ))
    myModel.Material(name='Sandwich')
    myModel.materials['Sandwich'].Elastic(type=ENGINEERING_CONSTANTS, table=((0.344, 0.262, 413.69, 0.45, 0.0001, 0.0001, 0.12, 203.4, 82.74), ))
    
    p = myModel.parts['Part-1']
    faces = p.faces
    p.Set(faces=faces, name='Set-Cone')
    
    
    
    # Create Cylinder Bottom and Top Edge for Boundary Condition and Load Application
    
    createGeometrySet_Edge(myRadius,0.0,0.0,'Set-Top_EDGE')
    createGeometrySet_Edge(my_Small_Radius,0.0,myHeight,'Set-Bottom_EDGE')
    
    #-------------------------------------------------------------------------- 
    
    # Define material orientation
    
    #-------------------------------------------------------------------------- 
    
#    def createMaterial_CompositeLayup():
    
    region = regionToolset.Region(faces=faces)
    orientation = myModel.parts['Part-1'].datums[myID]
    myModel.parts['Part-1'].MaterialOrientation(region=region, orientationType=SYSTEM, axis=AXIS_2, localCsys=orientation, fieldName='', additionalRotationType=ROTATION_NONE, angle=0.0, additionalRotationField='')
    #: Specified material orientation has been assigned to the selected regions.
    
    #-------------------------------------------------------------------------- 
    
    # Assembly 
    
    #-------------------------------------------------------------------------- 
    
    a = myModel.rootAssembly
    p = myModel.parts['Part-1']
    a.Instance(name='Part-1-1', part=p, dependent=ON)
        
    #--------------------------------------------------------------------------    
    
    # Reference Punkte erstellen
    
    #-------------------------------------------------------------------------- 
    
    a = myModel.rootAssembly
    r = a.referencePoints
    e = a.instances['Part-1-1'].edges
    a.features['Part-2-1'].suppress()
                   
    myRefpoint1,myVertice1 = createReferencePoint(0.0,0.0,0.0)
    myRefpoint2,myVertice2 = createReferencePoint(0.0,0.0,myHeight)
    #-------------------------------------------------------------------------- 
    
    # Create the rigid body for the top of the cylinder
    
    #-------------------------------------------------------------------------- 
    myEdgeTop = e.findAt((myRadius,0,0),)
    createRigidBodyInteraction(myVertice1,myEdgeTop,'Constraint-Top')
    #-------------------------------------------------------------------------- 
    
    # Create the rigid body for the bottom of the cylinder
    
    #-------------------------------------------------------------------------- 
    myEdgeBottom = e.findAt((my_Small_Radius,0,myHeight),)
    createRigidBodyInteraction(myVertice2,myEdgeBottom,'Constraint-Bottom')
    #-------------------------------------------------------------------------- 
    
    # Create the simulation step 
    
    #-------------------------------------------------------------------------- 
    
    createStaticStep('Step-1','Initial',300,0.001,0.001,1E-05)
    
    #-------------------------------------------------------------------------- 
    
    # Apply the loading (displacement controlled)
    
    #-------------------------------------------------------------------------- 
    
#   uncomment / comment for axial compression    
    #-------------------------------------------------------------------------- 
    #-------------------------------------------------------------------------- 
#    createDispBC_Edge_RP(myVertice1,'RP-1','BC Top','Initial',0,0,UNSET,0,0,0)
#    createDispBC_Edge_RP(myVertice2,'RP-2','BC Bottom','Initial',0,0,0,0,0,0)
#    createDispBC_Edge_RP(myVertice1,'RP-1','Disp','Step-1',UNSET,UNSET,disp,UNSET,UNSET,UNSET)
    
#   uncomment / comment for pure bending
    #-------------------------------------------------------------------------- 
    #--------------------------------------------------------------------------
    createDispBC_Edge_RP(myVertice1,'RP-1','BC Top','Initial',0,0,UNSET,UNSET,0,0)
    createDispBC_Edge_RP(myVertice2,'RP-2','BC Bottom','Initial',0,0,UNSET,UNSET,0,0)
    
    createDispBC_Edge_RP(myVertice1,'RP-1','Bend_RP-1','Step-1',UNSET,UNSET,UNSET,disp,UNSET,UNSET)
    createDispBC_Edge_RP(myVertice2,'RP-2','Bend_RP-2','Step-1',UNSET,UNSET,UNSET,-disp,UNSET,UNSET)
    
    #-------------------------------------------------------------------------- 
    
    # Define Datum plane and Edge for Perturbation application
    myID_YZ = createPlane_by_Prinicpal(YZPLANE,0.0)
    myID_XY = createPlane_by_Prinicpal(XYPLANE,myHeight/2.0)
    myID_XZ = createPlane_by_Prinicpal(XZPLANE,0)
    ##
    create_LRSM_Surface(myPerturbation)
    createGeometrySet_Face(0,(myRadius+my_Small_Radius)/2.0,myHeight/2.0,'Set-LRSM-reduced-stiffness-surface')
    createFace_Partition_byPlane(myID_YZ)
    createEdge_Set(0,(myRadius+my_Small_Radius)/2.0,myHeight/2.0,'TOP_EDGE')
    a=myModel.parts['Part-1']
    a.SetByBoolean(name='Set-Cone_no_LRSM', operation=DIFFERENCE, sets=(a.sets['Set-Cone'], a.sets['Set-LRSM-reduced-stiffness-surface'], ))
    createFace_Partition_byPlane(myID_XY)
    myID_XY = createPlane_by_Prinicpal(XYPLANE,myHeight/2.0+myPerturbation+1*Mesh_Size)
    createFace_Partition_byPlane(myID_XY)
    myID_XY = createPlane_by_Prinicpal(XYPLANE,myHeight/2.0-myPerturbation-1*Mesh_Size)
    createFace_Partition_byPlane(myID_XY)
    createFace_Partition_byPlane(myID_XZ)
    myID_YZ = createPlane_by_Prinicpal(YZPLANE,-myPerturbation-2*Mesh_Size)
    createFace_Partition_byPlane(myID_YZ)
    myID_YZ = createPlane_by_Prinicpal(YZPLANE,myPerturbation+2*Mesh_Size)
    createFace_Partition_byPlane(myID_YZ)
#    
    #--------------------------------------------------------------------------
    # Define LRSM surface and LRSM reference surface
    create_GeneralStiffness('Section-Reference',A11, A12, A22, A13, A23, A33, B11, B12, B13, D11, B12, B22, B23, D12, D22, B13, B23, B33, D13, D23, D33)
    create_GeneralStiffness('Section-Reduced_Stiffness_LRSM',A11/LRSM_Factor, A12/LRSM_Factor, A22/LRSM_Factor, A13/LRSM_Factor, A23/LRSM_Factor, A33/LRSM_Factor, 0, 0, 0, D11, 0, 0, 0, D12, D22, 0, 0, 0, D13, D23, D33)
#   uncomment / comment for cutout
#-------------------------------------------------------------------------- 
#-------------------------------------------------------------------------- 
#    p = myModel.parts['Part-1']
#    f = p.faces
#    p.RemoveFaces(faceList = f[28:30]+f[32:33]+f[35:36], deleteCells=False)
    #--------------------------------------------------------------------------
    
    # create FEA mesh
    
    createMesh('Part-1',Mesh_Size)
    
    #--------------------------------------------------------------------------
    
    # create material section
        
    #--------------------------------------------------------------------------
    a = myModel.rootAssembly
    region = p.sets['Set-LRSM-reduced-stiffness-surface']
    p.SectionAssignment(region=region, sectionName='Section-Reduced_Stiffness_LRSM', offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
    layupOrientation = myModel.parts['Part-1'].datums[myID]
    p = myModel.parts['Part-1']
    f = p.faces
    region = p.sets['Set-Cone_no_LRSM']
    myRegion=region
    
    compositeLayup = myModel.parts['Part-1'].CompositeLayup(name='CompositeLayup-1', description='', elementType=SHELL, offsetType=MIDDLE_SURFACE, symmetric=False, thicknessAssignment=FROM_SECTION)
    
    compositeLayup.Section(preIntegrate=OFF, integrationRule=SIMPSON, thicknessType=UNIFORM, poissonDefinition=DEFAULT, temperature=GRADIENT, useDensity=OFF)
    
    compositeLayup.ReferenceOrientation(orientationType=SYSTEM, localCsys=layupOrientation, fieldName='', additionalRotationType=ROTATION_NONE, angle=0.0, additionalRotationField='', axis=AXIS_2)    
    
    for j in range(0,myLayerNumber):
        if (myLayerOrientation[j] == 1):
            compositeLayup.CompositePly(suppressed=False, plyName='myPly-'+str(j), region=myRegion, material='Sandwich', thicknessType=SPECIFY_THICKNESS, thickness=myCore, orientationType=SPECIFY_ORIENT, orientationValue=myLayerOrientation[j])
        else:
            compositeLayup.CompositePly(suppressed=False, plyName='myPly-'+str(j), region=myRegion, material='CFRP_Cone', thicknessType=SPECIFY_THICKNESS, thickness=myLaminateThickness, orientationType=SPECIFY_ORIENT, orientationValue=myLayerOrientation[j])
    
    #-------------------------------------------------------------------------- 
            
    # Generate Job name and submit Job
    
    #--------------------------------------------------------------------------   
    
    myJobName = myString
    myJob = mdb.Job(name=myJobName, model=myJobName, description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', multiprocessingMode=DEFAULT, numCpus=myCpu, numDomains=myCpu, numGPUs=0)
    myJob.submit()
    myJob.waitForCompletion()
#    
    myPerturbation = [myPerturbation/Average_Radius]
    pert_L_v = pert_L_v + myPerturbation
    odb = session.openOdb(str(myString)+'.odb')
    time_v = []
    Force_v = []
    AxialDisp_v = []
    lastStep=odb.steps['Step-1']
    
    for x in range(len(lastStep.frames)):
        lastFrame = lastStep.frames[x]
        Time1=lastFrame.frameValue
        rForce = lastFrame.fieldOutputs['RM']
        center = odb.rootAssembly.nodeSets['RP-1']
        centerRForce = rForce.getSubset(region=center)
        time_vr = [Time1]
        time_v = time_v + time_vr
        
        
        
        rDisp = lastFrame.fieldOutputs['UR']
        center = odb.rootAssembly.nodeSets['RP-1']
        centerRDisp = rDisp.getSubset(region=center)
    
        for i in centerRDisp.values:
            AxialDisp_vr = [i.data[0]]
            AxialDisp_v = AxialDisp_v + AxialDisp_vr
    
        for i in centerRForce.values:
            Force_vr = [i.data[0]]
            Force_v = Force_v + Force_vr
#    
#    
  
    Force = np.max(Force_v) 
    somefloats_fr = [Force]
    somefloats_f = somefloats_f + somefloats_fr
    Index = np.argmax(Force_v)  
    Disp = AxialDisp_v[Index]
    somefloats_fr3 = [Disp]
    somefloats_f3 = somefloats_f3 + somefloats_fr3
    np.savetxt('N_'+str(myString_result)+'.txt',somefloats_f)
    np.savetxt('u_'+str(myString_result)+'.txt',somefloats_f3)
    np.savetxt('v_'+str(myString_result)+'.txt',pert_L_v)
    np.savetxt('Rf_'+str(myString)+'.txt',Force_v)
    np.savetxt('u_'+str(myString)+'.txt',AxialDisp_v)
  
#  