#Create Asymmetric geometry without notch
# *******************************************
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
from abaqusConstants import *
# ******************************************
Mdb()
session.viewports['Viewport: 1'].setValues(displayedObject=None)
# *******************************************
# a) sketch the symmetric design using line tool and constraints actautors
La=2+40
a0_up=4
a0_down=4
W=50
Wa0=10
Wa_sample=8
Ls=40
Ls_extra=2
Ws=6
x=0.75
y=1.2
n=3
x_overlap=0
x_overlap1=0
x_extra=2
x_extra1=2
tape=1#la valeur qu'on ajoute pour obtenir la largeur d'actauteur cote exterieur
v_etch=0.6#doesn't have importance in this file
Wa=Wa0+2*tape
t_etch=int(round((Wa)/(2*v_etch)))
t_over=3
if a0_up<=1:
	r_up=a0_up*float(1)/float(10)
if 1<a0_up<=14:
	r_up=a0_up*float(1)/float(100)
if 14<a0_up<25:
	r_up=a0_up*float(1)/float(30)
if 25<=a0_up<30:
	r_up=a0_up*float(1)/float(40)
if 30<=a0_up:
	r_up=a0_up*float(1)/float(50)
if a0_down<=1:
	r_down=a0_down*float(1)/float(10)
if 1<a0_down<=14:
	r_down=a0_down*float(1)/float(100)
if 14<a0_down<25:
	r_down=a0_down*float(1)/float(30)
if 25<=a0_down<30:
	r_down=a0_down*float(1)/float(40)
if 30<=a0_down:
	r_down=a0_down*float(1)/float(50)
B=0.07
Ea=207*pow(10,9)
Es=100*pow(10,9)
alpha_a=4e-06
alpha_s=0
Delta=-507.25
mu_s=0.3
mu_a=0.3
# c1=a_1+n
Wa_extra = 20
La_extra = 20
# W_free=40#should be higher than tape since already contains tape value
meshsize1=0.5
meshsize2=1
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
sheetSize=200.0)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.Line(point1=(0.0, 0.0), point2=(-La, -tape))
s1.Line(point1=(-La, -tape), point2=(-La, Wa0+tape))
s1.Line(point1=(-La, Wa0+tape), point2=(0.0, Wa0))
s1.Line(point1=(0.0, Wa0), point2=(0.0, 0.0))
s1.Line(point1=(Ws+x_extra+x_extra1 ,Wa0), point2=(Ws+x_extra+x_extra1+La, Wa0+tape))
s1.Line(point1=(Ws+x_extra+x_extra1+La, Wa0+tape), point2=(x_extra+x_extra1+Ws+La, -tape))
s1.Line(point1=(x_extra+x_extra1+Ws+La, -tape), point2=(x_extra+x_extra1+Ws, 0))
s1.Line(point1=(x_extra+x_extra1+Ws, 0), point2=(x_extra+x_extra1+Ws, Wa0))

# b) Create a  part named "Part-1"
p = mdb.models['Model-1'].Part(name='Part-1', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Part-1']
p.BaseShell(sketch=s1)
s1.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['Part-1']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
# Sketch Specimen part
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
sheetSize=200.0)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.Line(point1=(-x_overlap, float(Wa0-Wa_sample)/2), point2=(-x_overlap, Wa_sample+float(Wa0-Wa_sample)/2))
s1.Line(point1=(-x_overlap, Wa_sample+float(Wa0-Wa_sample)/2), point2=(+x_extra, Wa_sample+float(Wa0-Wa_sample)/2))
s1.Line(point1=(x_extra, Wa_sample+float(Wa0-Wa_sample)/2), point2=(x_extra,Wa_sample+float(Wa0-Wa_sample)/2+Ls+Ls_extra))
s1.Line(point1=(x_extra, Wa_sample + float(Wa0 - Wa_sample) / 2+Ls+Ls_extra),point2=(x_extra+Ws, Wa_sample+float(Wa0-Wa_sample)/2+Ls+Ls_extra))
s1.Line(point1=(x_extra+Ws, Wa_sample+float(Wa0-Wa_sample)/2+Ls+Ls_extra), point2=(x_extra+Ws,float(Wa0-Wa_sample)/2+Wa_sample))
s1.Line(point1=(x_extra+Ws,float(Wa0-Wa_sample)/2+Wa_sample), point2=(x_extra+Ws+x_extra1+x_overlap1, float(Wa0-Wa_sample)/2+Wa_sample))
s1.Line(point1=(x_extra+Ws+x_extra1+x_overlap1, float(Wa0-Wa_sample)/2+Wa_sample), point2=(x_extra+Ws+x_extra1+x_overlap1, float(Wa0-Wa_sample)/2))
s1.Line(point1=(x_extra+x_overlap1+x_extra1+Ws, float(Wa0-Wa_sample)/2), point2=(x_extra+Ws, float(Wa0-Wa_sample)/2))
s1.Line(point1=(x_extra+Ws,float(Wa0-Wa_sample)/2), point2=(Ws+x_extra, -Ls+float(Wa0-Wa_sample)/2-Ls_extra))
s1.Line(point1=(Ws+x_extra, -Ls+float(Wa0-Wa_sample)/2-Ls_extra), point2=(x_extra, -Ls+float(Wa0-Wa_sample)/2-Ls_extra))
s1.Line(point1=(x_extra,-Ls+float(Wa0-Wa_sample)/2-Ls_extra), point2=(x_extra, float(Wa0-Wa_sample)/2))
s1.Line(point1=(x_extra, float(Wa0-Wa_sample)/2), point2=(-x_overlap, float(Wa0-Wa_sample)/2))
#Define notch
s1.Line(point1=(float(Ws)/2+x_extra, float(Wa0-Wa_sample)/2+float(Wa_sample)/2-n), point2=(x_extra+float(Ws)/2+x, float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n+y))
s1.Line(point1=(x_extra+float(Ws)/2+x, float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n+y), point2=(x_extra+x+float(Ws)/2, float(Wa0-Wa_sample)/2+float(Wa_sample)/2+n-y))
s1.Line(point1=(x_extra+x+float(Ws)/2, float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n-y), point2=(x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+float(Wa_sample)/2+n))
s1.Line(point1=(x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+n+float(Wa_sample)/2), point2=(-x+x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+float(Wa_sample)/2+n-y))
s1.Line(point1=(-x+x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+n-y+float(Wa_sample)/2), point2=(-x+x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+float(Wa_sample)/2-n+y))
s1.Line(point1=(x_extra+float(Ws)/2-x, float(Wa0-Wa_sample)/2-n+y+float(Wa_sample)/2), point2=(x_extra+float(Ws)/2, float(Wa0-Wa_sample)/2+float(Wa_sample)/2-n))

# b) Create a  part named "Part-3"
p = mdb.models['Model-1'].Part(name='Part-3', dimensionality=TWO_D_PLANAR, type=DEFORMABLE_BODY)
p = mdb.models['Model-1'].parts['Part-3']
p.BaseShell(sketch=s1)
s1.unsetPrimaryObject()
p = mdb.models['Model-1'].parts['Part-3']
session.viewports['Viewport: 1'].setValues(displayedObject=p)
del mdb.models['Model-1'].sketches['__profile__']
# c) sketch contour
s = mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=400.0)
g, v, d, c = s.geometry, s.vertices, s.dimensions, s.constraints
s.setPrimaryObject(option=STANDALONE)
s.rectangle(point1=(-La-La_extra, float(Wa0-Wa_sample)/2-Ls-Wa_extra), point2=(x_extra1+x_extra+Ws+La+La_extra,Wa_sample+float(Wa0-Wa_sample)/2+Ls+Wa_extra))
s.Line(point1=(-La, float(Wa0-Wa_sample)/2-Ls), point2=(x_extra+x_extra1+Ws+La,float(Wa0-Wa_sample)/2-Ls))
s.Line(point1=(x_extra+x_extra1+Ws+La,float(Wa0-Wa_sample)/2-Ls), point2=(x_extra+x_extra1+Ws+La,Wa_sample+float(Wa0-Wa_sample)/2+Ls))
s.Line(point1=(x_extra+x_extra1+Ws+La,Wa_sample+float(Wa0-Wa_sample)/2+Ls), point2=(-La,Wa_sample+float(Wa0-Wa_sample)/2+Ls))
s.Line(point1=(-La,Wa_sample+float(Wa0-Wa_sample)/2+Ls), point2=(-La,float(Wa0-Wa_sample)/2-Ls))

#b) Create a  part named "Part-2"
p1 = mdb.models['Model-1'].Part(name='Part-2', dimensionality=TWO_D_PLANAR,
								type=DEFORMABLE_BODY)
p1 = mdb.models['Model-1'].parts['Part-2']
p1.BaseShell(sketch=s)
s.unsetPrimaryObject()
session.viewports['Viewport: 1'].setValues(displayedObject=p1)
del mdb.models['Model-1'].sketches['__profile__']
# Create a  partition for crack contour calculus and partition between sample and actuator
p1 = mdb.models['Model-1'].parts['Part-3']
f1, e1, d2 = p1.faces, p1.edges, p1.datums
t = p1.MakeSketchTransform(sketchPlane=f1.findAt(coordinates=(Ws - 0.002, Wa_sample - 0.002,
															  0.0), normal=(0.0, 0.0, 1.0)), sketchPlaneSide=SIDE1,
						   origin=(0, Wa_sample+float(Wa0-Wa_sample)/2, 0.0))
s1 = mdb.models['Model-1'].ConstrainedSketch(name='__profile__',
											 sheetSize=188.5, gridSpacing=4.71, transform=t)
g, v, d, c = s1.geometry, s1.vertices, s1.dimensions, s1.constraints
s1.setPrimaryObject(option=SUPERIMPOSE)
p1 = mdb.models['Model-1'].parts['Part-3']
p1.projectReferencesOntoSketch(sketch=s1, filter=COPLANAR_EDGES)
s1.Line(point1=(float(Ws)/2+x_extra, -float(Wa_sample)/2-n), point2=(float(Ws)/2+x_extra, -float(Wa_sample)/2 - a0_down-n))
s1.CircleByCenterPerimeter(center=(float(Ws)/2+x_extra,-float(Wa_sample)/2 -n- a0_down), point1=(-r_down +float(Ws)/2+x_extra, -float(Wa_sample)/2-n- a0_down))
s1.Line(point1=(float(Ws)/2+x_extra, -float(Wa_sample)/2+n), point2=(float(Ws)/2+x_extra, -float(Wa_sample)/2+n+a0_up))
s1.CircleByCenterPerimeter(center=(float(Ws)/2+x_extra,-float(Wa_sample)/2 +n+a0_up), point1=(-r_up +float(Ws)/2+x_extra, -float(Wa_sample)/2+n+a0_up))

p1 = mdb.models['Model-1'].parts['Part-3']
f1 = p1.faces
pickedFaces = f1.findAt(((Ws+x_extra - 0.002, Wa_sample - 0.002, 0.0),))
e1, d2 = p1.edges, p1.datums
p1.PartitionFaceBySketch(faces=pickedFaces, sketch=s1)
s1.unsetPrimaryObject()
del mdb.models['Model-1'].sketches['__profile__']
#----------------------------------------------
#Create material SiN
mdb.models['Model-1'].Material(name='SiN')
mdb.models['Model-1'].materials['SiN'].Elastic(table=((Ea, mu_a), ))
mdb.models['Model-1'].materials['SiN'].Expansion(table=((alpha_a, ), ))
#Create material of Specimen
mdb.models['Model-1'].Material(name='Specimen')
mdb.models['Model-1'].materials['Specimen'].Elastic(table=((Es, mu_s), ))
mdb.models['Model-1'].materials['Specimen'].Expansion(table=((alpha_s, ), ))
# Create a section
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-1', material='SiN',
thickness=B)
mdb.models['Model-1'].HomogeneousSolidSection(name='Section-2', material='Specimen',
thickness=B)
# Assign a section
p = mdb.models['Model-1'].parts['Part-1']
f = p.faces
faces = f.findAt(((-La/2, Wa_sample/2, 0.0), ),((Ws+x_extra+x_extra1+x_overlap1+La/2,Wa_sample/2,0.0 ), ))
region = p.Set(faces=faces, name='Actuators')#Create Set
p.SectionAssignment(region=region, sectionName='Section-1', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='', thicknessAssignment=FROM_SECTION)
p1 = mdb.models['Model-1'].parts['Part-2']
f1 = p1.faces
faces1=f1.findAt(((-La-La_extra/2, 0.0, 0.0), ))
faces2=f1.findAt(((Ws+La+La_extra/2+x_extra+x_extra1+x_overlap1, 0.0, 0.0), ))
faces3=f1.findAt(((x_extra+Ws/2, float(Wa0-Wa_sample)/2-Ls-Wa_extra/2, 0.0), ))
faces4=f1.findAt(((x_extra+float(Ws)/2,Wa0+Ls+Wa_extra/2-float(Wa0-Wa_sample)/2,0.0 ), ))
p1.Set(faces=faces1+faces2+faces3+faces4, name='UnderetchingZone')
region1 = p1.Set(faces=faces1, name='UnderetchingZone')
p1.SectionAssignment(region=region1, sectionName='Section-1', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
p1 = mdb.models['Model-1'].parts['Part-3']
f1 = p1.faces
faces1=f1.findAt(((x_extra, Wa_sample/2, 0.0), ))
faces2=f1.findAt(((x_extra+float(Ws)/2,float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n-a0_down-r_down/4,0.0), ))
faces3=f1.findAt(((x_extra+float(Ws)/2,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n+a0_up+r_up/4,0.0), ))
p1.Set(faces=faces1+faces2+faces3, name='Specimen')
region1 = p1.Set(faces=faces1+faces2+faces3, name='Specimen')
p1.SectionAssignment(region=region1, sectionName='Section-2', offset=0.0,
offsetType=MIDDLE_SURFACE, offsetField='',
thicknessAssignment=FROM_SECTION)
#Assembly and create Instance
a = mdb.models['Model-1'].rootAssembly
a.DatumCsysByDefault(CARTESIAN)
p = mdb.models['Model-1'].parts['Part-1']
a.Instance(name='Part-1-1', part=p, dependent=ON)
p = mdb.models['Model-1'].parts['Part-2']
a.Instance(name='Part-2-1', part=p, dependent=ON)
p = mdb.models['Model-1'].parts['Part-3']
a.Instance(name='Part-3-1', part=p, dependent=ON)
#Create merged part from two others parts
a = mdb.models['Model-1'].rootAssembly
a.InstanceFromBooleanMerge(name='Part-4', instances=(a.instances['Part-1-1'], a.instances['Part-2-1'],a.instances['Part-3-1'], ), keepIntersections=ON,
                           originalInstances=SUPPRESS, domain=GEOMETRY)
#Create Set of specimen
a = mdb.models['Model-1'].rootAssembly
e2 = a.instances['Part-2-1'].edges
e1 = a.instances['Part-4-1'].edges
edges1 = e1.findAt(((-float(La)/2, -Ls+float(Wa0-Wa_sample)/2, 0.0), ))
edges2=e1.findAt(((x_extra,-Ls/2+float(Wa0-Wa_sample)/2,0.0), ))
edges3=e1.findAt(((-float(La)/2,-float(tape)/2,0.0), ))
edges4=e1.findAt(((-La,-float(tape)-0.1,0.0), ))
# edges5=e1.findAt(((-La,0.0,0.0), ))
edges6=e1.findAt(((-La,Wa0+float(Ls)/2,0.0), ))
edges7=e1.findAt(((-float(La)/2,Wa_sample+Ls+float(Wa0-Wa_sample)/2,0.0), ))
# edges8=e1.findAt(((float(Ws)/2,Wa0+W_free,0.0), ))
edges9=e1.findAt(((Ws+float(La)/2,Wa_sample+Ls+float(Wa0-Wa_sample)/2,0.0), ))
edges10=e1.findAt(((Ws+La+x_extra+x_extra1,Wa_sample+float(Wa0-Wa_sample)/2+Ls/2,0.0), ))
edges11=e1.findAt(((Ws+La+x_extra+x_extra1,float(Wa0)/2,0.0), ))
edges12=e1.findAt(((Ws+La+x_extra1+x_extra,-float(Ls)/2+float(Wa0-Wa_sample)/2,0.0), ))
edges13=e1.findAt(((Ws+x_extra+x_extra1+float(La)/2,-Ls+float(Wa0-Wa_sample)/2,0.0), ))
# edges14=e1.findAt(((float(Ws)/2,-Ls,0.0), ))
edges15=e1.findAt(((-float(La)/2,Wa0+float(tape)/2,0.0), ))
edges16=e1.findAt(((Ws+float(La)/2+x_extra1+x_extra,Wa0+float(tape)/2,0.0), ))
edges17=e1.findAt(((Ws+float(La)/2+x_extra1+x_extra,-float(tape)/2,0.0), ))
edges18=e1.findAt(((Ws+x_extra,-float(Ls)/2+float(Wa0-Wa_sample)/2,0.0), ))
edges19=e1.findAt(((0.0+x_extra/4,float(Wa0-Wa_sample)/2,0.0), ))
edges20=e1.findAt(((x_extra,Wa_sample+float(Wa0-Wa_sample)/2+Ls/2,0.0), ))
edges21=e1.findAt(((float(Ws)/2+x+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n/100,0.0), ))
edges22=e1.findAt(((x_extra+Ws,Wa_sample+float(Wa0-Wa_sample)/2+Ls/2,0.0), ))
edges23=e1.findAt(((float(Ws)/2-float(x)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n-float(y)/2,0.0), ))
edges41=e1.findAt(((float(Ws)/2-float(x)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n+float(y)/2,0.0), ))
edges24=e1.findAt(((float(Ws)/2+float(x)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n+float(y)/2,0.0), ))
edges25=e1.findAt(((float(Ws)/2+float(x)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n-float(y)/2,0.0), ))
edges26=e1.findAt(((float(Ws)/2-x+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n/100,0.0), ))
edges27=e1.findAt(((Ws+x_extra+x_extra1/2,Wa_sample+float(Wa0-Wa_sample)/2,0.0), ))
edges28=e1.findAt(((Ws+x_extra+x_extra1/2,float(Wa0-Wa_sample)/2,0.0), ))
edges29=e1.findAt(((x_extra/2,float(Wa0-Wa_sample)/2,0.0), ))
edges40=e1.findAt(((x_extra/2,float(Wa0-Wa_sample)/2+Wa_sample,0.0), ))
edges30=e1.findAt(((0.0,float(Wa0-Wa_sample)/4,0.0), ))
edges31=e1.findAt(((0.0,Wa0-float(Wa0-Wa_sample)/4,0.0), ))
edges32=e1.findAt(((Ws+x_extra+x_extra1,Wa0-float(Wa0-Wa_sample)/4,0.0), ))
edges33=e1.findAt(((Ws+x_extra+x_extra1,float(Wa0-Wa_sample)/4,0.0), ))
edges34=e1.findAt(((x_extra/2,-Ls+float(Wa0-Wa_sample)/2,0.0), ))
edges35=e1.findAt(((Ws+x_extra+x_extra1/2,-Ls+float(Wa0-Wa_sample)/2,0.0), ))
edges36=e1.findAt(((Ws+x_extra+x_extra1+La/2,Ls+float(Wa0-Wa_sample)/2+Wa_sample,0.0), ))
edges37=e1.findAt(((-La/2,Ls+float(Wa0-Wa_sample)/2+Wa_sample,0.0), ))
edges38=e1.findAt(((Ws+x_extra+x_extra1+La/2,-Ls+float(Wa0-Wa_sample)/2,0.0), ))
edges39=e1.findAt(((-La/2,-Ls+float(Wa0-Wa_sample)/2,0.0), ))
# # edges26=e1.findAt(((float(Ws)/2,Wa0-n-float(a_1)/2,0.0), ))
a.Set(edges=edges1+edges2+edges3+edges4+edges6+edges7+edges9+edges10+edges12+edges13+
            edges15+edges16+edges17+edges18+edges19+edges20+edges21+edges22+edges23+edges24+edges25+edges26+edges27+
			edges28+edges29+edges30+edges31+edges32+edges33+edges34+edges35+edges36+edges37+edges38+edges39+edges40+edges41,name='Contour')

a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
faces = f1.findAt(((-La+0.002, Wa0-0.002, 0.0), ), ((x_extra+Ws-0.002, -Ls+float(Wa0-Wa_sample)/2-Ls_extra/2, 0.0), ),
((x_extra+Ws+0.002, Wa_sample-0.002, 0.0), ),((float(Ws)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n+a0_up,0.0 ), ),
				  ((float(Ws)/2+x_extra,float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n-a0_down,0.0),),
				   ((Ws+x_extra+x_extra1+La/2,float(Wa_sample)/2+float(Wa0-Wa_sample)/2,0.0), ))
f2 = a.instances['Part-4-1'].faces
faces2=f2.findAt(((-La-La_extra/2, 0.0, 0.0), ), ((x_extra+Ws+x_extra1+La+La_extra/2, 0.0, 0.0), ),((x_extra+Ws/2, Ls+float(Wa0-Wa_sample)/2+Ls_extra/2+Wa_sample, 0.0), ),
((x_extra+Ws/2, -Ls+float(Wa0-Wa_sample)/2-Wa_extra/2, 0.0), ),((x_extra+float(Ws)/2,Wa_sample+float(Wa0-Wa_sample)/2+Ls+Wa_extra/2,0.0 ), ))
a.Set(faces=faces+faces2, name='Whole')
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
faces1 = f1.findAt(((-La - La_extra / 2, 0.0, 0.0),))
faces2 = f1.findAt(((Ws + La + La_extra / 2 + x_extra + x_extra1 , 0.0, 0.0),))
faces3 = f1.findAt(((x_extra + Ws / 2, float(Wa0 - Wa_sample) / 2 - Ls - Wa_extra / 2, 0.0),))
faces4 = f1.findAt(((x_extra + float(Ws) / 2, Wa0 + Ls + Wa_extra / 2 - float(Wa0 - Wa_sample) / 2, 0.0),))
faces5 = f1.findAt(((Ws/ 2+x_extra, -Ls+float(Wa0 - Wa_sample) / 2-Ls_extra/2, 0.0),))
faces6 = f1.findAt(((Ws/ 2+x_extra, Ls+float(Wa0 - Wa_sample) / 2+Wa_sample+Ls_extra/2, 0.0),))
a.Set(faces=faces1 + faces2 + faces3 + faces4+faces5+faces6, name='UnderetchingZone')
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].edges
edges1 = e1.findAt(((0.0-x_overlap, Wa_sample-0.002, 0.0), ))
a.Set(edges=edges1, name='LeftSpecimen')
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].edges
edges1 = e1.findAt(((x_extra+x_extra1+Ws+x_overlap1, Wa_sample-0.002, 0.0), ))
a.Set(edges=edges1, name='RightSpecimen')
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
faces1 = f1.findAt(((La+Ws+x_extra1+x_extra, Wa0/2, 0.0), ))
a.Set(faces=faces1, name='RightActuator')
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
faces1 = f1.findAt(((-La, Wa0/2, 0.0), ))
a.Set(faces=faces1, name='LeftActuator')
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
faces1=f1.findAt(((x_extra, Wa_sample/2, 0.0), ))
faces2=f1.findAt(((x_extra+float(Ws)/2,float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n-a0_down-r_down/4,0.0), ))
faces3=f1.findAt(((x_extra+float(Ws)/2,float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n+a0_up+r_up/4,0.0), ))
faces4=f1.findAt(((x_extra+float(Ws)/2,Wa_sample+float(Wa0-Wa_sample)/2+Ls+Ls_extra/2,0.0), ))
faces5=f1.findAt(((x_extra+float(Ws)/2,+float(Wa0-Wa_sample)/2-Ls-Ls_extra/2,0.0), ))
a.Set(faces=faces1+faces2+faces3+faces4+faces5, name='Specimen')#Create Set Specimen

# Create crack seam
a = mdb.models['Model-1'].rootAssembly
a.makeIndependent(instances=(a.instances['Part-4-1'], ))
e1 = a.instances['Part-4-1'].edges
if a0_down!=0:
    edges1 = e1.findAt(((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - n - 0.002, 0.0),),
                       ((float(Ws) / 2 + x_extra,
                         float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n + r_down / 5,
                         0.0),))  # one point coordinates will be sufficient
    pickedRegions1 = a.Set(edges=edges1, name='cracktip1')
    mdb.models['Model-1'].rootAssembly.engineeringFeatures.assignSeam(
        regions=pickedRegions1)

if a0_up!=0:
    edges2 = e1.findAt(((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + n + 0.002, 0.0),),
                       ((float(Ws) / 2 + x_extra,
                         float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n - r_up / 5,
                         0.0),))  # one point coordinates will be sufficient
    pickedRegions2 = a.Set(edges=edges2, name='cracktip2')
    mdb.models['Model-1'].rootAssembly.engineeringFeatures.assignSeam(
        regions=pickedRegions2)
# Create crack
if a0_up!=0:
    a = mdb.models['Model-1'].rootAssembly
    f1 = a.instances['Part-4-1'].faces
    faces1 = f1.findAt(((float(Ws) / 2 + 0.5 * r_down + x_extra,
                         float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n + r_down * 0.1, 0.0),))
    crackFront = regionToolset.Region(faces=faces1)
    a = mdb.models['Model-1'].rootAssembly
    v1 = a.instances['Part-4-1'].vertices
    verts1 = v1.findAt(
        ((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n, 0.0),))
    crackTip = regionToolset.Region(vertices=verts1)
    v11 = a.instances['Part-4-1'].vertices
    v21 = a.instances['Part-4-1'].vertices
    a.engineeringFeatures.ContourIntegral(name='Crack-1', symmetric=OFF,
                                          crackFront=crackFront, crackTip=crackTip,
                                          extensionDirectionMethod=Q_VECTORS, qVectors=((v11.findAt(coordinates=(
            float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n + r_down, 0.0)),
                                                                                         v21.findAt(coordinates=(
                                                                                         float(Ws) / 2 + x_extra,
                                                                                         float(Wa_sample) / 2 + float(
                                                                                             Wa0 - Wa_sample) / 2 - a0_down - n,
                                                                                         0.0))),), midNodePosition=0.25,
                                          collapsedElementAtTip=SINGLE_NODE)

if a0_up!=0:
    a = mdb.models['Model-1'].rootAssembly
    f1 = a.instances['Part-4-1'].faces
    faces1 = f1.findAt(((float(Ws) / 2 + 0.5 * r_up + x_extra,
                         float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n - r_up * 0.1, 0.0),))
    crackFront = regionToolset.Region(faces=faces1)
    a = mdb.models['Model-1'].rootAssembly
    v1 = a.instances['Part-4-1'].vertices
    verts1 = v1.findAt(((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n, 0.0),))
    crackTip = regionToolset.Region(vertices=verts1)
    v11 = a.instances['Part-4-1'].vertices
    v21 = a.instances['Part-4-1'].vertices
    a.engineeringFeatures.ContourIntegral(name='Crack-2', symmetric=OFF,
                                          crackFront=crackFront, crackTip=crackTip,
                                          extensionDirectionMethod=Q_VECTORS, qVectors=((v11.findAt(coordinates=(
            float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n - r_up, 0.0)),
                                                                                         v21.findAt(coordinates=(
                                                                                         float(Ws) / 2 + x_extra,
                                                                                         float(Wa_sample) / 2 + float(
                                                                                             Wa0 - Wa_sample) / 2 + a0_up + n,
                                                                                         0.0))),), midNodePosition=0.25,
                                          collapsedElementAtTip=SINGLE_NODE)
# Create steps
mdb.models['Model-1'].StaticStep(name='beforeRelease', previous='Initial',
maxNumInc=1000, initialInc=0.001, minInc=1e-10)
mdb.models['Model-1'].StaticStep(name='Release', previous='beforeRelease',
maxNumInc=1000, initialInc=0.001, minInc=1e-10)
# Create boundaryConditions
a = mdb.models['Model-1'].rootAssembly
region = a.sets['Whole']
mdb.models['Model-1'].TemperatureBC(name='initialtemperature',
createStepName='Initial', region=region, distributionType=UNIFORM,
fieldName='', magnitude=0.0)
a = mdb.models['Model-1'].rootAssembly
region = a.sets['Whole']
mdb.models['Model-1'].EncastreBC(name='Fixedstructure',
createStepName='Initial', region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
region = a.sets['UnderetchingZone']
mdb.models['Model-1'].EncastreBC(name='FixedUnderEtchingZone',
createStepName='beforeRelease', region=region, localCsys=None)

a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].faces
faces1 = e1.findAt(((Ws/2+x_extra,-Ls+float(Wa0-Wa_sample)/2-x_overlap/2,0.0), ))
region = a.Set(faces=faces1, name='down')
mdb.models['Model-1'].EncastreBC(name='downOverlap', createStepName='Initial',
region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].faces
faces1 = e1.findAt(((Ws/2+x_extra,Ls+float(Wa0-Wa_sample)/2+x_overlap/2+Wa_sample,0.0), ))
region = a.Set(faces=faces1, name='up')
mdb.models['Model-1'].EncastreBC(name='upOverlap', createStepName='Initial',
region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].edges
edges1 = e1.findAt(((Ws+x_extra+x_extra1+La, Wa0-0.002, 0.0), ))
region = a.Set(edges=edges1, name='Right')
mdb.models['Model-1'].EncastreBC(name='RightEdge', createStepName='Initial',
region=region, localCsys=None)
a = mdb.models['Model-1'].rootAssembly
e1 = a.instances['Part-4-1'].edges
edges1 = e1.findAt(((-La, Wa0-0.002, 0.0), ))
region = a.Set(edges=edges1, name='Left')
mdb.models['Model-1'].EncastreBC(name='LeftEdge', createStepName='Initial',
region=region, localCsys=None)
# Create predefinedFields
a = mdb.models['Model-1'].rootAssembly
region = a.sets['Whole']
mdb.models['Model-1'].Temperature(name='Predefined Field-1',
createStepName='beforeRelease', region=region,
distributionType=UNIFORM,
crossSectionDistribution=CONSTANT_THROUGH_THICKNESS, magnitudes=(
Delta, ))
# deactivate fixed boundaryConditions in 'release' step
mdb.models['Model-1'].boundaryConditions['Fixedstructure'].deactivate(
'Release')
# Assign Mesh type
a.makeIndependent(instances=(a.instances['Part-4-1'], ))
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
pickedRegions = f1.findAt(((Ws+x_extra1+x_extra+La/2, Wa0/2, 0.0), ), ((-La/2,
Wa0/2, 0.0), ),((-x_overlap/2,Wa_sample/2+float(Wa0-Wa_sample)/2,0.0),),
						  ((Ws+x_extra+x_extra1+x_overlap1/2,float(Wa0-Wa_sample)/2+Wa_sample/2,0.0),))#both actuators are selected

a.setMeshControls(regions=pickedRegions, technique=STRUCTURED,elemShape=QUAD)
a = mdb.models['Model-1'].rootAssembly

e1 = a.instances['Part-4-1'].faces
faces1=e1.findAt(((-La-La_extra/2, 0.0, 0.0), ))
faces2=e1.findAt(((Ws+x_extra+x_extra1+La+La_extra/2, 0.0, 0.0), ))
faces3=e1.findAt(((Ws/2+x_extra,-Ls+float(Wa0-Wa_sample)/2-Wa_extra/2, 0.0), ))
faces4=e1.findAt(((float(Ws)/2+x_extra,Wa_sample+float(Wa0-Wa_sample)/2+Ls+Wa_extra/2,0.0 ), ))
a.setMeshControls(regions=faces1+faces2+faces3+faces4, elemShape=QUAD, algorithm=MEDIAL_AXIS)
pickedRegions = f1.findAt(((Ws+x_extra-0.002, Wa_sample+float(Wa0-Wa_sample)/2-0.002, 0.0), ))
a.setMeshControls(regions=pickedRegions, elemShape=QUAD, algorithm=MEDIAL_AXIS)
a = mdb.models['Model-1'].rootAssembly
f1 = a.instances['Part-4-1'].faces
pickedRegions = f1.findAt(((float(Ws)/2+x_extra, float(Wa_sample)/2+float(Wa0-Wa_sample)/2-n-a0_down, 0.0), ))# circle
a.setMeshControls(regions=pickedRegions, technique=SWEEP)
pickedRegions1 = f1.findAt(((float(Ws)/2+x_extra, float(Wa_sample)/2+float(Wa0-Wa_sample)/2+n+a0_up, 0.0), ))# circle
a.setMeshControls(regions=pickedRegions1, technique=SWEEP)
# Seeding
if a0_up != 0:
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(((float(Ws) / 2 + x_extra + float(r_up) / 2,
							  float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + n + a0_up + sqrt(3) * float(r_up) / 2,
							  0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=12, constraint=FINER)
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(((float(Ws) / 2 + x_extra - float(r_up) / 2,
							  float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + n + a0_up + sqrt(3) * float(r_up) / 2,
							  0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FINER)
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(
		((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n - 0.001, 0.0),),
		((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 + a0_up + n, 0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FINER)
if a0_down != 0:
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(((float(Ws) / 2 + x_extra + float(r_down) / 2,
							  float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - n - a0_down + sqrt(3) * float(
								  r_down) / 2, 0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=12, constraint=FINER)
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(((float(Ws) / 2 + x_extra - float(r_down) / 2,
							  float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - n - a0_down + sqrt(3) * float(
								  r_down) / 2, 0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=4, constraint=FINER)
	a = mdb.models['Model-1'].rootAssembly
	e1 = a.instances['Part-4-1'].edges
	pickedEdges = e1.findAt(
		((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n + 0.001, 0.0),),
		((float(Ws) / 2 + x_extra, float(Wa_sample) / 2 + float(Wa0 - Wa_sample) / 2 - a0_down - n, 0.0),))
	a.seedEdgeByNumber(edges=pickedEdges, number=1, constraint=FINER)
	a = mdb.models['Model-1'].rootAssembly
	a = mdb.models['Model-1'].rootAssembly
	partInstances = (a.instances['Part-4-1'],)
	a.seedPartInstance(regions=partInstances, size=meshsize1, deviationFactor=0.1,
					   minSizeFactor=0.1)
a3 = mdb.models['Model-1'].rootAssembly
e1 = a3.instances['Part-4-1'].edges
pickedEdges = e1.findAt(((-La, tape - Ls / 2, 0.0),), ((-La, Wa0 / 2, 0.0),), ((
																				   -La, Wa_sample + float(
																					   Wa0 - Wa_sample) / 2 + Ls / 2,
																				   0.0),),
						((-La / 2, Wa_sample + float(Wa0 - Wa_sample) / 2 + Ls, 0.0),),
						((x_extra + Ws / 2, Wa_sample + Ls + float(Wa0 - Wa_sample) / 2, 0.0),), ((
																									  Ws + La / 2 + x_extra + x_extra1,
																									  Wa_sample + Ls + float(
																										  Wa0 - Wa_sample) / 2,
																									  0.0),),
						((Ws + x_extra + x_extra1 + La, Wa_sample + float(Wa0 - Wa_sample) / 2 + Ls / 2, 0.0),),
						((Ws + La + x_extra + x_extra1, Wa_sample / 2 + float(Wa0 - Wa_sample) / 2, 0.0),), ((
																												 Ws + La + x_extra + x_extra1,
																												 -tape - Ls / 2 + float(
																													 Wa0 - Wa_sample) / 2,
																												 0.0),),
						((La / 2 + Ws + x_extra1 + x_extra, float(Wa0 - Wa_sample) / 2 - Ls, 0.0),),
						((Ws / 2 + x_extra, float(Wa0 - Wa_sample) / 2 - Ls, 0.0),), (
							(-La / 2, float(Wa0 - Wa_sample) / 2 - Ls, 0.0),),
						((x_extra + Ws / 2, float(Wa0 - Wa_sample) / 2 - Ls - Wa_extra, 0.0),),
						((x_extra + x_extra1 + Ws + La + La_extra, 0.0, 0.0),), (
							(Ws / 2 + x_extra, Wa_sample + Ls + float(Wa0 - Wa_sample) / 2 + Wa_extra, 0.0),),
						((-La - La_extra, 0.0, 0.0),))
a3.seedEdgeBySize(edges=pickedEdges, size=meshsize2, deviationFactor=0.1,
				  minSizeFactor=0.1, constraint=FINER)
#generate mesh
a = mdb.models['Model-1'].rootAssembly
partInstances =(a.instances['Part-4-1'], )
a.generateMesh(regions=partInstances)
#Create HistoryOutputRequest
mdb.models['Model-1'].historyOutputRequests['H-Output-1'].setValues(
contourIntegral='Crack-1', sectionPoints=DEFAULT, rebar=EXCLUDE,
numberOfContours=5, contourType=K_FACTORS, kFactorDirection=MERR)
mdb.models['Model-1'].HistoryOutputRequest(name='H-Output-2',
createStepName='beforeRelease', contourIntegral='Crack-2',
sectionPoints=DEFAULT, rebar=EXCLUDE, numberOfContours=5,
contourType=K_FACTORS, kFactorDirection=MERR)
regionDef=mdb.models['Model-1'].rootAssembly.sets['LeftSpecimen']
mdb.models['Model-1'].HistoryOutputRequest(name='Left Specimen Output',
createStepName='beforeRelease', variables=('S11', 'S22', 'S12', 'U1', 'U2', 'RF1', 'RF2'),
region=regionDef, sectionPoints=DEFAULT, rebar=EXCLUDE)#beforerelease step because it will start calculate from this step and go on to the next steps
regionDef1=mdb.models['Model-1'].rootAssembly.sets['RightSpecimen']
mdb.models['Model-1'].HistoryOutputRequest(name='Right Specimen Output',
createStepName='beforeRelease', variables=('S11', 'S22', 'S12', 'U1', 'U2', 'RF1', 'RF2','LE'),
region=regionDef1, sectionPoints=DEFAULT, rebar=EXCLUDE)
#Create fieldOutputRequests
mdb.models['Model-1'].fieldOutputRequests['F-Output-1'].setValues(variables=(
'S', 'E', 'PE', 'PEEQ', 'PEMAG', 'LE', 'U', 'RF', 'CF', 'TF',
'CSTRESS', 'CDISP', 'IVOL', 'PHILSM', 'PSILSM'))
#Create job
# job_name='OverCNTSymEs'+repr(int(Es/pow(10,9)))+'alpha_s'+repr(int(alpha_s/pow(10,-6)))+'La'+repr(La)+'aupDown'+repr(int(a0_up*100))+'Ls'+repr(int(Ls*100))+'Ws'+repr(int(Ws*100))+ 'Wa0'+repr(int(Wa0*100))+'y'+repr(int(y*100))+'x'+repr(int(x*100))+'extraW'+repr(Wa_extra)+'L'+repr(La_extra)+'Xextra'+repr(int(x_extra*100))
job_name='OvCNTSymEs'+repr(int(Es/pow(10,9)))+'La'+repr(La)+'aDownUp'+repr(int(a0_down*100))+repr(int(a0_up*100))+'Ls'+repr(int(Ls*100))+'Ws'+repr(int(Ws*100))+'Wa0'+repr(int(Wa0*100))+'Wsmpl'+repr(int(Wa_sample*100))+'n'+repr(int(n*100))+'x'+repr(int(x*100))+'extraW'+repr(Wa_extra)+'L'+repr(La_extra)+'X'+repr(int(x_extra*100))
mdb.Job(name=job_name, 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=2, numDomains=2, numGPUs=0)
#Submit job
mdb.jobs[job_name].submit(consistencyChecking=OFF)
mdb.jobs[job_name].waitForCompletion()