# -*- coding: mbcs -*-
# Do not delete the following import lines
from odbAccess import *
from abaqus import *
from abaqusConstants import *
from odbAccess import openOdb, Odb 
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 math
import numpy
import array

#### Parameter ####
## Geometrie ##
D_klein = 36.0						#kleiner Durchmesser
D_gross = (2*(D_klein)**3)**(1/3.)		#grosser Durchmesser
R_mitte = 187.368					#Raduis der Verjüngung in der Mitte
a = (D_gross - D_klein)/2			#Tiefe der Verjüngung
Sehne = 2*(a*(2*R_mitte-a))**(1/2.)	#Länge der Kreissehne = Länge des modellierten Miittelteils

## Vernetzung Hochbeanspruchtes Volumen ##
HB_t =	0.5	#Tiefe hochbeanspruchter Bereich
SA_a = 40	#Seedanzahl axial (gerade Anzahl wählen)
SA_r = 80	#Seedanzahl radial
SA_t = 20 	#Seedanzahl tangential (gerade Anzahl wählen)
SA_ra = int(round(D_gross))
Bias_r = 5
Bias_a = 1

## Geometrieparameter für die Innere Schalung
R_innen = R_mitte + HB_t
alpha_R=2*asin(Sehne/(2*R_mitte))
Sehne_innen = 2*R_innen*sin(alpha_R/2)
h_S_Si=(Sehne_innen-Sehne)/2
alpha=asin(h_S_Si/(R_innen-R_mitte))
x=cos(alpha)*(R_innen-R_mitte)


## Vernetzung Restgeometrie
SA_r_i = SA_t/2+1	#Seedanzahl radial (Anders geht's innen nicht auf von der Vergröberung)

Modell = mdb.models['Model-1']

#### Probengeometrie aufbauen ####

sketch_ungekerbteProbe = Modell.ConstrainedSketch(name='ungekerbteProbe', sheetSize=200.0)		#Zeichenebene											
g, v, d, c = sketch_ungekerbteProbe.geometry, sketch_ungekerbteProbe.vertices, sketch_ungekerbteProbe.dimensions, sketch_ungekerbteProbe.constraints		
#2
sketch_ungekerbteProbe.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))	#Drehachse 
sketch_ungekerbteProbe.FixedConstraint(entity=g[2])
#3
sketch_ungekerbteProbe.ConstructionLine(point1=(-100.0, 0.0), point2=(100.0, 0.0))
sketch_ungekerbteProbe.HorizontalConstraint(entity=g[3], addUndoState=False)
#4
sketch_ungekerbteProbe.Arc3Points(point1=(-Sehne/2, D_gross/2), point2=(Sehne/2, D_gross/2), point3=(0.0,D_klein/2))
#5
sketch_ungekerbteProbe.Line(point1=(-Sehne/2, D_gross/2), point2=((-Sehne/2)-h_S_Si, (D_gross/2)-x))						#Verbindung AC
#6
sketch_ungekerbteProbe.Line(point1=((-Sehne/2)-h_S_Si, (D_gross/2)-x), point2=((-Sehne/2)-h_S_Si,0))
#7
sketch_ungekerbteProbe.Line(point1=((-Sehne/2)-h_S_Si, 0.0), point2=(Sehne/2+h_S_Si, 0.0))
#8
sketch_ungekerbteProbe.Line(point1=(Sehne/2+h_S_Si, 0.0), point2=(Sehne/2+h_S_Si, (D_gross/2)-x))							
#9
sketch_ungekerbteProbe.Line(point1=(Sehne/2+h_S_Si, (D_gross/2)-x),point2=(Sehne/2, D_gross/2))							#Verbindung CA


sketch_ungekerbteProbe.sketchOptions.setValues(constructionGeometry=ON)
sketch_ungekerbteProbe.assignCenterline(line=g[3])
part_ungekerbteProbe = mdb.models['Model-1'].Part(name='ungekerbteProbe', dimensionality=THREE_D, 
        type=DEFORMABLE_BODY)
part_ungekerbteProbe.BaseSolidRevolve(sketch=sketch_ungekerbteProbe, angle=360.0, flipRevolveDirection=OFF)

sketch_ungekerbteProbe_1 = Modell.ConstrainedSketch(name='ungekerbteProbe_1', sheetSize=200.0)		#Zeichenebene											
g, v, d, c = sketch_ungekerbteProbe_1.geometry, sketch_ungekerbteProbe_1.vertices, sketch_ungekerbteProbe_1.dimensions, sketch_ungekerbteProbe_1.constraints		
sketch_ungekerbteProbe_1.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))	#Drehachse 
sketch_ungekerbteProbe_1.FixedConstraint(entity=g[2])
sketch_ungekerbteProbe_1.ConstructionLine(point1=(-100.0, 0.0), point2=(100.0, 0.0))
sketch_ungekerbteProbe_1.HorizontalConstraint(entity=g[3], addUndoState=False)
sketch_ungekerbteProbe_1.Arc3Points(point1=(-Sehne/2, D_gross/2), point2=(Sehne/2, D_gross/2), point3=(0.0,D_klein/2))
sketch_ungekerbteProbe_1.Line(point1=(-Sehne/2, D_gross/2), point2=((-Sehne/2)-h_S_Si, (D_gross/2)-x))						#Verbindung AC
#sketch_ungekerbteProbe_1.Line(point1=((-Sehne/2)-h_S_Si, (D_gross/2)-x), point2=((-Sehne/2)-h_S_Si,0))
#sketch_ungekerbteProbe_1.Line(point1=((-Sehne/2)-h_S_Si, 0.0), point2=(Sehne/2+h_S_Si, 0.0))
#sketch_ungekerbteProbe_1.Line(point1=(Sehne/2+h_S_Si, 0.0), point2=(Sehne/2+h_S_Si, (D_gross/2)-x))							
sketch_ungekerbteProbe_1.Line(point1=(Sehne/2+h_S_Si, (D_gross/2)-x),point2=(Sehne/2, D_gross/2))							#Verbindung CA
sketch_ungekerbteProbe_1.Arc3Points(point1=((-Sehne/2)-h_S_Si, (D_gross/2)-x), point2=((Sehne/2)+h_S_Si, (D_gross/2)-x), point3=(0.0,(D_klein/2)-(R_innen-R_mitte)))


#### Set erstellen ####

cell_ges = part_ungekerbteProbe.cells
part_ungekerbteProbe.Set(cells=cell_ges, name='Probe')

SUE=part_ungekerbteProbe.edges.findAt(((-Sehne/2)-h_S_Si,D_gross/2-x,0),)
sketch_face=part_ungekerbteProbe.faces.findAt(-Sehne/2-h_S_Si,0,0)

#### Partitionieren ####

# EbenenErstellung zur Abgrenzung des Hochbelasteten Bereichs # 
Ebene_w1 = part_ungekerbteProbe.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=Sehne/4)
Ebene_w1 = part_ungekerbteProbe.datums[Ebene_w1.id]
Ebene_w2 = part_ungekerbteProbe.DatumPlaneByPrincipalPlane(principalPlane=YZPLANE, offset=-Sehne/4)
Ebene_w2 = part_ungekerbteProbe.datums[Ebene_w2.id]
Ebene_w3 = part_ungekerbteProbe.DatumPlaneByPrincipalPlane(principalPlane=XYPLANE, offset=0.0)
Ebene_w3 = part_ungekerbteProbe.datums[Ebene_w3.id]

# Partitionieren der Mantelfläche
f, e, d2 = part_ungekerbteProbe.faces, part_ungekerbteProbe.edges, part_ungekerbteProbe.datums
t = part_ungekerbteProbe.MakeSketchTransform(sketchPlane=Ebene_w3, sketchUpEdge=d2[1], 
        sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', 
        sheetSize=200, gridSpacing=5., transform=t)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
part_ungekerbteProbe = mdb.models['Model-1'].parts['ungekerbteProbe']
part_ungekerbteProbe.projectReferencesOntoSketch(sketch=s, filter=COPLANAR_EDGES)
s.retrieveSketch(sketch=mdb.models['Model-1'].sketches['ungekerbteProbe_1'])
s.rotate(centerPoint=(0.0, 0.0), angle=90.0, objectList=(g[8], g[9], g[10], 
        g[11], g[12], g[13]))
cell=part_ungekerbteProbe.cells
f1, e1, d1 = part_ungekerbteProbe.faces, part_ungekerbteProbe.edges, part_ungekerbteProbe.datums
part_ungekerbteProbe.PartitionCellBySketch(sketchPlane=Ebene_w3, sketchUpEdge=d1[1], 
        cells=cell, sketch=s)

e = part_ungekerbteProbe.edges
cell_ges=part_ungekerbteProbe.cells
#pickedEdges=part_ungekerbteProbe.edges.findAt((0.0, D_klein/2-0.5, 0.0),)
#pickedSweep=part_ungekerbteProbe.edges.findAt((Sehne/2+h_S_Si, -D_gross/2+x, 0.0),)

SUE=e.findAt(((0.0, D_klein/2-0.5, 0.0),),)
Sweep=e.findAt(((Sehne/2+h_S_Si, -D_gross/2+x, 0.0),),)
pickedEdges=(e[1],)
part_ungekerbteProbe.PartitionCellBySweepEdge(sweepPath=Sweep[0], cells=cell_ges, edges=pickedEdges)

# Partitionieren der Stirnflächen
f, e, d1 = part_ungekerbteProbe.faces, part_ungekerbteProbe.edges, part_ungekerbteProbe.datums
SUE=e.findAt(((-Sehne/2-h_S_Si,D_gross/2-x,0),),)
sketch_face=f.findAt(((-Sehne/2-h_S_Si,0.0,0.0),),)
t = part_ungekerbteProbe.MakeSketchTransform(sketchPlane=sketch_face[0], sketchUpEdge=SUE[0], 
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
s1 = mdb.models['Model-1'].ConstrainedSketch(name='sketch_01', 
        sheetSize=200, gridSpacing=5., transform=t)
s1.rectangle(point1=(-D_klein/4, -D_klein/4), point2=(D_klein/4,D_klein/4))
s1.Line(point1=(D_klein/4, D_klein/4), point2=((D_gross/2-x)*sin(45*pi/180), (D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(-D_klein/4, D_klein/4), point2=(-(D_gross/2-x)*sin(45*pi/180), (D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(-D_klein/4, -D_klein/4), point2=(-(D_gross/2-x)*sin(45*pi/180), -(D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(D_klein/4, -D_klein/4), point2=((D_gross/2-x)*sin(45*pi/180), -(D_gross/2-x)*sin(45*pi/180)))
Face_1=f.findAt(((-Sehne/2-h_S_Si,0.0,0.0),),)
SUE=e.findAt(((-Sehne/2-h_S_Si,D_gross/2-x,0),),)
part_ungekerbteProbe.PartitionFaceBySketch(sketchUpEdge=SUE[0], faces=Face_1[0], sketch=s1)


f, e, d1 = part_ungekerbteProbe.faces, part_ungekerbteProbe.edges, part_ungekerbteProbe.datums
SUE=e.findAt(((Sehne/2+h_S_Si,D_gross/2-x,0),),)
sketch_face=f.findAt(((Sehne/2+h_S_Si,0.0,0.0),),)
t = part_ungekerbteProbe.MakeSketchTransform(sketchPlane=sketch_face[0], sketchUpEdge=SUE[0], 
    sketchPlaneSide=SIDE1, origin=(0.0, 0.0, 0.0))
s1 = mdb.models['Model-1'].ConstrainedSketch(name='sketch_01', 
        sheetSize=200, gridSpacing=5., transform=t)
s1.rectangle(point1=(-D_klein/4, -D_klein/4), point2=(D_klein/4,D_klein/4))
s1.Line(point1=(D_klein/4, D_klein/4), point2=((D_gross/2-x)*sin(45*pi/180), (D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(-D_klein/4, D_klein/4), point2=(-(D_gross/2-x)*sin(45*pi/180), (D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(-D_klein/4, -D_klein/4), point2=(-(D_gross/2-x)*sin(45*pi/180), -(D_gross/2-x)*sin(45*pi/180)))
s1.Line(point1=(D_klein/4, -D_klein/4), point2=((D_gross/2-x)*sin(45*pi/180), -(D_gross/2-x)*sin(45*pi/180)))
Face_1=f.findAt(((Sehne/2+h_S_Si,0.0,0.0),),)
SUE=e.findAt(((Sehne/2+h_S_Si,D_gross/2-x,0),),)
part_ungekerbteProbe.PartitionFaceBySketch(sketchUpEdge=SUE[0], faces=Face_1[0], sketch=s1)

pickedCells = part_ungekerbteProbe.cells
e, d1 = part_ungekerbteProbe.edges, part_ungekerbteProbe.datums
pickedEdges=(e[3], e[6], e[7], e[10])
part_ungekerbteProbe.PartitionCellByExtrudeEdge(line=d1[1], cells=pickedCells, edges=pickedEdges, sense=REVERSE)

# Partitionieren in 3 Teile # 
cell_ges = part_ungekerbteProbe.cells
part_ungekerbteProbe.PartitionCellByDatumPlane(datumPlane=Ebene_w1, cells=cell_ges)
cell_ges = part_ungekerbteProbe.cells
part_ungekerbteProbe.PartitionCellByDatumPlane(datumPlane=Ebene_w2, cells=cell_ges)

#### Material, Section, SectionAssignment ####

Modell.Material(name='S42CrMo4 QT').Elastic(table=((210000.0, 0.3), ))												#Material Definition Material
																												#Elastic elastisches Verhalten E-Modul; Poisonzahl
Modell.HomogeneousSolidSection(name='Probe', 														#HomogeneousSolidSection Erstellung Section aus Material von oben
	material='S42CrMo4 QT', thickness=None)

region = part_ungekerbteProbe.sets['Probe']
part_ungekerbteProbe.SectionAssignment(region=region, sectionName='Probe', 									#SectionAssignment Zuweisung Section der Parts
	offset=0.0, offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
	
## Material Orientierung zuweisen
f, e, c = part_ungekerbteProbe.faces, part_ungekerbteProbe.edges, part_ungekerbteProbe.cells

#Mantelfläche in der Mitte
cell=c.findAt(((0.0,D_klein/2, 0.0),),)
region=regionToolset.Region(cells=cell)
Testregion=part_ungekerbteProbe.Set(cells=cell, name='Testregion')
side1Face=f.findAt(((0.0, 0.0 , D_klein/2),),)
normalAxisRegion=part_ungekerbteProbe.Surface(side1Faces=side1Face, name='Mantel_Mitte')
edge=e.findAt(((0.0,D_klein/2, 0.0),),)
primaryAxisRegion=part_ungekerbteProbe.Set(edges=edge, name='axialeKante_Mitte')

part_ungekerbteProbe.MaterialOrientation(region=region, orientationType=DISCRETE, axis=AXIS_1,
    normalAxisDefinition=SURFACE, normalAxisRegion=normalAxisRegion, flipNormalDirection=False, 
    normalAxisDirection=AXIS_1, primaryAxisDefinition=EDGE, primaryAxisRegion=primaryAxisRegion,
    primaryAxisDirection=AXIS_2, flipPrimaryDirection=False, additionalRotationType=ROTATION_NONE, 
    angle=0.0, additionalRotationField='', stackDirection=STACK_3)
	
# Mantelflächen an den Seiten
cellS1=c.findAt(((-3*Sehne/8, 0.0, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16)),),)
region=regionToolset.Region(cells=cellS1)
Testregion=part_ungekerbteProbe.Set(cells=cellS1, name='cellS1')
side2Face=f.getClosest(coordinates=(((-3*Sehne/8, 0.0, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16))),))
side2Facev=f.findAt((((side2Face[0][1])),))
normalAxisRegion=part_ungekerbteProbe.Surface(side1Faces=side2Facev, name='Mantel_Außen_1')
primaryAxisRegion_1=e.getClosest(coordinates=(((-3*Sehne/8, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16),0.0)),))
primaryAxisRegionv=e.findAt((((primaryAxisRegion_1[0][1])),))
primaryAxisRegion=part_ungekerbteProbe.Set(edges=primaryAxisRegionv, name='axialeKante_Außen_1')

part_ungekerbteProbe.MaterialOrientation(region=region, orientationType=DISCRETE, axis=AXIS_1,
    normalAxisDefinition=SURFACE, normalAxisRegion=normalAxisRegion, flipNormalDirection=False, 
    normalAxisDirection=AXIS_1, primaryAxisDefinition=EDGE, primaryAxisRegion=primaryAxisRegion,
    primaryAxisDirection=AXIS_2, flipPrimaryDirection=False, additionalRotationType=ROTATION_NONE, 
    angle=0.0, additionalRotationField='', stackDirection=STACK_3)

cellS2=c.findAt(((3*Sehne/8, 0.0, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16)),),)
region=regionToolset.Region(cells=cellS2)
Testregion=part_ungekerbteProbe.Set(cells=cellS2, name='cellS2')

side3Face=f.getClosest(coordinates=(((3*Sehne/8, 0.0, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16))),))
side3Facev=f.findAt((((side3Face[0][1])),))
normalAxisRegion=part_ungekerbteProbe.Surface(side1Faces=side3Facev, name='Mantel_Außen_2')
primaryAxisRegion_2=e.getClosest(coordinates=(((3*Sehne/8, 0.7+D_klein/2+((D_gross/2-D_klein/2)/(Sehne/2)*3*Sehne/16),0.0)),))
primaryAxisRegionv=e.findAt((((primaryAxisRegion_2[0][1])),))
primaryAxisRegion=part_ungekerbteProbe.Set(edges=primaryAxisRegionv, name='axialeKante_Außen_2')

part_ungekerbteProbe.MaterialOrientation(region=region, orientationType=DISCRETE, axis=AXIS_1,
    normalAxisDefinition=SURFACE, normalAxisRegion=normalAxisRegion, flipNormalDirection=False, 
    normalAxisDirection=AXIS_1, primaryAxisDefinition=EDGE, primaryAxisRegion=primaryAxisRegion,
    primaryAxisDirection=AXIS_2, flipPrimaryDirection=False, additionalRotationType=ROTATION_NONE, 
    angle=0.0, additionalRotationField='', stackDirection=STACK_3)


#### Netzerstellung ####
## Seeding der einzelnen Kanten
# Axiale Kante in der Mitte
p = mdb.models['Model-1'].parts['ungekerbteProbe']
e = p.edges
EdgeMa1 =e.findAt(((0.0, (D_klein/2), 0.0),))
p.seedEdgeByNumber(edges=EdgeMa1, number=SA_a, constraint=FINER)


# Kanten an der ersten Stirnfläche
EdgeS1r1=e.findAt(((Sehne/2+h_S_Si, 0.0 , -D_gross/2+x),),)
EdgeS1r2=e.findAt(((Sehne/2+h_S_Si, -D_gross/2+x , 0.0),),)
EdgeS1r3=e.findAt(((Sehne/2+h_S_Si, 0.0 , D_gross/2-x),),)
EdgeS1r4=e.findAt(((Sehne/2+h_S_Si, (D_gross/2-x)*cos(-35*pi/180),(D_gross/2-x)*sin(-35*pi/180)),),) 
EdgeS1r5=e.findAt(((Sehne/2+h_S_Si, (D_gross/2-x)*cos(35*pi/180),(D_gross/2-x)*sin(35*pi/180)),),) 
Edge=(EdgeS1r1[0], EdgeS1r2[0], EdgeS1r3[0])
Edgek=(EdgeS1r4[0], EdgeS1r5[0])
p.seedEdgeByNumber(edges=Edge, number =SA_r/4, constraint=FIXED)
p.seedEdgeByNumber(edges=Edgek, number=SA_r/8, constraint=FIXED)


EdgeS1t1=e.findAt(((Sehne/2+h_S_Si, 0.0 , -D_klein/4),),)
EdgeS1t2=e.findAt(((Sehne/2+h_S_Si, -D_klein/4, 0.0),),)
EdgeS1t3=e.findAt(((Sehne/2+h_S_Si, 0.0, D_klein/4),),)
EdgeS1t4=e.findAt(((Sehne/2+h_S_Si, D_klein/4, 0.0),),)
EdgeS1t5=e.findAt(((Sehne/2+(h_S_Si/2), (D_gross/2-(x/2)), 0.0),),)
Edge=(EdgeS1t1[0], EdgeS1t2[0], EdgeS1t3[0], EdgeS1t4[0])
p.seedEdgeByNumber(edges=Edge, number =SA_t, constraint=FIXED)
p.seedEdgeByNumber(edges=EdgeS1t5, number=1, constraint=FIXED)


EdgeS1q1=e.findAt(((Sehne/2+h_S_Si, (D_gross/2-2*x)*sin(45*pi/180), -(D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS1q2=e.findAt(((Sehne/2+h_S_Si, -(D_gross/2-2*x)*sin(45*pi/180), -(D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS1q3=e.findAt(((Sehne/2+h_S_Si, -(D_gross/2-2*x)*sin(45*pi/180), (D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS1q4=e.findAt(((Sehne/2+h_S_Si, (D_gross/2-2*x)*sin(45*pi/180), (D_gross/2-2*x)*sin(45*pi/180)),),)
Edge=(EdgeS1q1[0],EdgeS1q2[0],EdgeS1q3[0],EdgeS1q4[0], )
#p.seedEdgeByNumber(edges=Edge, number =SA_t/2, constraint=FIXED)
p.seedEdgeByBias(biasMethod=SINGLE, end2Edges=(Edge), ratio=5.0, number=SA_t/2, constraint=FIXED)

# Kanten an der zweiten Stirnfläche
EdgeS2r1=e.findAt(((-Sehne/2-h_S_Si, 0.0 , -D_gross/2+x),),)
EdgeS2r2=e.findAt(((-Sehne/2-h_S_Si, -D_gross/2+x , 0.0),),)
EdgeS2r3=e.findAt(((-Sehne/2-h_S_Si, 0.0 , D_gross/2-x),),)
EdgeS2r4=e.findAt(((-Sehne/2-h_S_Si, (D_gross/2-x)*cos(-35*pi/180),(D_gross/2-x)*sin(-35*pi/180)),),) 
EdgeS2r5=e.findAt(((-Sehne/2-h_S_Si, (D_gross/2-x)*cos(35*pi/180),(D_gross/2-x)*sin(35*pi/180)),),) 
Edge=(EdgeS2r1[0], EdgeS2r2[0], EdgeS2r3[0])
Edgek=(EdgeS2r4[0], EdgeS2r5[0])
p.seedEdgeByNumber(edges=Edge, number =SA_r/4, constraint=FIXED)
p.seedEdgeByNumber(edges=Edgek, number=SA_r/8, constraint=FIXED)


EdgeS2t1=e.findAt(((-Sehne/2-h_S_Si, 0.0 , -D_klein/4),),)
EdgeS2t2=e.findAt(((-Sehne/2-h_S_Si, -D_klein/4, 0.0),),)
EdgeS2t3=e.findAt(((-Sehne/2-h_S_Si, 0.0, D_klein/4),),)
EdgeS2t4=e.findAt(((-Sehne/2-h_S_Si, D_klein/4, 0.0),),)
EdgeS2t5=e.findAt(((-Sehne/2-(h_S_Si/2), (D_gross/2-(x/2)), 0.0),),)
Edge=(EdgeS2t1[0], EdgeS2t2[0], EdgeS2t3[0], EdgeS2t4[0])
p.seedEdgeByNumber(edges=Edge, number =SA_t, constraint=FIXED)
p.seedEdgeByNumber(edges=EdgeS2t5, number=1, constraint=FIXED)


EdgeS2q1=e.findAt(((-Sehne/2-h_S_Si, (D_gross/2-2*x)*sin(45*pi/180), -(D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS2q2=e.findAt(((-Sehne/2-h_S_Si, -(D_gross/2-2*x)*sin(45*pi/180), -(D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS2q3=e.findAt(((-Sehne/2-h_S_Si, -(D_gross/2-2*x)*sin(45*pi/180), (D_gross/2-2*x)*sin(45*pi/180)),),)
EdgeS2q4=e.findAt(((-Sehne/2-h_S_Si, (D_gross/2-2*x)*sin(45*pi/180), (D_gross/2-2*x)*sin(45*pi/180)),),)
Edge=(EdgeS2q1[0],EdgeS2q2[0],EdgeS2q3[0],EdgeS2q4[0], )
#p.seedEdgeByNumber(edges=Edge, number =SA_t/2, constraint=FIXED)
p.seedEdgeByBias(biasMethod=SINGLE, end2Edges=(Edge), ratio=5.0, number=SA_t/2, constraint=FIXED)

# Axiale Kanten im mittleren Segment
#Edgema2=e.findAt(((-Sehne/4-1, D_klein/2+((Sehne/4+1)*((D_gross/2-D_klein/2)/Sehne/2)), 0.0),),)
#Edgema2=e.findAt(((-Sehne/4, D_klein/2 + (((D_gross/2-D_klein/2)/(Sehne/2))*(Sehne/4)), 0.0),),)
#Edgema2=e.findAt(((-Sehne/2, D_gross/2, 0.0),),)
#Edge=(Edgema2[0], )
#p.seedEdgeByNumber(edges=Edge, number=2, constraint=FIXED)

# radiale Kanten im mittleren Segment
EdgeMr1=e.getClosest(coordinates=(((-Sehne/4, 0.0, D_gross/2)),))
EdgeMr2=e.getClosest(coordinates=(((-Sehne/4, 0.0, D_gross/2)),))
EdgeMr1v=e.findAt((((EdgeMr1[0][1])),))
EdgeMr2v=e.findAt((((EdgeMr2[0][1])),))
Edge=(EdgeMr1v[0], EdgeMr2v[0], )
p.seedEdgeByNumber(edges=Edge, number=80, constraint=FIXED)

## Vernetzung
elemType1 = mesh.ElemType(elemCode=C3D20, elemLibrary=STANDARD)
elemType2 = mesh.ElemType(elemCode=C3D15, elemLibrary=STANDARD)
elemType3 = mesh.ElemType(elemCode=C3D10, elemLibrary=STANDARD)
part_ungekerbteProbe.setElementType(regions=region, elemTypes=(elemType1, elemType2, elemType3))
part_ungekerbteProbe.setMeshControls(regions=part_ungekerbteProbe.cells, algorithm=MEDIAL_AXIS)
part_ungekerbteProbe.generateMesh()

### Erstellen der Randbedingungen ###

## Assembly
Assembly1 = mdb.models['Model-1'].rootAssembly
a1 = mdb.models['Model-1'].rootAssembly
a1.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['ungekerbteProbe']
a1.Instance(name='ungekerbteProbe-1', part=p, dependent=ON)
Assembly1.ReferencePoint(point=((-Sehne/2-h_S_Si)-10, 0.0, 0.0))
Assembly1.ReferencePoint(point=((Sehne/2+h_S_Si)+10, 0.0, 0.0))

Assembly1.ReferencePoint(point=(-141, 0.0, 0.0))
Assembly1.ReferencePoint(point=(710, 0.0, 0.0))

## Constraint
r1 = Assembly1.referencePoints.findAt(((-Sehne/2-h_S_Si)-10, 0.0, 0.0),)
region1=Assembly1.Set(referencePoints=(r1, ), name='Festlager_Moment')

s1 = Assembly1.instances['ungekerbteProbe-1'].faces.getByBoundingBox((-Sehne/2-h_S_Si)-1, (-D_gross/2)-1, (-D_gross/2)-1, (-Sehne/2)+1, (D_gross/2)+1, (D_gross/2)+1)		#getByBoundingCylinder ((Punkt1 Achse), (Punkt2 Achse), radius))	
region2=Assembly1.Surface(side1Faces=s1, name='Festlager_Moment_Stirnflaeche')

mdb.models['Model-1'].Coupling(name='Festlager_Moment', controlPoint=region1, 
        surface=region2, influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, 
        localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)

r2 = Assembly1.referencePoints.findAt(((Sehne/2+h_S_Si)+10, 0.0, 0.0),)
region3=Assembly1.Set(referencePoints=(r2, ), name='Loslager_Moment')

s2 = Assembly1.instances['ungekerbteProbe-1'].faces.getByBoundingBox((Sehne/2)-1, (-D_gross/2)-1, (-D_gross/2)-1, (Sehne/2+h_S_Si)+1, (D_gross/2)+1, (D_gross/2)+1)		#getByBoundingCylinder ((Punkt1 Achse), (Punkt2 Achse), radius))	
region4=Assembly1.Surface(side1Faces=s2, name='Loslager_Moment_Stirnflaeche')

mdb.models['Model-1'].Coupling(name='Loslager_Moment', controlPoint=region3, 
        surface=region4, influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, 
        localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
		
#r3 = Assembly1.referencePoints.findAt((-141, 0.0, 0.0),)
#region5=Assembly1.Set(referencePoints=(r3, ), name='Festlager_Kraft')

#s3 = Assembly1.instances['ungekerbteProbe-1'].faces.getByBoundingBox((-Sehne/2-h_S_Si)-1, (-D_gross/2)-1, (-D_gross/2)-1, (-Sehne/2)+1, (D_gross/2)+1, (D_gross/2)+1)		#getByBoundingCylinder ((Punkt1 Achse), (Punkt2 Achse), radius))	
#region6=Assembly1.Surface(side1Faces=s3, name='Festlager_Kraft_Stirnflaeche')

#mdb.models['Model-1'].Coupling(name='Festlager_Kraft', controlPoint=region5, 
#        surface=region6, influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, 
#        localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)
		
#r4 = Assembly1.referencePoints.findAt((710, 0.0, 0.0),)
#region7=Assembly1.Set(referencePoints=(r4, ), name='Loslager_Kraft')

#s4 = Assembly1.instances['ungekerbteProbe-1'].faces.getByBoundingBox((Sehne/2)-1, (-D_gross/2)-1, (-D_gross/2)-1, (Sehne/2+h_S_Si)+1, (D_gross/2)+1, (D_gross/2)+1)		#getByBoundingCylinder ((Punkt1 Achse), (Punkt2 Achse), radius))	
#region8=Assembly1.Surface(side1Faces=s4, name='Loslager_Kraft_Stirnflaeche')

#mdb.models['Model-1'].Coupling(name='Loslager_Kraft', controlPoint=region7, 
#        surface=region8, influenceRadius=WHOLE_SURFACE, couplingType=KINEMATIC, 
#        localCsys=None, u1=ON, u2=ON, u3=ON, ur1=ON, ur2=ON, ur3=ON)


## Boundary Conditions
mdb.models['Model-1'].EncastreBC(name='Festlager_Moment', createStepName='Initial', 
									region=region1, localCsys=None) 

## Lastaufbringung
mdb.models['Model-1'].StaticStep(name='Last_Moment', previous='Initial')
Assembly1 = mdb.models['Model-1'].rootAssembly
mdb.models['Model-1'].Moment(name='Biegemoment_z', createStepName='Last_Moment', 
        region=region3, cm3=1000000.0, distributionType=UNIFORM, field='', 
        localCsys=None)

#mdb.models['Model-1'].StaticStep(name='Last_wechsel', previous='Last_Moment')
#Assembly1 = mdb.models['Model-1'].rootAssembly
#mdb.models['Model-1'].EncastreBC(name='Festlager_Kraft', createStepName='Last_wechsel',
#									region=region5, localCsys=None)
#mdb.models['Model-1'].boundaryConditions['Festlager_Moment'].deactivate(
#        'Last_wechsel')

#mdb.models['Model-1'].StaticStep(name='Last_Kraft', previous='Last_wechsel')
#Assembly1 = mdb.models['Model-1'].rootAssembly
#mdb.models['Model-1'].ConcentratedForce(name='Kraft_y', createStepName='Last_Kraft', 
#        region=region7, cf2=(1000000/(710+141)), distributionType=UNIFORM, field='', 
#        localCsys=None)
#mdb.models['Model-1'].loads['Biegemoment_z'].deactivate('Last_Kraft')


### Sets für die Auswertung anlegen
Nodes = part_ungekerbteProbe.nodes
Liste1=[]
zKriterium=0




### Job erstellen ###
mdb.Job(name='Biegung_Moment', model='Model-1', 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='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, 
        numGPUs=0)
Job1=mdb.jobs['Biegung_Moment'].submit(consistencyChecking=OFF)
mdb.jobs['Biegung_Moment'].waitForCompletion()