MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
FoamGeometryConstruction_Periodic.py
Go to the documentation of this file.
1 
13 import os
14 import os.path
15 import numpy as np
16 import re
17 import time
18 import random
19 import math
20 
21 mypath=os.getcwd()
22 
30 def main(MU,SIGMA,NumOfCells,filenameOut,packing,alternativePackingAlgorithm,
31  tesselation,visualizeTesselation,geometry,statistics,hypermesh,deleteFiles,
32  dx,dy,dz):
33 
34  if packing:
35  NumSpheres=NumOfCells
36  Rad=abs((np.random.normal(MU, SIGMA,NumSpheres)))
37  Rads1=list(range(NumSpheres))
38  t=0
39  for i in range(NumSpheres):
40  c=abs(Rad[t])
41  Rads1[t]=c.astype(np.float)
42  t=t+1
43  Rads1=sorted(Rads1)
44  v=0.00
45  for i in range(NumSpheres):
46  v=v+((2.00*Rads1[i])**3.00)
47  centers = [[0 for i in range(3)] for j in range(NumSpheres)]
48  v=v*1.40
49  lc=v**(1.00/3.00)
50  K=0
51  while K==0:
52  j=-1
53  h=0
54  timeout = time.time() + 10
55  while NumSpheres>=j and h==0:
56  if time.time()>timeout:
57  h=1
58  break
59  j=j+1
60  if j==NumSpheres:
61  K=1
62  break
63  PickCenterX,PickCenterY,PickCenterZ=\
64  lc*random.random(),lc*random.random(),lc*random.random()
65  while (lc-Rads1[j]>=PickCenterX and lc-Rads1[j]>=PickCenterY
66  and lc-Rads1[j]>=PickCenterZ and Rads1[j]<PickCenterX
67  and Rads1[j]<PickCenterY and Rads1[j]<PickCenterZ):
68  PickCenterX,PickCenterY,PickCenterZ=\
69  lc*random.random(),lc*random.random(),lc*random.random()
70  centers[j][0],centers[j][1],centers[j][2]=\
71  PickCenterX,PickCenterY,PickCenterZ
72  KeepCentreX, KeepCentreY, KeepCentreZ, KeepR=\
73  PickCenterX, PickCenterY, PickCenterZ, Rads1[j]
74  if j>0:
75  for t in range(0,j):
76  if (((( ((KeepCentreX-centers[t][0])**2.00)+
77  ((KeepCentreY-centers[t][1])**2.00)+
78  ((KeepCentreZ-centers[t][2])**2.00))**0.50)-
79  (KeepR+Rads1[t]))<0.000) and t!=j:
80  centers[j][0],centers[j][0],centers[j][0]=0,0,0
81  j=j-1
82  break
83  mypath=os.getcwd()
84  myfile=os.path.join(mypath,'Project01.rco')
85  f=open(myfile,'w')
86  for i in range(NumSpheres):
87  f.write('{0:f}{1}{2:f}{3}{4:f}{5}{6:f}\n'.format(
88  centers[i][0],' ',centers[i][1],' ',centers[i][2],' ',
89  2.0*Rads1[i]))
90  f.close()
91  MAXcenters=max(centers)
92  Mincenters=min(centers)
93  EdgeCubeSize=[math.ceil(MAXcenters[0]-Mincenters[0]),
94  math.ceil(MAXcenters[1]-Mincenters[1]),
95  math.ceil(MAXcenters[2]-Mincenters[2])]
96  EdgeRVESize=int(3.0*max(EdgeCubeSize)) #For NEPER: Size of edge of RVE
97  if alternativePackingAlgorithm:
98  os.system('./spherepack')
99  if tesselation:
100  myfile12=os.path.join(mypath,'Project01.rco')
101  CentersRads=np.loadtxt(myfile12,usecols = (0,1,2,3))
102  Centers=CentersRads[:,:3] #All centers of spheres
103  Rads=CentersRads[:,3] #All radii of spheres
104  Rads=Rads/2
105  MaxCenters=np.amax(Centers, axis=0)
106  L2=int(0.5+MaxCenters[0])
107  print (L2)
108  X=L2+Centers[:,0]
109  Y=L2+Centers[:,1]
110  Z=L2+Centers[:,2]
111  Centers27=np.array([
112  X,X+L2,X-L2,X+L2,X,X+L2,X,X-L2,X-L2,X-L2,X+L2,X,X,X-L2,X+L2,
113  X+L2,X,X,X-L2,X,X,X-L2,X+L2,X+L2,X-L2,X-L2,X+L2,Y,Y+L2,Y-L2,
114  Y+L2,Y+L2,Y,Y-L2,Y,Y-L2,Y+L2,Y-L2,Y-L2,Y+L2,Y,Y,Y,Y+L2,Y,Y,
115  Y-L2,Y,Y+L2,Y-L2,Y+L2,Y-L2,Y+L2,Y-L2,Z,Z+L2,Z-L2,Z,Z+L2,Z+L2,
116  Z-L2,Z-L2,Z,Z,Z,Z+L2,Z-L2,Z+L2,Z-L2,Z,Z,Z+L2,Z,Z,Z-L2,Z+L2,
117  Z+L2,Z-L2,Z+L2,Z-L2,Z-L2])
118  # Creation of input .txt files for neper
119  myfile3=os.path.join(mypath,'Centers.txt')
120  ff=open(myfile3,'w')
121  for i in range(0,27):
122  for j in range(0,NumOfCells):
123  ff.write('{0:f}\t{1:f}\t{2:f}\n'.format(Centers27[i,j],
124  Centers27[i+27,j],Centers27[i+54,j]))
125  ff.close()
126  myfile4=os.path.join(mypath,'Rads.txt')
127  fff=open(myfile4,'w')
128  for i in range(0,27):
129  for j in range(0,NumOfCells):
130  fff.write('{0:f}\n'.format(Rads[j]))
131  fff.close()
132  commandTessellation="neper -T -n {0:d} -domain \
133  'cube({1:d},{2:d},{3:d})' -morpho @Centers.txt -weight @Rads.txt \
134  -regularization 1 -mloop 15 \
135  -o RVE27 -format tess,geo \
136  -statcell vol -statedge length -statface area \
137  -statver x".format((27*NumSpheres),EdgeRVESize,EdgeRVESize,
138  EdgeRVESize)
139  os.system(commandTessellation)
140  if visualizeTesselation: # needs POV-Ray
141  commandVisualization="neper -V RVE27.tess -datacellcol ori \
142  -datacelltrs 0.5 -showseed all -dataseedrad @Rads.txt \
143  -dataseedtrs 1.0 -print RVE27"
144  os.system(commandVisualization)
145 
148  if geometry:
149  myfile12=os.path.join(mypath,'RVE27.stver')
150  verX=np.loadtxt(myfile12)
151  NumOfNodes=len(verX)
152  myfile12=os.path.join(mypath,'RVE27.stedge')
153  EdgesLength=np.loadtxt(myfile12)
154  NumOfEdges=len(EdgesLength)
155  myfile12=os.path.join(mypath,'RVE27.stface')
156  SurfacesArea=np.loadtxt(myfile12)
157  NumOfSurfaces=len(SurfacesArea)
158  myfile12=os.path.join(mypath,'RVE27.stcell')
159  CellVolumes=np.loadtxt(myfile12)
160  NumOfVolumes=len(CellVolumes)
161 
162  myfile12=os.path.join(mypath,'RVE27.geo')
163  text_file = open(myfile12, "r") lines = text_file.readlines()
164 
165  t=0
166  Nodes=list(range(NumOfNodes))
167  for i in range(0,NumOfNodes):
168  currentline = lines[i].split("{")
169  a=currentline[1].split("}")
170  b=a[0].split(",")
171  c=np.array(b)
172  Nodes[t]=np.absolute(c.astype(np.float))
173  t=t+1
174 
175  Edges=list(range(NumOfEdges))
176  r=0
177  t=0
178  for i in range(NumOfNodes,NumOfNodes+NumOfEdges):
179  currentline = lines[i].split("{")
180  a=currentline[1].split("}")
181  b=a[0].split(",")
182  c=np.array(b)
183  Edges[t]=np.absolute(c.astype(np.float))
184  t=t+1
185 
186  Faces=list(range(NumOfSurfaces))
187  r=0
188  t=0
189  for i in range(0,NumOfSurfaces):
190  j=NumOfNodes+NumOfEdges+(2*i)
191  currentline = lines[j].split("{")
192  a=currentline[1].split("}")
193  b=a[0].split(",")
194  c=np.array(b)
195  Faces[t]=np.absolute(c.astype(np.float))
196  t=t+1
197 
198  Volumes=list(range(NumOfVolumes))
199  #print NumOfNodes+NumOfEdges+(2*NumOfSurfaces)
200  a=list(lines[NumOfNodes+NumOfEdges+(2*NumOfSurfaces)])
201  a[0:21]=''
202  o=len(a)
203  a[o-3:]=''
204  lines[NumOfNodes+NumOfEdges+(2*NumOfSurfaces)]="".join(a)
205  currentline=lines[NumOfNodes+NumOfEdges+(2*NumOfSurfaces)].split(",")
206  c=np.array(currentline)
207  Volumes[0]=np.absolute(c.astype(np.float))
208  t=1
209  for i in range(1,NumOfVolumes):
210  j=NumOfNodes+NumOfEdges+(2*NumOfSurfaces)+i
211  currentline = lines[j].split(" = {")
212  a=currentline[2].split("}")
213  b=a[0].split(",")
214  c=np.array(b)
215  Volumes[t]=np.absolute(c.astype(np.float))
216  t=t+1
217 
218  MAX0=list(range(0,NumOfCells))
219  for i in range(NumOfCells):
220  MAX0[i]=max(Volumes[i])
221  MaxIndexOfFaces=max(MAX0)
222  MAX1=list(range(0,int(MaxIndexOfFaces)))
223  for i in range(int(MaxIndexOfFaces)):
224  MAX1[i]=max(Faces[i])
225  MaxIndexOfEdges=max(MAX1)
226  MAX2=list(range(0,int(MaxIndexOfEdges)))
227  for i in range(int(MaxIndexOfEdges)):
228  MAX2[i]=max(Edges[i])
229  MaxIndexOfNodes=max(MAX2)
230  text_file.close()
231 
233  myfile13=os.path.join(mypath,filenameOut+".geo")
234  text_GEO = open(myfile13, "w")
235  for i in range(int(MaxIndexOfNodes)):
236  text_GEO.write('{0}'.format(lines[i]))
237  for i in range(int(NumOfNodes),int(NumOfNodes+MaxIndexOfEdges)):
238  text_GEO.write('{0}'.format(lines[i]))
239  j=int(NumOfNodes+NumOfEdges)
240  for i in range(int(NumOfNodes+NumOfEdges),
241  int(NumOfNodes+NumOfEdges+MaxIndexOfFaces)):
242  text_GEO.write('{0}{1}'.format(lines[j],lines[j+1]))
243  j=2+j
244  text_GEO.write('{0}{1}{2}\n'.format(' Surface Loop (1) = {',
245  lines[int(NumOfNodes+NumOfEdges+(2*NumOfSurfaces))],'};'))
246  NumOfCells=NumOfVolumes/27
247  j=int(NumOfNodes+NumOfEdges+(2*NumOfSurfaces)+1)
248  for i in range(int(NumOfCells-1)):
249  text_GEO.write('{0}'.format(lines[j]))
250  j+=1
251  text_GEO.write('{0}{1}{2}{3}{4}'.format(
252  'Volume (',NumOfCells,') = {',NumOfCells,'};'))
253  text_GEO.close()
254 
256  myfile13=os.path.join(mypath,filenameOut+".gnu")
257  text_gnu = open(myfile13, "w")
258  for i in range(int(NumOfNodes),int(NumOfNodes+MaxIndexOfEdges)):
259  a=re.findall("[-+]?\d+[\.]?\d*", lines[i])
260  b=int(a[1])-1
261  c=int(a[2])-1
262  text_gnu.write('{0}\t{1}\t{2}\n'.format(
263  (Nodes[b][0]-L2)/L2,(Nodes[b][1]-L2)/L2,(Nodes[b][2]-L2)/L2))
264  text_gnu.write('{0}\t{1}\t{2}\n\n\n'.format(
265  (Nodes[c][0]-L2)/L2,(Nodes[c][1]-L2)/L2,(Nodes[c][2]-L2)/L2))
266  text_gnu.close()
267  # End of Geometry construction
268 
271  if statistics:
272  myfile14=os.path.join(mypath,'EdgeLengths_PeriodicRVE.txt')
273  EdgesLengthsOfMiddleRVE= open(myfile14, "w")
274  for i in range(int(MaxIndexOfEdges)):
275  EdgesLengthsOfMiddleRVE.write('{0}\n'.format(EdgesLength[i]))
276  EdgesLengthsOfMiddleRVE.close()
277  myfile15=os.path.join(mypath,'FaceAreas_PeriodicRVE.txt')
278  AreaOfFaceOfMiddleRVE= open(myfile15, "w")
279  for i in range(int(MaxIndexOfFaces)):
280  AreaOfFaceOfMiddleRVE.write('{0}\n'.format(SurfacesArea[i]))
281  AreaOfFaceOfMiddleRVE.close()
282  myfile16=os.path.join(mypath,'CellVolumes_PeriodicRVE.txt')
283  VolumeOfCellOfMiddleRVE= open(myfile16, "w")
284  for i in range(int(NumOfCells)):
285  VolumeOfCellOfMiddleRVE.write('{0}\n'.format(CellVolumes[i]))
286  VolumeOfCellOfMiddleRVE.close()
287 
290  if hypermesh:
291  myfile17=os.path.join(mypath,'HyperMeshInput.cmf')
292  HyperMeshInput= open(myfile17, "w")
293  HyperMeshInput.write('{0}\n'.format('*cleanuptoleranceset(0.001)'))
294  HyperMeshInput.write('{0}\n'.format('*toleranceset(0.001)'))
295  for i in range(int(MaxIndexOfNodes)):
296  HyperMeshInput.write('{0}{1}{2}{3}{4}{5}{6}\n'.format(
297  '*createnode(',Nodes[i][0],',',Nodes[i][1],',',Nodes[i][2],
298  ',0,0,0)'))
299  for i in range(int(MaxIndexOfEdges)):
300  HyperMeshInput.write('{0}\n{1}{2}{3}{4}\n'.format(
301  '*linecreatefromnodes(1,0,150,5,179)','*createlist(nodes,1) ',
302  int(Edges[i][0]),' ',int(Edges[i][1])))
303  for i in range(int(MaxIndexOfFaces)):
304  HyperMeshInput.write('{0}\n'.format('*surfacemode(4)'))
305  HyperMeshInput.write('{0}'.format('*createmark(lines,1) '))
306  for j in range(len(Faces[i])):
307  HyperMeshInput.write('{0}\t'.format(int(Faces[i][j])))
308  HyperMeshInput.write('\n{0}\n{1}\n'.format(
309  '*createplane(1,1,0,0,0,0,0)','*splinesurface(lines,1,1,1,1)'))
310  HyperMeshInput.write('{0}\n'.format('*createmark(lines,1) "all"'))
311  HyperMeshInput.write('{0}\n'.format('*deletemark(lines,1)'))
312  HyperMeshInput.write('{0}\n'.format('*settopologyedgedisplay(2,0)'))
313  HyperMeshInput.write('{0}\n'.format('*plot()'))
314  HyperMeshInput.write('{0}\n'.format('*settopologyedgedisplay(3,0)'))
315  HyperMeshInput.write('{0}\n'.format('*plot()'))
316  HyperMeshInput.write('{0}\n'.format('*settopologyedgedisplay(1,0)'))
317  HyperMeshInput.write('{0}\n'.format('*plot()'))
318  Tol=0.001
319  for i in range(NumOfCells):
320  HyperMeshInput.write('{0}\n'.format(
321  '*createmark(surfaces,1) "displayed"'))
322  HyperMeshInput.write('{0}{1}{2}{3}{4}\n'.format(
323  '*selfstitchcombine(1,82,',Tol,',',Tol,')'))
324  Tol=0.0001+Tol
325  HyperMeshInput.close()
326 
329  if deleteFiles:
330  os.chdir(mypath)
331  os.remove('Project01.prj')
332  os.remove('Project01.rco')
333  os.remove('Project01.rst')
334  os.remove('Project01.sco')
335  os.remove('Rads.txt')
336  os.remove('Centers.txt')
337  os.remove('RVE27.geo')
338  os.remove('RVE27.stcell')
339  os.remove('RVE27.stedge')
340  os.remove('RVE27.stface')
341  os.remove('RVE27.stver')
342  os.remove('RVE27.tess')
343 
FoamConstruction.FoamGeometryConstruction_Periodic.main
def main(MU, SIGMA, NumOfCells, filenameOut, packing, alternativePackingAlgorithm, tesselation, visualizeTesselation, geometry, statistics, hypermesh, deleteFiles, dx, dy, dz)
Main function.
Definition: FoamGeometryConstruction_Periodic.py:32