Skip to content

Commit

Permalink
Generate cohesive elements for ABAQUS mesh
Browse files Browse the repository at this point in the history 
Generate cohesive elements for ABAQUS mesh at the grain boundaries in 3D
  • Loading branch information
Nicolò Grilli authored Jul 13, 2020
0 parents commit b687246
Show file tree
Hide file tree
Showing 15 changed files with 11,765 additions and 0 deletions.
10,744 changes: 10,744 additions & 0 deletions Job-1.inp

Large diffs are not rendered by default.

39 changes: 39 additions & 0 deletions README.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
Nicolo Grilli
University of Oxford
AWE project 2020
23 Giugno 2020

Code to generate cohesive elements
at the grain boundaries in 3D

Copy 'Job-1.inp' input file in this folder
Elset must be named 'GRAIN1', 'GRAIN2', et cetera

Run
python3.6 main.py

Old bulk element list from abaqus input file will be in:
'Job-1-bulk-elems.inp'

Old node list from abaqus input file will be in:
'Job-1-node.inp'

Old element sets with the generate keyword will be in:
'Job-1-elset.inp'

New node list will be in file:
'Job-1-node-new.inp'

New bulk element list will be in file:
'Job-1-elems-new.inp'

Interface (zero thickness) element list will be in file:
'Job-1-int-elems.inp'

Maximum triple junctions are handled
A node cannot belong to more than 3 grains
Only C3D8 elements are handled
It is fundamental that nodes at the grain boundary
belongs to bulk elements in two/three grains
and are not already duplicated

Binary file added TripleJunctionTest.cae
Binary file not shown.
203 changes: 203 additions & 0 deletions TripleJunctionTest.jnl
Original file line number Diff line number Diff line change
@@ -0,0 +1,203 @@
# -*- coding: mbcs -*-
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
mdb.models['Model-1'].ConstrainedSketch(name='__profile__', sheetSize=200.0)
mdb.models['Model-1'].sketches['__profile__'].rectangle(point1=(0.0, 0.0),
point2=(10.0, 10.0))
mdb.models['Model-1'].Part(dimensionality=THREE_D, name='Part-1', type=
DEFORMABLE_BODY)
mdb.models['Model-1'].parts['Part-1'].BaseSolidExtrude(depth=0.5, sketch=
mdb.models['Model-1'].sketches['__profile__'])
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].Material(name='GRAIN1')
mdb.models['Model-1'].materials['GRAIN1'].Depvar(n=125)
mdb.models['Model-1'].materials['GRAIN1'].UserMaterial(mechanicalConstants=(
5.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0))
mdb.models['Model-1'].Material(name='GRAIN2')
mdb.models['Model-1'].materials['GRAIN2'].Depvar(n=125)
mdb.models['Model-1'].materials['GRAIN2'].UserMaterial(mechanicalConstants=(
5.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 2.0))
mdb.models['Model-1'].Material(name='GRAIN3')
mdb.models['Model-1'].materials['GRAIN3'].Depvar(n=125)
mdb.models['Model-1'].materials['GRAIN3'].UserMaterial(mechanicalConstants=(
5.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 3.0))
mdb.models['Model-1'].HomogeneousSolidSection(material='GRAIN1', name='GRAIN1',
thickness=None)
mdb.models['Model-1'].HomogeneousSolidSection(material='GRAIN2', name='GRAIN2',
thickness=None)
mdb.models['Model-1'].HomogeneousSolidSection(material='GRAIN3', name='GRAIN3',
thickness=None)
mdb.models['Model-1'].parts['Part-1'].DatumPointByCoordinate(coords=(5.0, 5.0,
0.5))
mdb.models['Model-1'].ConstrainedSketch(gridSpacing=0.7, name='__profile__',
sheetSize=28.3, transform=
mdb.models['Model-1'].parts['Part-1'].MakeSketchTransform(
sketchPlane=mdb.models['Model-1'].parts['Part-1'].faces[4],
sketchPlaneSide=SIDE1,
sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges[4],
sketchOrientation=RIGHT, origin=(5.0, 5.0, 0.5)))
mdb.models['Model-1'].parts['Part-1'].projectReferencesOntoSketch(filter=
COPLANAR_EDGES, sketch=mdb.models['Model-1'].sketches['__profile__'])
mdb.models['Model-1'].sketches['__profile__'].Line(point1=(0.0, 0.0), point2=(
5.0, 0.0))
mdb.models['Model-1'].sketches['__profile__'].HorizontalConstraint(
addUndoState=False, entity=
mdb.models['Model-1'].sketches['__profile__'].geometry[6])
mdb.models['Model-1'].sketches['__profile__'].CoincidentConstraint(
addUndoState=False, entity1=
mdb.models['Model-1'].sketches['__profile__'].vertices[6], entity2=
mdb.models['Model-1'].sketches['__profile__'].geometry[3])
mdb.models['Model-1'].sketches['__profile__'].EqualDistanceConstraint(
addUndoState=False, entity1=
mdb.models['Model-1'].sketches['__profile__'].vertices[1], entity2=
mdb.models['Model-1'].sketches['__profile__'].vertices[2], midpoint=
mdb.models['Model-1'].sketches['__profile__'].vertices[6])
mdb.models['Model-1'].parts['Part-1'].PartitionFaceBySketch(faces=
mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#10 ]',
), ), sketch=mdb.models['Model-1'].sketches['__profile__'], sketchUpEdge=
mdb.models['Model-1'].parts['Part-1'].edges[4])
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].ConstrainedSketch(gridSpacing=0.7, name='__profile__',
sheetSize=28.28, transform=
mdb.models['Model-1'].parts['Part-1'].MakeSketchTransform(
sketchPlane=mdb.models['Model-1'].parts['Part-1'].faces[4],
sketchPlaneSide=SIDE1,
sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges[11],
sketchOrientation=RIGHT, origin=(5.0, 5.0, 0.5)))
mdb.models['Model-1'].parts['Part-1'].projectReferencesOntoSketch(filter=
COPLANAR_EDGES, sketch=mdb.models['Model-1'].sketches['__profile__'])
mdb.models['Model-1'].sketches['__profile__'].Line(point1=(0.0, 0.0), point2=(
5.0, 5.0))
mdb.models['Model-1'].sketches['__profile__'].CoincidentConstraint(
addUndoState=False, entity1=
mdb.models['Model-1'].sketches['__profile__'].vertices[7], entity2=
mdb.models['Model-1'].sketches['__profile__'].geometry[7])
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].ConstrainedSketch(gridSpacing=0.7, name='__profile__',
sheetSize=28.28, transform=
mdb.models['Model-1'].parts['Part-1'].MakeSketchTransform(
sketchPlane=mdb.models['Model-1'].parts['Part-1'].faces[4],
sketchPlaneSide=SIDE1,
sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges[11],
sketchOrientation=RIGHT, origin=(5.0, 5.0, 0.5)))
mdb.models['Model-1'].parts['Part-1'].projectReferencesOntoSketch(filter=
COPLANAR_EDGES, sketch=mdb.models['Model-1'].sketches['__profile__'])
mdb.models['Model-1'].sketches['__profile__'].Line(point1=(0.0, 0.0), point2=(
5.0, 5.0))
mdb.models['Model-1'].sketches['__profile__'].CoincidentConstraint(
addUndoState=False, entity1=
mdb.models['Model-1'].sketches['__profile__'].vertices[7], entity2=
mdb.models['Model-1'].sketches['__profile__'].geometry[7])
mdb.models['Model-1'].parts['Part-1'].PartitionFaceBySketch(faces=
mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#10 ]',
), ), sketch=mdb.models['Model-1'].sketches['__profile__'], sketchUpEdge=
mdb.models['Model-1'].parts['Part-1'].edges[11])
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].ConstrainedSketch(gridSpacing=0.7, name='__profile__',
sheetSize=28.28, transform=
mdb.models['Model-1'].parts['Part-1'].MakeSketchTransform(
sketchPlane=mdb.models['Model-1'].parts['Part-1'].faces[0],
sketchPlaneSide=SIDE1,
sketchUpEdge=mdb.models['Model-1'].parts['Part-1'].edges[4],
sketchOrientation=RIGHT, origin=(3.333333, 4.333333, 0.5)))
mdb.models['Model-1'].parts['Part-1'].projectReferencesOntoSketch(filter=
COPLANAR_EDGES, sketch=mdb.models['Model-1'].sketches['__profile__'])
mdb.models['Model-1'].sketches['__profile__'].Line(point1=(-0.666667, 1.666667)
, point2=(4.333333, -3.333333))
mdb.models['Model-1'].sketches['__profile__'].CoincidentConstraint(
addUndoState=False, entity1=
mdb.models['Model-1'].sketches['__profile__'].vertices[7], entity2=
mdb.models['Model-1'].sketches['__profile__'].geometry[2])
mdb.models['Model-1'].parts['Part-1'].PartitionFaceBySketch(faces=
mdb.models['Model-1'].parts['Part-1'].faces.getSequenceFromMask(('[#1 ]',
), ), sketch=mdb.models['Model-1'].sketches['__profile__'], sketchUpEdge=
mdb.models['Model-1'].parts['Part-1'].edges[4])
del mdb.models['Model-1'].sketches['__profile__']
mdb.models['Model-1'].parts['Part-1'].PartitionCellByExtrudeEdge(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#1 ]',
), ), edges=(mdb.models['Model-1'].parts['Part-1'].edges[1], ), line=
mdb.models['Model-1'].parts['Part-1'].edges[10], sense=REVERSE)
mdb.models['Model-1'].parts['Part-1'].PartitionCellByExtrudeEdge(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#1 ]',
), ), edges=(mdb.models['Model-1'].parts['Part-1'].edges[11], ), line=
mdb.models['Model-1'].parts['Part-1'].edges[17], sense=REVERSE)
mdb.models['Model-1'].parts['Part-1'].PartitionCellByExtrudeEdge(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#2 ]',
), ), edges=(mdb.models['Model-1'].parts['Part-1'].edges[11], ), line=
mdb.models['Model-1'].parts['Part-1'].edges[1], sense=FORWARD)
mdb.models['Model-1'].parts['Part-1'].Set(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#1 ]',
), ), name='GRAIN1')
mdb.models['Model-1'].parts['Part-1'].SectionAssignment(offset=0.0,
offsetField='', offsetType=MIDDLE_SURFACE, region=
mdb.models['Model-1'].parts['Part-1'].sets['GRAIN1'], sectionName='GRAIN1',
thicknessAssignment=FROM_SECTION)
mdb.models['Model-1'].parts['Part-1'].Set(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#4 ]',
), ), name='GRAIN2')
mdb.models['Model-1'].parts['Part-1'].SectionAssignment(offset=0.0,
offsetField='', offsetType=MIDDLE_SURFACE, region=
mdb.models['Model-1'].parts['Part-1'].sets['GRAIN2'], sectionName='GRAIN2',
thicknessAssignment=FROM_SECTION)
mdb.models['Model-1'].parts['Part-1'].Set(cells=
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#2 ]',
), ), name='GRAIN3')
mdb.models['Model-1'].parts['Part-1'].SectionAssignment(offset=0.0,
offsetField='', offsetType=MIDDLE_SURFACE, region=
mdb.models['Model-1'].parts['Part-1'].sets['GRAIN3'], sectionName='GRAIN3',
thicknessAssignment=FROM_SECTION)
mdb.models['Model-1'].parts['Part-1'].setElementType(elemTypes=(ElemType(
elemCode=C3D8, elemLibrary=STANDARD, secondOrderAccuracy=OFF,
distortionControl=DEFAULT), ElemType(elemCode=C3D6, elemLibrary=STANDARD),
ElemType(elemCode=C3D4, elemLibrary=STANDARD)), regions=(
mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask(('[#7 ]',
), ), ))
mdb.models['Model-1'].parts['Part-1'].setMeshControls(algorithm=ADVANCING_FRONT
, regions=mdb.models['Model-1'].parts['Part-1'].cells.getSequenceFromMask((
'[#7 ]', ), ), technique=SWEEP)
mdb.models['Model-1'].parts['Part-1'].seedPart(deviationFactor=0.1,
minSizeFactor=0.1, size=0.5)
mdb.models['Model-1'].parts['Part-1'].generateMesh()
# Save by engs1992 on 2020_06_05-16.07.50; build 2017 2016_09_27-22.54.59 126836
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
mdb.models['Model-1'].rootAssembly.DatumCsysByDefault(CARTESIAN)
mdb.models['Model-1'].rootAssembly.Instance(dependent=ON, name='Part-1-1',
part=mdb.models['Model-1'].parts['Part-1'])
mdb.models['Model-1'].StaticStep(name='Step-1', nlgeom=ON, previous='Initial')
mdb.models['Model-1'].rootAssembly.Set(faces=
mdb.models['Model-1'].rootAssembly.instances['Part-1-1'].faces.getSequenceFromMask(
('[#800 ]', ), ), name='Set-1')
mdb.models['Model-1'].DisplacementBC(amplitude=UNSET, createStepName='Step-1',
distributionType=UNIFORM, fieldName='', fixed=OFF, localCsys=None, name=
'BC-1', region=mdb.models['Model-1'].rootAssembly.sets['Set-1'], u1=0.0,
u2=0.0, u3=0.0, ur1=UNSET, ur2=UNSET, ur3=UNSET)
mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF,
explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF,
memory=90, memoryUnits=PERCENTAGE, model='Model-1', modelPrint=OFF,
multiprocessingMode=DEFAULT, name='Job-1', nodalOutputPrecision=SINGLE,
numCpus=1, numGPUs=0, queue=None, resultsFormat=ODB, scratch='', type=
ANALYSIS, userSubroutine='', waitHours=0, waitMinutes=0)
# Save by engs1992 on 2020_06_05-16.09.17; build 2017 2016_09_27-22.54.59 126836
2 changes: 2 additions & 0 deletions __init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# this file tells python
# that all classes in this folder can be imported
97 changes: 97 additions & 0 deletions abaqusparser.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,97 @@
# Nicolo Grilli
# University of Oxford
# AWE project 2020
# 23 Giugno 2020

# class representing abaqus input file
# and functions for parsing

class AbaqusParser:

# prefissofile is the prefix of abaqus .inp file
def __init__(self,prefissofile):
self.prefissofile = prefissofile
self.filename = prefissofile + '.inp'

# read element sets from file
# only generate keyword is accepted
def ReadElsets(self):
fid = open(self.filename,'r')
elsetsfilename = self.prefissofile + '-elsets.inp'
fout = open(elsetsfilename,'w')

flagprintnextline = False # next line should be printed or not

for line in fid:
if (flagprintnextline):
fout.write(line)
flagprintnextline = False
if (line[0:6] == '*Elset'): # this is an elset
posgrain = line.find('GRAIN') # -1 if this is not a grain
posgenerate = line.find('generate') # generate is only format accepted
if (posgrain >= 0): # this is a grain
if (posgenerate >= 0):
fout.write(line)
flagprintnextline = True # print next line containing the elements start/end
else:
print('generate keyword is not present' + '\n')
print('generate keyword is the only accepted format' + '\n')
exit()

fid.close()
fout.close()

return 1

# read nodes from file
def ReadNodes(self):
fid = open(self.filename,'r')
nodesfilename = self.prefissofile + '-node.inp'
fout = open(nodesfilename,'w')

flagfoundnodes = False # block with nodes has been found

for line in fid:
if (flagfoundnodes): # this line may contain nodes
if (line[0] == '*'): # reached the end of nodes
flagfoundnodes = False
break # needed because you may find also *Node Output
else:
fout.write(line)
if (line[0:5] == '*Node'): # these are the nodes
flagfoundnodes = True

fid.close()
fout.close()

return 1

# read bulk elements from file
def ReadBulkElems(self):
fid = open(self.filename,'r')
elemsfilename = self.prefissofile + '-bulk-elems.inp'
fout = open(elemsfilename,'w')

flagfoundelems = False # block with elements has been found

for line in fid:
if (flagfoundelems): # this line may contain elements
if (line[0] == '*'): # reached the end of elements
flagfoundelems = False
break # needed because you may find also *Element Output
else:
fout.write(line)
if (line[0:8] == '*Element'): # these are the elements
flagfoundelems = True

fid.close()
fout.close()

return 1







48 changes: 48 additions & 0 deletions cohelement.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
# Nicolo Grilli
# University of Oxford
# AWE project 2020
# 5 giugno 2020

# class representing a COH3D8 element
# at the boundary between grains

import numpy as np
from interfelement import InterfElement

class CohElement:

def __init__(self,elem1,elem2):
# a cohesive elements is at the interface
# of two interface elements
self.elem1 = elem1
self.elem2 = elem2
# nodes
self.nodes = np.zeros(shape=(8))
# index must follow the indices of bulk elements
self.index = 0

# assign nodes based on the input interface elements
def AssignNodes(self):
# follow the order of the first interface element
# so first interface element face will be counterclockwise
# and second interface element face will be clockwise
temp1order = np.int_(self.elem1.cohnodeorder)
temp1corr = self.elem1.corrnodes
tempnodes = np.zeros(shape=(8))
for n in range(0,4):
tempnodes[n] = int(self.elem1.nodes[temp1order[n]-1])
for search in range(0,4):
if (temp1corr[search][0] == temp1order[n]):
tempindex2 = int(temp1corr[search][1])
tempnodes[n+4] = int(self.elem2.nodes[tempindex2-1])
self.nodes = np.int_(tempnodes)
return 1









Loading

0 comments on commit b687246

Please sign in to comment.