#######################################
# Knickstab VSG Parameterrechnung#
#######################################



"""
Dieses Skript...
... modelliert einen Stab
... führt eine Buckling Analyse durch
... verschiebt eine Ecke aus der Ebene heraus
... und belastet optional um das Eigengewicht

"""

hoehe = 2000.0		#x Ausdehnung
breite = 250.0		#y Ausdehnung
dglas = [10, 10]			#z Ausdehnung Glas
dfolie = [1.52]			#z Ausdehnung Folie
dicke = sum(dglas) + sum(dfolie)       #z Ausdehnung Gesamt


meshgroesse = 50
anzElemDickeFolie = 2    #Langosch orientieren später Konervenzstudie durchführen
anzElemDickeGlas = 1    #Langosch orientieren später Konervenzstudie durchführen pro Halbe Scheibendicke
emodulGlas = 70000
mueGlas = 0.23
SchubmodulFolie = 1		# Langosch G[0;1;10;∞]
mueFolie = 0.39
emodulFolie = 2*(1+mueFolie)*SchubmodulFolie

mInc = 0.05		#maximale Inkrementgröße
initInc = 0.05   #Initiale Inkrementgröße
workspace= 'C:\Users\my-yi\Desktop\MA Glasbau\Versuche'
beElemType='C3D8R'

#Buckle Step bezogen
eNull = hoehe/400          #Amplitude der Imperfektion
numEigForm = 4     #Anzahl der angefragten Eigenformen
eigForm = 1        #ausgewählte Eigenform


###########################################################################################################################################
########################### SKRIPT BODY ###################################################################################################
###########################################################################################################################################

modelname = 'Stuetze_'+str(dglas[0])

for i in range(1,len(dglas)):
    modelname=modelname + '-' + str(int(dfolie[i-1]/0.38))+'f'
    modelname=modelname + '-' + str(int(dglas[i]))

modelname = modelname + '_G'+str(SchubmodulFolie)


partname = 'Part_'+modelname
delta = dicke/anzElemDickeGlas/10

# Do not delete the following import lines
from abaqus import *
from abaqusConstants import *
import __main__
import section
import regionToolset
import displayGroupMdbToolset as dgm
import part
import material
import assembly
import step
import interaction
import load
import mesh
import optimization
import job
import sketch
import visualization
import xyPlot
import displayGroupOdbToolset as dgo
import connectorBehavior
import os

from odbAccess import *
import numpy as np
#############################################################################

#!!!!!!!!!!!!!!WIIIIIICHTIG!!!!!!!!!!!!!!!!!!!
#                   | 
#                   |
#                   V   
#Anpassen des Workdirectorys
os.chdir(r"C:\Users\my-yi\Desktop\MA Glasbau\Versuche")#


# Erstellen des Modells
mdb.Model(name=modelname, modelType=STANDARD_EXPLICIT)
a = mdb.models[modelname].rootAssembly

# Erstellen des Parts
s1 = mdb.models[modelname].ConstrainedSketch(name='__profile__', 
	sheetSize=200.0)
s1.rectangle(point1=(0, 0), point2=(breite, hoehe))
p = mdb.models[modelname].Part(name=partname, dimensionality=THREE_D, 
	type=DEFORMABLE_BODY)
p.BaseSolidExtrude(sketch=s1, depth=dicke)
p = mdb.models[modelname].parts[partname]
del mdb.models[modelname].sketches['__profile__']

# Erstellen der Datumsebenen relativ zu den Flächen mit Berücksichtigung der Dicken
f = p.faces.findAt((delta,delta,0),)


#2Scheiben und 3Scheiben (erstellt für 4 und 5 Scheiben die letzten 2 planes nicht)
datumDicken = []
datumDicken.append(dglas[0])

for i in range(1, len(dglas)):
    if i % 2 != 0:
        datumDicken.append(datumDicken[-1] + dfolie[i - 1])
    else:
        datumDicken.append(datumDicken[-1] + dglas[i - 1])

if len(dglas) % 2 != 0:
    datumDicken.append(datumDicken[-1] + dfolie[-1])
    datumDicken.append(datumDicken[-1] + dglas[-1])
else:
    datumDicken.append(datumDicken[-1] + dfolie[-1])

for i in range(len(datumDicken)-1):
    p.DatumPlaneByOffset(plane=f, flip=SIDE2, offset=float(datumDicken[i]))


"""
#4 Scheiben (erstellt 1 Plane am Ende des Parts und 1 Plane dahinter zu viel)
datumDicken = []
datumDicken.append(dglas[0])

for i in range(1, len(dglas)):
    if i % 2 != 0:
        datumDicken.append(datumDicken[-1] + dfolie[i - 1])
    else:
        datumDicken.append(datumDicken[-1] + dglas[i - 1])

for i in range(len(dglas)):
    if i % 2 != 0:
        datumDicken.append(datumDicken[-1] + dfolie[-1])
    else:
        datumDicken.append(datumDicken[-1] + dglas[-1])

for i in range(len(datumDicken)-1):
    p.DatumPlaneByOffset(plane=f, flip=SIDE2, offset=float(datumDicken[i]))

"""
# Erstellen der Partition basierend auf den zuvor erstellten Datumsebenen
for i in range(len(datumDicken)-1):
    c = p.cells
    # Partitionierung der Zellen mit der ersten Datumsebene (f[4])
    pickedCells = c.findAt((0,0,datumDicken[i]-0.01),)
    d1 = p.datums
    p.PartitionCellByDatumPlane(datumPlane=d1[i+2], cells=pickedCells)

#Erstellen einer Partition im Schwerpunkt der Glasscheiben zur Einleitung der Last im Fall ohne Schlupfendbehinderung
for i in range(len(dglas)):
    j=i*2
    #Partitions an Unterseite
    f = p.faces
    pickedFaces = f.findAt((breite/2,0,(datumDicken[j]-0.01)),)
    v, e, d = p.vertices, p.edges, p.datums
    p.PartitionFaceByShortestPath(faces=pickedFaces, point1=p.InterestingPoint(
        edge=e.findAt((0,0,datumDicken[j]-0.01),), rule=MIDDLE), point2=p.InterestingPoint(edge=e.findAt((breite,0,datumDicken[j]-0.01),), 
        rule=MIDDLE))
    #Partitions an Oberseite
    pickedFaces = f.findAt((breite/2,hoehe,(datumDicken[j]-0.01)),)
    v, e, d = p.vertices, p.edges, p.datums
    p.PartitionFaceByShortestPath(faces=pickedFaces, point1=p.InterestingPoint(
        edge=e.findAt((0,hoehe,datumDicken[j]-0.01),), rule=MIDDLE), point2=p.InterestingPoint(edge=e.findAt((breite,hoehe,datumDicken[j]-0.01),), 
        rule=MIDDLE))


# Material erstellen Glas
mdb.models[modelname].Material(name='Glas')
mdb.models[modelname].materials['Glas'].Elastic(table=((emodulGlas, mueGlas), ))

# Material erstellen Folie elastischer Schubmodul
mdb.models[modelname].Material(name='Folie')
mdb.models[modelname].materials['Folie'].Elastic(table=((emodulFolie, mueFolie), ))

#Section Erstellung
mdb.models[modelname].HomogeneousSolidSection(name='GlasSection', 
	material='Glas', thickness=None)

mdb.models[modelname].HomogeneousSolidSection(name='FolieSection', 
material='Folie', thickness=None)

#Section Zuweisung
#Zuweisung der Glassection
for i in range(len(dglas)):
    j=i*2
    c = p.cells
    c1 = c.findAt((0,0,datumDicken[j]-0.01), )
    q = c1.index
    cells = c[q:q+1]
    region = p.Set(cells=cells, name='Scheibe-'+str(i))
    p.SectionAssignment(region=region, sectionName='GlasSection', offset=0.0, 
    offsetType=MIDDLE_SURFACE, offsetField='', 
    thicknessAssignment=FROM_SECTION)

#Zuweisung der Foliensection
for i in range(len(dfolie)):
    j=i*2
    c = p.cells
    c1 = c.findAt((0,0,datumDicken[j]+0.01), )
    q = c1.index
    cells = c[q:q+1]
    region = p.Set(cells=cells, name='Folie-'+str(i))
    p.SectionAssignment(region=region, sectionName='FolieSection', offset=0.0, 
    offsetType=MIDDLE_SURFACE, offsetField='', 
    thicknessAssignment=FROM_SECTION)
    

#Instance Erstellen
a.DatumCsysByDefault(CARTESIAN)
a.Instance(name=partname, part=p, dependent=ON)

#Reference Point erstellen und auslesen der ID zur Ansprache in Coupling
a.ReferencePoint(point=(breite/2, hoehe+100, dicke/2))
idRF1 = mdb.models[modelname].rootAssembly.referencePoints.values()[0]

a.ReferencePoint(point=(breite/2, -100, dicke/2))
idRF2 = mdb.models[modelname].rootAssembly.referencePoints.values()[0]

#Refpoint umbenennen
mdb.models[modelname].rootAssembly.features.changeKey(fromName='RP-1', 
    toName='Loslager+Lastpunkt')
mdb.models[modelname].rootAssembly.features.changeKey(fromName='RP-2', 
    toName='Festlager')
"""
# 3 Scheiben
#Coupling Reference Point mit Glasmitten oben bei 2 Scheiben Fehler bei edgeCurr

refPoints1=(idRF1, )
region1=a.Set(referencePoints=refPoints1, name='Loslager+Lastpunkt')

edgesTop1 = p.edges.findAt((breite/2,hoehe,dglas[0]/2),)
for i in range(0,len(dglas)):
    j=i*2
    p = mdb.models[modelname].parts[partname]
    s = p.edges
    edgeCurr = s.getByBoundingBox(-1,hoehe-1,datumDicken[j]-dglas[i]/2-1,breite+1,hoehe+1,datumDicken[j]-dglas[i]/2+1)
    idx = edgeCurr[0].index
    p.Set(edges=s[idx:idx+1], name='slaveTop'+str(i+1))

inst = a.instances[partname]
combined_edgesTop = inst.edges[0:0]
for i in range(0, len(dglas)):
    edge_set = inst.sets['slaveTop' + str(i + 1)]
    combined_edgesTop += edge_set.edges  

a.Set(edges=combined_edgesTop, name='Combined_slaveTop')
region2=a.sets['Combined_slaveTop']
    
mdb.models[modelname].Coupling(name='Loslager+Last', 
    controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, 
    couplingType=KINEMATIC, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, 
    ur2=ON, ur3=ON)
"""
#Coupling Reference Point mit Glasmitten oben

refPoints1=(idRF1, )
region1=a.Set(referencePoints=refPoints1, name='Loslager+Lastpunkt')

edgesTop1 = p.edges.findAt((breite/2,hoehe,dglas[0]/2),)
for i in range(0,len(dglas)):
    j=i*2
    p = mdb.models[modelname].parts[partname]
    s = p.edges
    edgeCurr = s.getByBoundingBox(-1,hoehe-1,datumDicken[j]-dglas[i]/2-1,breite+1,hoehe+1,datumDicken[j]-dglas[i]/2+1)
    idx = edgeCurr[0].index
    p.Set(edges=s[idx:idx+1], name='slaveTop'+str(i+1))

inst = a.instances[partname]
combined_edgesTop = inst.edges[0:0]
for i in range(0, len(dglas)):
    edge_set = inst.sets['slaveTop' + str(i + 1)]
    combined_edgesTop += edge_set.edges  

a.Set(edges=combined_edgesTop, name='Combined_slaveTop')
region2=a.sets['Combined_slaveTop']
    
mdb.models[modelname].Coupling(name='Loslager+Last', 
    controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, 
    couplingType=KINEMATIC, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, 
    ur2=ON, ur3=ON)
#Coupling Reference Point mit Glasmitten unten

refPoints2=(idRF2, )
region1=a.Set(referencePoints=refPoints2, name='Festlager')

edgesBottom1 = p.edges.findAt((breite/2,0,dglas[0]/2),)
for i in range(0,len(dglas)):
    j=i*2
    p = mdb.models[modelname].parts[partname]
    s = p.edges
    edgeCurr = s.getByBoundingBox(-1,-1,datumDicken[j]-dglas[i]/2-1,breite+1,1,datumDicken[j]-dglas[i]/2+1)
    idx = edgeCurr[0].index
    p.Set(edges=s[idx:idx+1], name='slaveBottom'+str(i+1))

inst = a.instances[partname]
combined_edgesBottom = inst.edges[0:0]
for i in range(0, len(dglas)):
    edge_set = inst.sets['slaveBottom' + str(i + 1)]
    combined_edgesBottom += edge_set.edges  

a.Set(edges=combined_edgesBottom, name='Combined_slaveBottom')
region2=a.sets['Combined_slaveBottom']
    
mdb.models[modelname].Coupling(name='Festlager', 
    controlPoint=region1, surface=region2, influenceRadius=WHOLE_SURFACE, 
    couplingType=KINEMATIC, localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, 
    ur2=ON, ur3=ON)

#BC Erstellen (Euler Fall II)

#BC Festlager
mdb.models[modelname].rootAssembly.Set(name='Set-2', referencePoints=(
    mdb.models[modelname].rootAssembly.referencePoints[5], ))
mdb.models[modelname].DisplacementBC(amplitude=UNSET, createStepName=
    'Initial', distributionType=UNIFORM, fieldName='', localCsys=None, name=
    'Festlager', region=mdb.models[modelname].rootAssembly.sets['Set-2']
    , u1=SET, u2=SET, u3=SET, ur1=UNSET, ur2=SET, ur3=UNSET)

#BC Loslager

mdb.models[modelname].rootAssembly.Set(name='Set-3', referencePoints=(
    mdb.models[modelname].rootAssembly.referencePoints[4], ))
mdb.models[modelname].DisplacementBC(amplitude=UNSET, createStepName=
    'Initial', distributionType=UNIFORM, fieldName='', localCsys=None, name=
    'Loslager', region=mdb.models[modelname].rootAssembly.sets['Set-3'], 
    u1=SET, u2=UNSET, u3=SET, ur1=UNSET, ur2=SET, ur3=UNSET)

###########################################################################################################################################################################################################
# Mesh erstellen
region = (p.cells,)
elemType = mesh.ElemType(
    elemCode=C3D20R,
    elemLibrary=STANDARD,
    kinematicSplit=AVERAGE_STRAIN,
    secondOrderAccuracy=OFF,
    hourglassControl=ENHANCED,
    distortionControl=DEFAULT
)
p.setElementType(regions=region, elemTypes=(elemType,))
p.seedPart(size=meshgroesse, deviationFactor=0.1, minSizeFactor=0.1)

# Seed-Operationen für Kanten basierend auf den verschiedenen Dicken
edges = e.findAt(((0, 0, float(dfolie[0] / 2)),), ((hoehe, 0, float(dfolie[0] / 2)),),
                 ((hoehe, breite, float(dfolie[0] / 2)),), ((0, breite, float(dfolie[0] / 2)),))
p.seedEdgeByNumber(edges=edges, number=anzElemDickeFolie, constraint=FINER)

for d in dglas:
    edges = e.findAt(((0, 0, float(d / 2)),), ((hoehe, 0, float(d / 2)),),
                     ((hoehe, breite, float(d / 2)),), ((0, breite, float(d / 2)),))
    p.seedEdgeByNumber(edges=edges, number=anzElemDickeGlas, constraint=FINER)

edges = e.findAt(((0, 0, float(dicke - dfolie[0] / 2)),), ((hoehe, 0, float(dicke - dfolie[0] / 2)),),
                 ((hoehe, breite, float(dicke - dfolie[0] / 2)),), ((0, breite, float(dicke - dfolie[0] / 2)),))
p.seedEdgeByNumber(edges=edges, number=anzElemDickeGlas, constraint=FINER)

# Mesh generieren
p.generateMesh()
a.regenerate()