new folder

Lattice_Cell_script/MicroSample.py

+#coding-Utf-9-*-
+import numpy as np
+import os
+import subprocess
+import pickle
+import string
+import time
+import sys
+#from MF_inverse import *
+from Auto_geometry_parametric import *
+import math as m
+import vmf_test
+import vmf_hollow_OD
+import vmf_hollow_ID
+import vmf_bisection
+import pandas as pd
+
+from scipy.spatial import distance
+from scipy.optimize import minimize
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
+import matplotlib.pyplot as plt 
+
+from Auto_geometry_parametric import *
+
+###################################################
+vf_tpye=["a", "b"]
+GeoID = []
+GeoID.append(["2", "3", "4", "5", "6", "7", "8", "9", "91"])
+GeoID.append(["2", "3", "4", "5", "6", "7", "8", "9", "91", "92"])
+cell_num = np.array([1,1,1])
+
+####################################################
+Nsim=200
+i=0
+j=1
+vfr= vf_tpye[i]
+cell_name= GeoID[i][j] 
+
+rho=str(128)       # defaut
+length = cell_num[0]
+width = cell_num[1]
+height = cell_num[2]
+
+
+####################################################
+if cell_name!="8":
+  VmFr = 0.3*np.random.uniform()
+if cell_name=="8":
+  VmFr=0.15*np.random.uniform()
+sb = np.random.uniform(0.5, 2.5)  #Set size of unit cell
+Radius = np.random.uniform()*sb
+Geo, Radi = Geo_generator(cell_name, rho, vfr, length, width, height, sb, Radius, VmFr)
+V_R_S = [VmFr, Radi, sb]
+
+print("Create geo file ", Geo, " cell type ", cell_name, " of volume fraction type", vfr, "with as results: volume fraction ", VmFr," Radius ", Radi)  
+#runInPath(folder, fname, Geo, V_R_S)
+
+#os.system("rm E_* *.geo *.msh")

Lattice_Cell_script/Readme.txt

+generator of lattice cells
+--------------------------
+
+1) Compile the f code with ./compile.sh
+2) geometry_parametric.py allows defining manually a cell, instruction in TUTORIAL.pdf
+3) MicroSample.py is a automatic generator of random geometries (as an example, one geo of type 3 and random volume fraction is generated)

Lattice_Cell_script/compile.sh

+f2py3 -c vmf_bisection.f95 -m vmf_bisection
+f2py3 -c vmf_hollow_ID.f95 -m vmf_hollow_ID
+f2py3 -c vmf_hollow_OD.f95 -m vmf_hollow_OD
+f2py3 -c vmf_test.f95 -m vmf_test
+

Lattice_Cell_script/utils_geo.py

+import numpy as np
+import math as m
+
+import vmf_test
+import vmf_hollow_OD
+import vmf_hollow_ID
+import vmf_bisection
+
+from scipy.spatial import distance
+
+from scipy.optimize import minimize
+
+from mpl_toolkits.mplot3d import Axes3D
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
+import matplotlib.pyplot as plt
+
+import time
+import sys
+import os
+
+def voxelized(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B):
+######
+# Create voxelized mesh -- CALCULATES VOLUME FRACTION
+######
+#####Voxelized mesh#####
+#number of voxels in the unit cell
+#to avoid loosing resolution in any direction, int(x*lenght_cell)
+
+    nx=int(int(rho)*x_cell)
+    ny=int(int(rho)*y_cell)
+    nz=int(int(rho)*z_cell)
+    
+    #Vox = open ("VOX"+name+"_unitcell_size_"+str(sb)+"_length_"+str(length)+"_width_"+str(width)+"_height_"+str(height)+".geo", 'w')
+    #Vox = open ("vox"+"_"+name+"_"+rho+".geo",'w')
+    Vox = open ("vox"+"_"+rho+".geo",'w')
+    Phase_map = open ("phase_map"+"_"+rho+".geo",'w')
+    Vox.write(str(x_cell*length) + " " + str(y_cell*width) + " " + str(z_cell*height) + " " + str(nx*length) + " " + str(ny*width) + " " + str(nz*height) + "\n")
+    Phase_map.write(str(x_cell*length) + " " + str(y_cell*width) + " " + str(z_cell*height) + " " + str(nx*length) + " " + str(ny*width) + " " + str(nz*height) + "\n")
+    start_time = time.time()     
+    
+    if name=="92":
+        vox_mesh_fort_OD=vmf_hollow_OD.vox_mesh_fort(A,nx,ny,nz,sb,len(A))
+        vox_mesh_fort_ID=vmf_hollow_ID.vox_mesh_fort2(A,nx,ny,nz,sb,len(A))
+        vox_mesh_fort= vox_mesh_fort_OD-vox_mesh_fort_ID    
+    else:
+        if vfr=="a" and name!="8":
+            #"6.a. Set volume fraction from 0 to 1"
+            vox_mesh_fort=vmf_bisection.vox_mesh_fort(A,nx,ny,nz,x_cell,y_cell,z_cell,vfraction,len(A))
+        if vfr=="a" and name=="8":
+            # "6.a. Set volume fraction from 0 to 0.15"
+            vox_mesh_fort=vmf_bisection.vox_mesh_fort(A,nx,ny,nz,x_cell,y_cell,z_cell,vfraction,len(A))    	
+        if vfr=="b":
+            vox_mesh_fort=vmf_test.vox_mesh_fort(A,nx,ny,nz,x_cell,y_cell,z_cell,len(A))         
+
+        
+    
+    print("Temps d execution : %s secondes ---" % (time.time() - start_time))
+    
+    start_time = time.time()
+        
+    for i in range (nz*height):
+        for j in range (ny*width):
+            for k in range  (nx*length):
+                x=k%nx
+                y=j%ny
+                z=i%nz
+                #Vox.write(str(vox_mesh[(x,y,z)]) + " ")
+                if vox_mesh_fort[x][y][z]!=0:    
+                    Vox.write(str(1) + " ")
+                else:
+                    Vox.write(str(0) + " ") 
+                Phase_map.write(str(vox_mesh_fort[(x,y,z)]) + " ")
+            Vox.write("\n")
+            Phase_map.write("\n")
+    Vox.close()
+    Phase_map.close()
+    print("youhou")
+    print("Temps d execution : %s secondes ---" % (time.time() - start_time))
+
+
+
+def mesh(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B, r0):
+#####-------------------------------------
+
+#to create the geo file
+
+#####-------------------------------------
+    # print ("specify radius of struts, from 0 to "+ str(sb))
+    # r0=input()
+    
+    
+    if name=="92":
+        Union1 = "BooleanUnion("+str(len(A)+2)+") = {"
+        Intersection1 ="BooleanIntersection("+str(len(A)+3)+")={ Volume{"+str(len(A)+2)+"}; Delete;}{Volume{1}; Delete;};"
+        Union2 = "BooleanUnion("+str(2*len(A)+5)+") = {"
+        Intersection2 ="BooleanIntersection("+str(2*len(A)+6)+")={ Volume{"+str(2*len(A)+5)+"}; Delete;}{Volume{"+str(len(A)+4)+"}; Delete;};"
+        Difference="BooleanDifference("+str(2*len(A)+7)+") = { Volume{"+str(len(A)+3)+"}; Delete; }{ Volume{"+str(2*len(A)+6)+"}; Delete; };"
+        #dimension of the lattice/ NUMBER OF UNIT CELL IN EACH DIRECTION
+        #print ("number of unit cells in each direction")
+        #print ("lenght")
+        #length = int(input())
+        #print ("width")
+        #width = int(input())
+       # print("height")
+       # height = int(input())
+        
+        #Lattice    = open(name+"_unitcell_size_"+str(sb)+"_length_"+str(length)+"_width_"+str(width)+"_height_"+str(height)+".geo", 'w')
+        #Lattice = open (name+".geo",'w')
+        Lattice = open (name+"_"+rho+".geo",'w')
+        Lattice.write( "SetFactory('OpenCASCADE');" + "\n" 
+                        + "Box(1) = {0,0,0,"+str(x_cell)+","+str(y_cell)+","+str(z_cell)+"};"+ "\n" )
+        for i in range (len(A)):
+            a=A[i][0][0]
+            b=A[i][0][1]
+            c=A[i][0][2]
+            d=A[i][1][0]-A[i][0][0]
+            e=A[i][1][1]-A[i][0][1]
+            f=A[i][1][2]-A[i][0][2]
+            r=A[i][2][0]
+            Lattice.write("Cylinder("+str(i+2)+") = {"+ str(a)+","+ str(b)+","+ str(c)+","+ str(d)+","+ str(e)+","+ str(f)+","+ str(r)+",2*Pi};"+ "\n")
+            if i in range (0,8):	                 	
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")  
+            if i != len(A)-1:
+                Union1+="Volume{"+str(i+2)+"};"
+            else:
+                Union1+="Delete; }{ Volume{"+str(i+2)+"}; Delete; };"
+        Lattice.write(Union1+"\n"+Intersection1+"\n")   
+        Lattice.write("Box("+str(len(A)+4)+") = {0,0,0,"+str(sb)+","+str(sb)+","+str(sb)+"};"+ "\n" )       
+        
+        for i in range (len(A)):
+            a=A[i][0][0]
+            b=A[i][0][1]
+            c=A[i][0][2]
+            d=A[i][1][0]-A[i][0][0]
+            e=A[i][1][1]-A[i][0][1]
+            f=A[i][1][2]-A[i][0][2]
+            r=A[i][3][0]
+            Lattice.write("Cylinder("+str(len(A)+i+5)+") = {"+ str(a)+","+ str(b)+","+ str(c)+","+ str(d)+","+ str(e)+","+ str(f)+","+ str(r)+",2*Pi};"+ "\n")
+            if i in range (0,8):	                 	
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(len(A)+i+5)+"}; }"+ "\n")          
+            if i != len(A)-1:
+                Union2+="Volume{"+str(len(A)+i+5)+"};"
+            else:
+                Union2+="Delete; }{ Volume{"+str(len(A)+i+5)+"}; Delete; };"
+        
+        Lattice.write(Union2+"\n"+Intersection2+"\n")
+        
+        
+        Lattice.write(Difference)
+     
+        #----------------------------------------------
+        
+        #lattice creation
+    
+        number = 2*len(A)+7
+        
+        Union = "BooleanUnion{"
+        
+        for i in range (length):
+            for j in range (width):
+                for h in range (height):
+                    if i==0 and j==0 and h==0:
+                        continue
+                    else:
+                        number+=1
+                        Union = Union + "Volume{"+str(number)+"};"
+                        Lattice.write("\n"+"Translate {"+str(i*sb)+","+str(j*sb) +"," + str(h*sb)+"} {Duplicata { Volume{"+str(2*len(A)+7)+"}; }}")
+            print(i)
+        if length!=1 or width!=1 or height!=1:
+            Lattice.write("\n"+Union+ "Delete;}{Volume{"+str(2*len(A)+7)+"}; Delete;}")
+    
+    
+    ### for geometry 8    
+    elif name =="8":
+        Union = "BooleanUnion("+str(len(A)+2)+") = {"
+        Intersection ="BooleanIntersection("+str(len(A)+3)+")={ Volume{"+str(len(A)+2)+"}; Delete;}{Volume{1}; Delete;};"
+        Lattice = open (name+"_"+rho+".geo",'w')
+        Lattice.write( "SetFactory('OpenCASCADE');" + "\n" )
+        for i in range (1,25):
+            Lattice.write("Disk("+str(i)+") = {"+str(sb/2)+","+str(sb/2)+","+str(sb/2)+", "+str(sb)+", "+str(sb)+"};"+ "\n")
+            if i in range(1,9):
+                Lattice.write("Rotate {{1, 0, 0}, {"+str(sb/2)+","+str(sb/2)+","+str(sb/2)+"}, Pi/2} { Surface{"+str(i)+"}; }"+ "\n")
+            if i in range (9,17):
+                Lattice.write("Rotate {{0, 1, 0}, {"+str(sb/2)+","+str(sb/2)+","+str(sb/2)+"}, Pi/2} { Surface{"+str(i)+"}; }"+ "\n")
+        for i in range (len(B)):
+            a=B[i][0][0]
+            b=B[i][0][1]
+            c=B[i][0][2]
+            d=B[i][1][0]-B[i][0][0]
+            e=B[i][1][1]-B[i][0][1]
+            f=B[i][1][2]-B[i][0][2]
+            r=r0
+            Lattice.write("Cylinder("+str(i+25)+") = {"+ str(a)+","+ str(b)+","+ str(c)+","+ str(d)+","+ str(e)+","+ str(f)+","+ str(r)+",2*Pi};"+ "\n")
+        counter=1000
+        for i in range (1,25):
+            Lattice.write("BooleanFragments("+str(counter)+")={ Surface{"+str(i)+"}; Delete;}{Volume{"+str(i+48)+"}; Delete;};"+ "\n")
+            counter=counter+1
+        
+        Union = "BooleanUnion("+str(counter)+") = {"
+        Box = "Box("+(str(counter+1))+") = {0,0,0,"+str(x_cell)+","+str(y_cell)+","+str(z_cell)+"};"+ "\n" 
+        Intersection ="BooleanIntersection("+str(counter+2)+")={ Volume{"+str(counter)+"}; Delete;}{Volume{"+(str(counter+1))+"}; Delete;};" 
+        for i in range (25, counter):
+            if (i < (49)) or (i>=1000 and i < (counter-1)) :
+                Union+="Volume{"+str(i)+"};"
+            elif (i == counter-1):
+                Union+="Delete; }{ Volume{"+str(counter-1)+"}; Delete; };"
+    
+        Lattice.write(Union+"\n"+Box+"\n"+Intersection)   
+            
+    else:
+        Union = "BooleanUnion("+str(len(A)+2)+") = {"
+        Intersection ="BooleanIntersection("+str(len(A)+3)+")={ Volume{"+str(len(A)+2)+"}; Delete;}{Volume{1}; Delete;};"
+        
+    # Writing .geo file 
+    
+        Lattice = open (name+"_"+rho+".geo",'w')
+        Lattice.write( "SetFactory('OpenCASCADE');" + "\n" 
+                        + "Box(1) = {0,0,0,"+str(x_cell)+","+str(y_cell)+","+str(z_cell)+"};"+ "\n" )    
+    ##for the other geometries 
+            
+        for i in range (len(A)):
+            a=A[i][0][0]
+            b=A[i][0][1]
+            c=A[i][0][2]
+            d=A[i][1][0]-A[i][0][0]
+            e=A[i][1][1]-A[i][0][1]
+            f=A[i][1][2]-A[i][0][2]
+            r=r0
+            Lattice.write("Cylinder("+str(i+2)+") = {"+ str(a)+","+ str(b)+","+ str(c)+","+ str(d)+","+ str(e)+","+ str(f)+","+ str(r)+",2*Pi};"+ "\n")
+            
+            #cylinder rotations
+            if name =="2":
+            	if i in range (0,5):
+            	    Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, -Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")
+            	elif i in range (5,10):
+            	    Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")
+    	        elif i in range (10,18):
+    	            Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")
+            if name =="3" and (i in range (0,4)): 
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")
+            if name =="4" and (i in range (0,12)):
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")
+            if name =="5" and (i in range (0,18)): 
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")
+            if name =="7":
+            	if i in range (0,10): 	
+            	    Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")
+            	elif i in range (10,14):
+            	    Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")
+    	        elif i in range (14,22):
+    	            Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, -Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")  
+            if (name =="9" or name=="91") and (i in range (0,8)):	                 	
+                Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi} {Volume{"+str(i+2)+"}; }"+ "\n")  
+            if name=="91" and i in range (8,10):
+            	if i==8:
+    	            Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")              
+            	if i==9:  
+    	            Lattice.write("Rotate {{"+ str(d)+","+ str(e)+","+ str(f)+"}, {"+ str(a)+","+ str(b)+","+ str(c)+"}, -Pi/2} {Volume{"+str(i+2)+"}; }"+ "\n")
+    	                          
+            if i != len(A)-1:
+                Union+="Volume{"+str(i+2)+"};"
+            else:
+                Union+="Delete; }{ Volume{"+str(i+2)+"}; Delete; };"
+        
+        Lattice.write(Union+"\n"+Intersection)    
+        
+        #lattice creation    
+        number = len(A)+3
+        
+        Union = "BooleanUnion{"
+        
+        for i in range (length):
+            for j in range (width):
+                for h in range (height):
+                    if i==0 and j==0 and h==0:
+                        continue
+                    else:
+                        number+=1
+                        Union = Union + "Volume{"+str(number)+"};"
+                        Lattice.write("\n"+"Translate {"+str(i*x_cell)+","+str(j*y_cell) +"," + str(h*z_cell)+"} {Duplicata { Volume{"+str(len(A)+3)+"}; }}")
+            print(i)
+        if length!=1 or width!=1 or height!=1:
+            Lattice.write("\n"+Union+ "Delete;}{Volume{"+str(len(A)+3)+"}; Delete;}")
+    
+    # periodic mesh conditions
+    Rx= x_cell*length
+    Ry= y_cell*width
+    Rz= z_cell*height
+    
+    
+        
+    Lattice.write ("\n"+"e=1e-6;"
+                    +"\n" + "Rx="+ str(Rx) + ";"
+                    +"\n" + "Ry="+ str(Ry) + ";"
+                    +"\n" + "Rz="+ str(Rz) + ";"
+                    +"\n" + "bottomx() = Surface In BoundingBox{-e,-e,-e, e,"+str(Ry)+"+e,"+str(Rz)+"+e};"
+                    +"\n" + "bottomy() = Surface In BoundingBox{-e,-e,-e, "+str(Rx)+"+e,e,"+str(Rz)+"+e};"
+                    +"\n" + "bottomz() = Surface In BoundingBox{-e,-e,-e, "+str(Rx)+"+e,"+str(Ry)+"+e,e};"
+                    +"\n" + "For i In {0:#bottomx()-1}"
+                    +"\n" + "  bb() = BoundingBox Surface { bottomx(i) };"
+                    +"\n" + "  topx() = Surface In BoundingBox{bb(0)-e+"+str(Rx)+", bb(1)-e, bb(2)-e, bb(3)+e+"+str(Rx)+", bb(4)+e, bb(5)+e };"
+                    +"\n" + "  For j In {0:#topx()-1}"
+                    +"\n" + "    bb2() = BoundingBox Surface { topx(j) };"
+                    +"\n" + "    bb2(0) -= "+str(Rx)+";"
+                    +"\n" + "    bb2(3) -= "+str(Rx)+";"
+                    +"\n" + "    If(Fabs(bb2(0)-bb(0)) < e && Fabs(bb2(1)-bb(1)) < e &&"
+                    +"\n" + "       Fabs(bb2(2)-bb(2)) < e && Fabs(bb2(3)-bb(3)) < e &&"
+                    +"\n" + "       Fabs(bb2(4)-bb(4)) < e && Fabs(bb2(5)-bb(5)) < e)"
+                    +"\n" + "      Periodic Surface {topx(j)} = {bottomx(i)} Translate {"+str(Rx)+",0,0};"
+                    +"\n" + "    EndIf"
+                    +"\n" + "   EndFor"
+                    +"\n" + " EndFor"
+                    +"\n" + "For i In {0:#bottomy()-1}"
+                    +"\n" + "  cc() = BoundingBox Surface { bottomy(i) };"
+                    +"\n" + "  topy() = Surface In BoundingBox{cc(0)-e, cc(1)-e+"+str(Ry)+", cc(2)-e, cc(3)+e, cc(4)+e+"+str(Ry)+", cc(5)+e};"
+                    +"\n" + "  For j In {0:#topy()-1}"
+                    +"\n" + "    cc2() = BoundingBox Surface { topy(j) };"
+                    +"\n" + "    cc2(1) -= "+str(Ry)+";"
+                    +"\n" + "    cc2(4) -= "+str(Ry)+";"
+                    +"\n" + "    If(Fabs(cc2(0)-cc(0)) < e && Fabs(cc2(1)-cc(1)) < e &&"
+                    +"\n" + "       Fabs(cc2(2)-cc(2)) < e && Fabs(cc2(3)-cc(3)) < e &&"
+                    +"\n" + "       Fabs(cc2(4)-cc(4)) < e && Fabs(cc2(5)-cc(5)) < e)"
+                    +"\n" + "      Periodic Surface {topy(j)} = {bottomy(i)} Translate {0,"+str(Ry)+",0};"
+                    +"\n" + "    EndIf"
+                    +"\n" + "  EndFor"
+                    +"\n" + "EndFor"
+                    +"\n" + "For i In {0:#bottomz()-1}"
+                    +"\n" + "  dd() = BoundingBox Surface { bottomz(i) };"
+                    +"\n" + "  topz() = Surface In BoundingBox{ dd(0)-e, dd(1)-e, dd(2)-e+"+str(Rz)+", dd(3)+e, dd(4)+e, dd(5)+e+"+str(Rz)+" };"
+                    +"\n" + "  For j In {0:#topz()-1}"
+                    +"\n" + "    dd2() = BoundingBox Surface { topz(j) };"
+                    +"\n" + "    dd2(2) -= "+str(Rz)+";"
+                    +"\n" + "    dd2(5) -= "+str(Rz)+";"
+                    +"\n" + "    If(Fabs(dd2(0)-dd(0)) < e && Fabs(dd2(1)-dd(1)) < e &&"
+                    +"\n" + "       Fabs(dd2(2)-dd(2)) < e && Fabs(dd2(3)-dd(3)) < e &&"
+                    +"\n" + "       Fabs(dd2(4)-dd(4)) < e && Fabs(dd2(5)-dd(5)) < e)"
+                    +"\n" + "      Periodic Surface {topz(j)} = {bottomz(i)} Translate {0,0,"+str(Rz)+"};"
+                    +"\n" + "    EndIf"
+                    +"\n" + "  EndFor"
+                    +"\n" + "EndFor")
+    Lattice.close()
+

Lattice_Cell_script/vmf_bisection.f95

+subroutine vox_mesh_fort(vmf,A, nx, ny, nz,x_cell,y_cell,z_cell,vfraction,lA)
+
+implicit none
+
+integer, intent(in) :: nx
+integer, intent(in) :: ny
+integer, intent(in) :: nz
+integer, intent(in) :: lA
+double precision, intent(in) :: vfraction
+double precision, intent(in) :: x_cell
+double precision, intent(in) :: y_cell
+double precision, intent(in) :: z_cell
+double precision, intent(in), dimension(lA,3,3) :: A
+double precision, intent(out), dimension(nx,ny,nz) :: vmf
+integer i
+integer j
+integer k
+double precision phys_value
+double precision phys_value_temp
+integer l
+double precision, dimension (3) :: X
+double precision, dimension (3) :: XA
+double precision, dimension (3) :: XB
+double precision, dimension (3) :: XC
+double precision alpha
+double precision d
+double precision lim_int
+double precision lim_ext
+double precision density
+double precision material
+double precision rlow
+double precision rhigh
+double precision radius
+double precision tol
+double precision resto
+integer whole
+material=0
+whole=0
+rlow=0
+rhigh=0.5
+tol=1.E-10
+resto=1
+radius=0
+do while (resto .ge. tol)
+material=0
+whole=0
+radius=((rlow+rhigh)/2.0)
+!print *, "radius1=",radius
+!print *, "density1=",density
+iloop: do i=1,nx,1
+jloop: do j=1,ny,1
+kloop: do k=1,nz,1
+
+X(1)=(i-1)*x_cell/nx+x_cell/(2*nx)
+X(2)=(j-1)*y_cell/ny+y_cell/(2*ny)
+X(3)=(k-1)*z_cell/nz+z_cell/(2*nz)
+phys_value=0
+l=1
+do while (phys_value .ne. 1 .and. l .le. lA)
+lim_int= radius-x_cell/(2*nx)
+lim_ext= radius+x_cell/(2*nx)
+XA(1)=X(1)
+XA(2)=X(2)
+XA(3)=X(3)
+XB(1)=A(l,1,1)
+XB(2)=A(l,1,2)
+XB(3)=A(l,1,3)
+XC(1)=A(l,2,1)
+XC(2)=A(l,2,2)
+XC(3)=A(l,2,3)
+alpha = -((XC(1)-XB(1))*(XB(1)-XA(1))+(XC(2)-XB(2))*(XB(2)-XA(2))+(XC(3)-XB(3))*(XB(3)-XA(3))) &
+/((XC(1)-XB(1))**2+(XC(2)-XB(2))**2+(XC(3)-XB(3))**2)
+if (alpha .ge. 0 .and. alpha .le. 1) then
+d=sqrt((XB(1)-XA(1)+alpha*(XC(1)-XB(1)))**2+(XB(2)-XA(2)+alpha &
+*(XC(2)-XB(2)))**2+(XB(3)-XA(3)+alpha*(XC(3)-XB(3)))**2)
+elseif (alpha .lt. 0) then
+d=sqrt((XB(1)-XA(1))**2+(XB(2)-XA(2))**2+(XB(3)-XA(3))**2)
+else
+d=sqrt((XC(1)-XA(1))**2+(XC(2)-XA(2))**2+(XC(3)-XA(3))**2)
+endif
+if (d .le. lim_int) then 
+phys_value =1
+elseif (d .gt. lim_int .and. d .le. lim_ext) then
+phys_value_temp = 1/(lim_int-lim_ext)*d-lim_ext/(lim_int-lim_ext)
+if (phys_value_temp .gt. phys_value) then
+phys_value=phys_value_temp
+endif
+endif
+l=l+1
+end do
+vmf(i,j,k)=phys_value
+material=material+phys_value
+whole=whole+1
+end do kloop
+end do jloop
+end do iloop
+density=material/whole
+resto= abs(vfraction-density)
+if (density .gt. vfraction) then
+rhigh= radius
+elseif (density .lt. vfraction) then
+rlow= radius
+endif
+!print *, "radius2=",rhigh
+!print *, "radius1=",rlow
+print *, "vf_calc=",density
+end do
+!print *, "radius2=",rhigh
+!print *, "radius1=",rlow
+print *, "vf_input=",vfraction
+print *, "vf_final=",density
+print *, "RADIUSf=",radius
+open(99, file = 'Radius.dat', status = 'old')
+write(99,*) radius
+ close(99)
+end subroutine vox_mesh_fort

Lattice_Cell_script/vmf_hollow_ID.f95

+subroutine vox_mesh_fort2(vmf,A, nx, ny, nz,sb,lA)
+
+implicit none
+
+integer, intent(in) :: nx
+integer, intent(in) :: ny
+integer, intent(in) :: nz
+integer, intent(in) :: lA
+double precision, intent(in) :: sb
+double precision, intent(in), dimension(lA,4,3) :: A
+double precision, intent(out), dimension(nx,ny,nz) :: vmf
+integer i
+integer j
+integer k
+double precision phys_value
+double precision phys_value_temp
+integer l
+double precision, dimension (3) :: X
+double precision, dimension (3) :: XA
+double precision, dimension (3) :: XB
+double precision, dimension (3) :: XC
+double precision alpha
+double precision d
+double precision lim_int
+double precision lim_ext
+double precision density
+double precision material
+integer whole
+material=0
+whole=0
+iloop: do i=1,nx,1
+jloop: do j=1,ny,1
+kloop: do k=1,nz,1
+print * ,i,j,k
+X(1)=(i-1)*sb/nx+sb/(2*nx)
+X(2)=(j-1)*sb/ny+sb/(2*ny)
+X(3)=(k-1)*sb/nz+sb/(2*nz)
+phys_value=0
+l=1
+do while (phys_value .ne. 1 .and. l .le. lA)
+lim_int= A(l,4,1)-sb/(2*nx)
+lim_ext= A(l,4,1)+sb/(2*nx)
+XA(1)=X(1)
+XA(2)=X(2)
+XA(3)=X(3)
+XB(1)=A(l,1,1)
+XB(2)=A(l,1,2)
+XB(3)=A(l,1,3)
+XC(1)=A(l,2,1)
+XC(2)=A(l,2,2)
+XC(3)=A(l,2,3)
+alpha = -((XC(1)-XB(1))*(XB(1)-XA(1))+(XC(2)-XB(2))*(XB(2)-XA(2))+(XC(3)-XB(3))*(XB(3)-XA(3))) &
+/((XC(1)-XB(1))**2+(XC(2)-XB(2))**2+(XC(3)-XB(3))**2)
+if (alpha .ge. 0 .and. alpha .le. 1) then
+d=sqrt((XB(1)-XA(1)+alpha*(XC(1)-XB(1)))**2+(XB(2)-XA(2)+alpha &
+*(XC(2)-XB(2)))**2+(XB(3)-XA(3)+alpha*(XC(3)-XB(3)))**2)
+elseif (alpha .lt. 0) then
+d=sqrt((XB(1)-XA(1))**2+(XB(2)-XA(2))**2+(XB(3)-XA(3))**2)
+else
+d=sqrt((XC(1)-XA(1))**2+(XC(2)-XA(2))**2+(XC(3)-XA(3))**2)
+endif
+if (d .le. lim_int) then 
+phys_value =1
+elseif (d .gt. lim_int .and. d .le. lim_ext) then
+phys_value_temp = 1/(lim_int-lim_ext)*d-lim_ext/(lim_int-lim_ext)
+if (phys_value_temp .gt. phys_value) then
+phys_value=phys_value_temp
+endif
+endif
+l=l+1
+end do 
+vmf(i,j,k)=phys_value
+material=material+phys_value
+whole=whole+1
+end do kloop
+end do jloop
+end do iloop
+density=material/whole
+print *, "density=",density
+open(99, file = 'Density.dat', status = 'old')
+write(99,*) density
+ close(99)
+
+end subroutine vox_mesh_fort2

Lattice_Cell_script/vmf_hollow_OD.f95

+subroutine vox_mesh_fort(vmf,A, nx, ny, nz,sb,lA)
+
+implicit none
+
+integer, intent(in) :: nx
+integer, intent(in) :: ny
+integer, intent(in) :: nz
+integer, intent(in) :: lA
+double precision, intent(in) :: sb
+double precision, intent(in), dimension(lA,4,3) :: A
+double precision, intent(out), dimension(nx,ny,nz) :: vmf
+integer i
+integer j
+integer k
+double precision phys_value
+double precision phys_value_temp
+integer l
+double precision, dimension (3) :: X
+double precision, dimension (3) :: XA
+double precision, dimension (3) :: XB
+double precision, dimension (3) :: XC
+double precision alpha
+double precision d
+double precision lim_int
+double precision lim_ext
+double precision density
+double precision material
+integer whole
+material=0
+whole=0
+iloop: do i=1,nx,1
+jloop: do j=1,ny,1
+kloop: do k=1,nz,1
+print * ,i,j,k
+X(1)=(i-1)*sb/nx+sb/(2*nx)
+X(2)=(j-1)*sb/ny+sb/(2*ny)
+X(3)=(k-1)*sb/nz+sb/(2*nz)
+phys_value=0
+l=1
+do while (phys_value .ne. 1 .and. l .le. lA)
+lim_int= A(l,3,1)-sb/(2*nx)
+lim_ext= A(l,3,1)+sb/(2*nx)
+XA(1)=X(1)
+XA(2)=X(2)
+XA(3)=X(3)
+XB(1)=A(l,1,1)
+XB(2)=A(l,1,2)
+XB(3)=A(l,1,3)
+XC(1)=A(l,2,1)
+XC(2)=A(l,2,2)
+XC(3)=A(l,2,3)
+alpha = -((XC(1)-XB(1))*(XB(1)-XA(1))+(XC(2)-XB(2))*(XB(2)-XA(2))+(XC(3)-XB(3))*(XB(3)-XA(3))) &
+/((XC(1)-XB(1))**2+(XC(2)-XB(2))**2+(XC(3)-XB(3))**2)
+if (alpha .ge. 0 .and. alpha .le. 1) then
+d=sqrt((XB(1)-XA(1)+alpha*(XC(1)-XB(1)))**2+(XB(2)-XA(2)+alpha &
+*(XC(2)-XB(2)))**2+(XB(3)-XA(3)+alpha*(XC(3)-XB(3)))**2)
+elseif (alpha .lt. 0) then
+d=sqrt((XB(1)-XA(1))**2+(XB(2)-XA(2))**2+(XB(3)-XA(3))**2)
+else
+d=sqrt((XC(1)-XA(1))**2+(XC(2)-XA(2))**2+(XC(3)-XA(3))**2)
+endif
+if (d .le. lim_int) then 
+phys_value =1
+elseif (d .gt. lim_int .and. d .le. lim_ext) then
+phys_value_temp = 1/(lim_int-lim_ext)*d-lim_ext/(lim_int-lim_ext)
+if (phys_value_temp .gt. phys_value) then
+phys_value=phys_value_temp
+endif
+endif
+l=l+1
+end do 
+vmf(i,j,k)=phys_value
+material=material+phys_value
+whole=whole+1
+end do kloop
+end do jloop
+end do iloop
+density=material/whole
+print *, "density=",density
+open(99, file = 'Density.dat', status = 'old')
+write(99,*) density
+ close(99)
+end subroutine vox_mesh_fort

Lattice_Cell_script/vmf_test.f95

+subroutine vox_mesh_fort(vmf,A, nx, ny, nz,x_cell,y_cell,z_cell,lA)
+
+implicit none
+
+integer, intent(in) :: nx
+integer, intent(in) :: ny
+integer, intent(in) :: nz
+integer, intent(in) :: lA
+double precision, intent(in) :: x_cell
+double precision, intent(in) :: y_cell
+double precision, intent(in) :: z_cell
+double precision, intent(in), dimension(lA,3,3) :: A
+double precision, intent(out), dimension(nx,ny,nz) :: vmf
+integer i
+integer j
+integer k
+double precision phys_value
+double precision phys_value_temp
+integer l
+double precision, dimension (3) :: X
+double precision, dimension (3) :: XA
+double precision, dimension (3) :: XB
+double precision, dimension (3) :: XC
+double precision alpha
+double precision d
+double precision lim_int
+double precision lim_ext
+double precision density
+double precision material
+double precision radius
+integer whole
+radius= A(1,3,1)
+material=0
+whole=0
+iloop: do i=1,nx,1
+jloop: do j=1,ny,1
+kloop: do k=1,nz,1
+print * ,i,j,k
+X(1)=(i-1)*x_cell/nx+x_cell/(2*nx)
+X(2)=(j-1)*y_cell/ny+y_cell/(2*ny)
+X(3)=(k-1)*z_cell/nz+z_cell/(2*nz)
+phys_value=0
+l=1
+do while (phys_value .ne. 1 .and. l .le. lA)
+lim_int= A(l,3,1)-x_cell/(2*nx)
+lim_ext= A(l,3,1)+x_cell/(2*nx)
+XA(1)=X(1)
+XA(2)=X(2)
+XA(3)=X(3)
+XB(1)=A(l,1,1)
+XB(2)=A(l,1,2)
+XB(3)=A(l,1,3)
+XC(1)=A(l,2,1)
+XC(2)=A(l,2,2)
+XC(3)=A(l,2,3)
+alpha = -((XC(1)-XB(1))*(XB(1)-XA(1))+(XC(2)-XB(2))*(XB(2)-XA(2))+(XC(3)-XB(3))*(XB(3)-XA(3))) &
+/((XC(1)-XB(1))**2+(XC(2)-XB(2))**2+(XC(3)-XB(3))**2)
+if (alpha .ge. 0 .and. alpha .le. 1) then
+d=sqrt((XB(1)-XA(1)+alpha*(XC(1)-XB(1)))**2+(XB(2)-XA(2)+alpha &
+*(XC(2)-XB(2)))**2+(XB(3)-XA(3)+alpha*(XC(3)-XB(3)))**2)
+elseif (alpha .lt. 0) then
+d=sqrt((XB(1)-XA(1))**2+(XB(2)-XA(2))**2+(XB(3)-XA(3))**2)
+else
+d=sqrt((XC(1)-XA(1))**2+(XC(2)-XA(2))**2+(XC(3)-XA(3))**2)
+endif
+if (d .le. lim_int) then 
+phys_value =1
+elseif (d .gt. lim_int .and. d .le. lim_ext) then
+phys_value_temp = 1/(lim_int-lim_ext)*d-lim_ext/(lim_int-lim_ext)
+if (phys_value_temp .gt. phys_value) then
+phys_value=phys_value_temp
+endif
+endif
+l=l+1
+end do
+!print *, "lim_int=" ,lim_int
+!print *, "lim_ext=", lim_ext 
+vmf(i,j,k)=phys_value
+material=material+phys_value
+whole=whole+1
+end do kloop
+end do jloop
+end do iloop
+density=material/whole
+!print *, "lim_int=" ,lim_int
+!print *, "lim_ext=", lim_ext
+print *, "densityf=",density
+print *, "density=",density
+open(99, file = 'Density.dat', status = 'old')
+write(99,*) density
+ close(99)
+print *, "radiusf=", radius
+open(99, file = 'Radius.dat', status = 'old')
+write(99,*) radius
+ close(99)
+end subroutine vox_mesh_fort