-
Notifications
You must be signed in to change notification settings - Fork 6
/
hMesh.py
330 lines (234 loc) · 11 KB
/
hMesh.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
"""
# Mesh
This method sets the mesh controls, the element types, and specifies a global seed to make the mesh.
The method makes sure that the number of nodes is between 60,000 and 75,000 inclusive.
makeMesh(globalSeed) returns nodes and elements.
*Input:
globalSeed: Numeric value for determining seed. Default is 20.
*Output:
nodes: Integer value for number of nodes.
elements: Integer value for number of elements.
EXAMPLE:
numElm, numNodes = makeMesh(20) # Global seed is 20.
numElm, numNodes = makeMesh() # Global seed is 20.
"""
from abaqus import mdb
import abaqusConstants as aq
import mesh
def makeMesh2D(globalSeed=20):
import mesh
part.setMeshControls(elemShape=aq.TRI, regions=(allSet.faces)) # Set to triangle shape
part.setElementType(elemTypes=(mesh.ElemType(elemCode=aq.DC2D8, elemLibrary=aq.STANDARD),
mesh.ElemType(elemCode=aq.DC2D6, elemLibrary=aq.STANDARD)), regions=(allSet)) #
nodes, elements = 0, 0 # Do I need this line?
while True:
part.seedPart(deviationFactor=0.1, minSizeFactor=0.1, size=globalSeed)
part.generateMesh()
modelObject.rootAssembly.regenerate()
nodes = len(modelObject.rootAssembly.instances['PART1'].nodes)
elements = len(modelObject.rootAssembly.instances['PART1'].elements)
if nodes < 60000:
globalSeed = globalSeed * 0.80
elif nodes > 95000:
globalSeed = globalSeed * 1.20
else:
break
return elements, nodes
# WARNING: Slow convergence happens FAR too often.
def makeMesh3D(modelObject, modelRootAssembly, meshSeed=35, df=0.05):
import mesh
import random
modelRootAssembly.setMeshControls(elemShape=aq.TET, regions=(modelRootAssembly.sets['assemblyAll'].cells), technique=aq.FREE)
modelRootAssembly.setElementType(elemTypes=(mesh.ElemType(elemCode=aq.C3D20R, elemLibrary=aq.STANDARD),
mesh.ElemType(elemCode=aq.C3D15, elemLibrary=aq.STANDARD), mesh.ElemType(elemCode=aq.C3D10, elemLibrary=aq.STANDARD)),
regions=(modelRootAssembly.sets['assemblyAll']))
modelRootAssembly.setElementType(elemTypes=(mesh.ElemType(elemCode=aq.DC3D20, elemLibrary=aq.STANDARD),
mesh.ElemType(elemCode=aq.DC3D15, elemLibrary=aq.STANDARD), mesh.ElemType(elemCode=aq.DC3D10, elemLibrary=aq.STANDARD)),
regions=(modelRootAssembly.sets['assemblyAll']))
nodes, elements = 0, 0
countTot = 0
#dfs = [0.03, 0.035, 0.04, 0.045, 0.05, 0.055, 0.06, 0.065, 0.07, 0.075]
#meshs = [20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100]
#df = random.choice(dfs)
#meshSeed = random.choice(meshs)
#totalCount = 0
while True:
#totalCount = totalCount + 1
#print("number of mesh iterations " + str(totalCount))
modelRootAssembly.seedPartInstance(deviationFactor=df, minSizeFactor=0.1,
regions=(modelRootAssembly.instances['matrixFull-1'], ), size=meshSeed)
modelRootAssembly.generateMesh(regions=( modelRootAssembly.instances['matrixFull-1'], ))
modelObject.rootAssembly.regenerate()
nodes = len(modelObject.rootAssembly.instances['matrixFull-1'].nodes)
elements = len(modelObject.rootAssembly.instances['matrixFull-1'].elements)
break
"""
if elements < 35000:
newMeshSeed = [i for i in meshs if i <= 40]
newdfs = [i for i in dfs if i <= 0.05]
df = random.choice(newdfs)
meshSeed = random.choice(newMeshSeed)
elif elements < 75000:
newMeshSeed = [i for i in meshs if i <= 65]
newdfs = [i for i in dfs if i <= 0.07]
df = random.choice(newdfs)
meshSeed = random.choice(newMeshSeed)
elif elements > 300000:
newMeshSeed = [i for i in meshs if i >= 65]
newdfs = [i for i in dfs if i >= 0.05]
df = random.choice(newdfs)
meshSeed = random.choice(newMeshSeed)
elif elements > 200000:
newMeshSeed = [i for i in meshs if i >= 25]
newdfs = [i for i in dfs if i >= 0.04]
df = random.choice(newdfs)
meshSeed = random.choice(newMeshSeed)
elif elements > 162000:
meshSeed = meshSeed * 1.1
df = df*1.1
else:
break
countTot = countTot + 1
if countTot == 7:
countTot = 0
df = random.choice(dfs)
meshSeed = random.choice(meshs)
#if nodes > 100000:
# meshSeeds = [25, 27.5, 30, 32.5, 35, 37.5, 40]
# meshSeed = random.choice(meshSeeds)
# df = random.choice(dfs)
"""
return elements, nodes, df, round(meshSeed)
### FIX
# Scoping and references!!!
# Need for extracting data
def makeElementSet(fullMatrixPart, modelRootAssembly):
bottomFace = modelRootAssembly.instances['matrixFull-1'].faces.pointsOn[3]
topFace = modelRootAssembly.instances['matrixFull-1'].faces.pointsOn[1]
fullMatrixPart.Set(faces=fullMatrixPart.faces.findAt(topFace), name='top')
fullMatrixPart.Set(faces=fullMatrixPart.faces.findAt(bottomFace), name='bot')
########### Info for testing Mesh
"""
import numpy
from abaqus import mdb
import abaqusConstants as aq
from hAssembly import *
from hCoordinates import *
from hJob import *
from hMaterial import *
from hMesh import *
from hModel import *
from hPart import *
from hProperty import *
from hStep import *
materials = getMaterialList()
matrix_materials = materials.keys()
epoxyResin = materials[matrix_materials[0]]
siliconRubber = materials[matrix_materials[1]]
polystyrene = materials[matrix_materials[2]]
esbr = materials[matrix_materials[3]]
epoxyResinFillers = epoxyResin['fillers'] # 4
siliconRubberFillers = siliconRubber['fillers'] # 3
polystyreneFillers = polystyrene['fillers'] # 2
esbrFillers = esbr['fillers'] # 2
epFillers = epoxyResinFillers.keys()
srFillers = siliconRubber['fillers'].keys() # 3
polyFillers = polystyrene['fillers'].keys() # 2
eFillers = esbr['fillers'].keys() # 2
modelObject, modelName = createModel(2)
num1 = 3
num2 = 1
whichMat = esbrFillers
side, radius, portions, dP, dM, cP, cM = defExperiment(modelObject, matrix_materials[num1], eFillers[num2])
key1 = matrix_materials[num1]
key1
#key2 = eFillers[2]
key2 = eFillers[num2]
key2
radius
portions
side
phrs = whichMat[key2]['phr']
phr = phrs[4]
radius, number = invPHRAlternate3D(phr, dP, radius, dM, side)
calcPHR = round(calculatePHR3D(number, dP, radius, dM, side))
#radius, number = invVolumeAlternate3D(phr, radius, side)
#calcPHR = round(calculateVolume(number, radius, side), 3)
#delta = 0.06
#delta = 0.2
delta = 0.15
interfaceConductivity= numpy.random.sample(1) * (cP-cM) + cM # Between cM and cP
interfaceConductivity = int(interfaceConductivity[0])
#interfaceConductivity = int((cP+cM)/2.0)
# Define interface materials
defineMaterial(modelObject, "Interface", dM, interfaceConductivity)
intPortionLimit = getInterfacePortionLimit(side, radius, number, delta)
minimumPortion = 0.15
intPortionLimit - minimumPortion
## THIS IS A PROBLEM. We need Interface limit to be the highest possible portion selected,
# and the minimum must be around 0.15. Interface Limit must exceed 0.15 then.
interfacePortion = numpy.random.sample(1) * (intPortionLimit-minimumPortion) + minimumPortion # random delta to limit inclusive
interfacePortion = round(interfacePortion[0], 3)
interfacePortion
xVals, yVals, zVals = getPoints3dDeterministic(side, radius, number)
part = createMatrix(modelObject, side, False) # Create the matrix
edges1, vertices1, face1 = assignGeomSequence(part) # Create references to important sets in our matrix
part2 = createSphereParticle(modelObject, radius, side) # Create Particle part
edges2, vertices2, face2 = assignGeomSequence(part2) # Create references to important sets in particle
part3 = createSphereParticle(modelObject, (radius + radius*interfacePortion), side, "Interface") # Create interface part
edges3, vertices3, face3 = assignGeomSequence(part3) # Create references to important sets in interface
matrixSet, particleSet, interfaceSet = create3DInitialSets(part, part2, side, part3) # create sets for particle, matrix, and interface
createSection(modelObject, part, key1, matrixSet) # Create section for matrix material
createSection(modelObject, part2, key2, particleSet) # Create section for filler material
createSection(modelObject, part3, "Interface", interfaceSet) # Create section for interface
modelRootAssembly, fullMatrixPart = create3DMatrixInclusions(modelObject, part, part2, number, xVals, yVals, zVals, part3) # Create assembly and return references to assembly sets
assemblyTop, assemblyBottom, assemblyAll = define3DAssemblySets(modelRootAssembly, side)
temp1, temp2 = 328.15, 298.15 # Assign heat temperature to be used in experiment
heatStep3D(modelObject, assemblyBottom, assemblyTop, temp1, temp2) # apply heat BC
limitOutputHFL(modelObject, assemblyBottom, assemblyTop) # Limit ODB Output
elements, nodes, df, meshSeed = makeMesh3D(modelObject, modelRootAssembly, 10, 0.03750) # Draw mesh and return number of nodes and elements
makeElementSet(fullMatrixPart, modelRootAssembly)
print(delta)
print (intPortionLimit)
print(interfacePortion)
print(df)
print(meshSeed)
print(calcPHR)
calcPHR - phr
del mdb.models['Model-2']
"""
""" delta = 25 seems very safe for low errors
MESH MAP delta 0.2 with 64 particle possibility
epoxyResin Alumina seed 1 control 0.05 particles 27 elements 0.8 millionish 0 errors intertfacePortion 0.058
epoxyResin Alumina seed 0.95 control 0.04 particles 27 elements 1 millionish 3 errors interfacePortion 0.145
epoxyResin Alumina seed 0.9 control 0.0375 particles 27 elements 1.2 millionish 5 errors interfacePortion 0.145
epoxyResin Alumina seed 0.85 control 0.035 particles 27 elements 1.4 millionish 5 errors interfacePortion 0.145
epoxyResin Alumina seed 0.80 control 0.0325 particles 27 elements 1.5 millionish 5 errors interfacePortion 0.145
epoxyResin Alumina seed 0.80 control 0.0325 particles 27 elements 1.5 millionish 5 errors interfacePortion 0.145
epoxyResin Alumina seed 0.75 control 0.030 particles 27 elements 1.8 millionish 14 errors interfacePortion 0.145
epoxyResin Alumina seed 1 control 0.05 particles 27 elements 0.8 millionish 33 errors interfacePortion 0.109
epoxyResin Alumina seed 0.95 control 0.045 particles 27 elements 1 millionish 3 errors interfacePortion 0.109
epoxyResin Alumina seed 0.90 control 0.045 particles 27 elements 1.1 millionish 10 errors interfacePortion 0.109
epoxyResin Alumina seed 0.85 control 0.0425 particles 27 elements 1.3 millionish 28 errors interfacePortion 0.109
epoxyResin Alumina seed 1 control 0.05 particles 27 elements 0.8 millionish 2 errors interfacePortion 0.131
epoxyResin Alumina seed 0.95 control 0.045 particles 27 elements 1 millionish 2 errors interfacePortion 0.131
epoxyResin Alumina seed 0.90 control 0.0425 particles 27 elements 1.2 millionish 0 errors interfacePortion 0.131
epoxyResin Alumina seed 1 control 0.05 particles 27 elements 1 millionish 106700 errors interfacePortion 0.08
delta = 18 (safer.. also make interface portion 0.15)
epoxyResin Alumina seed 1 control 0.05 particles 27 elements 0.85 millionish 3 errors interfacePortion 0.1
"""
"""
elements = 100000
if elements < 35000:
print("less than 35000")
elif elements < 75000:
print("less than 75000 but greater than 35000")
elif elements > 300000:
print("greater than 300000")
elif elements > 200000:
print("greater than 200000 but less than 300000")
elif elements > 162000:
print("greater than 162000 but less than 200000")
else:
print("between 75000 and 162000")
"""