diff --git a/Lattice_Cell_script/Auto_geometry_parametric.py b/Lattice_Cell_script/Auto_geometry_parametric.py
new file mode 100755
index 0000000000000000000000000000000000000000..ebc8723232dc9bb60eb46498a312ab5eb8bb3466
--- /dev/null
+++ b/Lattice_Cell_script/Auto_geometry_parametric.py
@@ -0,0 +1,593 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 14 10:39:21 2020
+
+@author: Lucia
+"""
+import numpy as np
+import math
+
+import vmf_test
+import vmf_hollow_OD
+import vmf_hollow_ID
+import vmf_bisection
+from utils_geo import *
+
+# f2py -c name.f95 -m name
+
+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
+file = open('Radius.dat','a+')
+file.close()
+file = open('Density.dat','a+')
+file.close()
+
+
+def Geo_generator(name, rho, vfr, length, width, height, sb, Radius, vfraction):
+ if vfr=="a":
+ if name =="92":
+ print ("NOTE: geo 92 does not work with option a")
+ sys.exit()
+ print ("you chose to introduce volume fraction")
+
+ if vfr=="b":
+ print ("you chose to introduce radius")
+ '''if(name=="2"):
+ NumVolume = 29
+ if(length+width+height>3):
+ NumVolume = 1
+ elif(name=="3"):
+ NumVolume = 11
+ elif(name=="9"):
+ NumVolume = 19
+ elif(name=="7"):
+ NumVolume = 43'''
+
+ A=[]
+ B=[]
+ x_cell=0
+ y_cell=0
+ z_cell=0
+ r=0
+ if name =="2":
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ r=Radius
+ A=A+[[[sb/2,0,sb/2],[0,sb/2,0],[r,0,0]], #0
+ [[sb/2,0,sb/2],[sb,sb/2,sb],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[0,1.5*sb,0],[r,0,0]], #15
+ [[sb/2,sb,sb/2],[sb,1.5*sb,sb],[r,0,0]], #16
+ [[sb,sb/2,0],[0,sb/2,sb],[r,0,0]], #17
+
+ [[sb/2,0,sb/2],[sb,-sb/2,0],[r,0,0]], #2
+ [[sb/2,0,sb/2],[0,-sb/2,sb],[r,0,0]], #3
+ [[0,sb/2,0],[sb,sb/2,sb],[r,0,0]], #4
+ [[sb/2,sb,sb/2],[sb,sb/2,0],[r,0,0]], #13
+ [[sb/2,sb,sb/2],[0,sb/2,sb],[r,0,0]], #14
+
+ [[0,sb/2,0],[-sb,sb/2,sb],[r,0,0]], #5
+ [[0,sb/2,0],[-sb/2,sb,sb/2],[r,0,0]], #6
+ [[sb,sb/2,sb],[0,sb/2,2*sb],[r,0,0]], #10
+ [[sb,sb/2,sb],[1.5*sb,sb,sb/2],[r,0,0]], #11
+ [[sb,sb/2,sb],[2*sb,sb/2,0],[r,0,0]], #12
+ [[sb,sb/2,0],[sb/2,0,-sb/2],[r,0,0]], #20
+ [[sb,sb/2,0],[0,sb/2,-sb],[r,0,0]], #21
+ [[0,sb/2,sb],[sb/2,0,1.5*sb],[r,0,0]], #22
+
+ [[0,sb/2,0],[sb,sb/2,-sb],[r,0,0]], #7
+ [[0,sb/2,0],[sb/2,sb,-sb/2],[r,0,0]], #8
+ [[sb,sb/2,sb],[sb/2,sb,1.5*sb],[r,0,0]], #9
+ [[sb,sb/2,0],[1.5*sb,0,sb/2],[r,0,0]], #18
+ [[sb,sb/2,0],[2*sb,sb/2,sb],[r,0,0]], #19
+ [[0,sb/2,sb],[sb,sb/2,2*sb],[r,0,0]], #23
+ [[0,sb/2,sb],[-sb/2,0,sb/2],[r,0,0]], #24
+ [[0,sb/2,sb],[-sb,sb/2,0],[r,0,0]]] #25
+
+ elif name =="3":
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ r=Radius
+ A=A+[[[0,0,0],[sb/2,sb/2,sb/2],[r,0,0]], #0
+ [[0,sb,0],[sb/2,sb/2,sb/2],[r,0,0]], #2
+ [[sb/2,sb/2,sb/2],[0,sb,sb],[r,0,0]], #5
+ [[sb/2,sb/2,sb/2],[0,0,sb],[r,0,0]], #7
+
+ [[sb,0,0],[sb/2,sb/2,sb/2],[r,0,0]], #1
+ [[sb,sb,0],[sb/2,sb/2,sb/2],[r,0,0]], #3
+ [[sb/2,sb/2,sb/2],[sb,sb,sb],[r,0,0]], #4
+ [[sb/2,sb/2,sb/2],[sb,0,sb],[r,0,0]]] #6
+
+
+ elif name== "4":
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ r=Radius
+
+ A=A+[[[sb/4,sb/4,sb/4],[0,0,0],[r,0,0]], #0
+ [[sb/4,sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #1
+ [[3.*sb/4,sb/4,sb/4],[sb,sb/2,sb/2],[r,0,0]], #7
+ [[3.*sb/4,3*sb/4,sb/4],[sb,sb/2,sb/2],[r,0,0]], #11
+ [[sb/4,3.*sb/4,sb/4],[0,sb,0],[r,0,0]], #12
+ [[sb/4,3.*sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #13
+ [[sb/4,sb/4,3.*sb/4],[0,sb/2,sb/2],[r,0,0]], #19
+ [[3*sb/4,sb/4,3.*sb/4],[sb,0,sb],[r,0,0]], #20
+ [[3*sb/4,sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #21
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb,sb,sb],[r,0,0]], #24
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #25
+ [[sb/4,3.*sb/4,3.*sb/4],[0,sb/2,sb/2],[r,0,0]], #31
+
+ [[sb/4,sb/4,sb/4],[sb/2,0,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/4],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3.*sb/4,sb/4,sb/4],[sb,0,0],[r,0,0]], #4
+ [[3.*sb/4,sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #5
+ [[3.*sb/4,sb/4,sb/4],[sb/2,0,sb/2],[r,0,0]], #6
+ [[3.*sb/4,3*sb/4,sb/4],[sb,sb,0],[r,0,0]], #8
+ [[3.*sb/4,3*sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #9
+ [[3.*sb/4,3*sb/4,sb/4],[sb/2,sb,sb/2],[r,0,0]], #10
+ [[sb/4,3.*sb/4,sb/4],[sb/2,sb,sb/2],[r,0,0]], #14
+ [[sb/4,3.*sb/4,sb/4],[0,sb/2,sb/2],[r,0,0]], #15
+ [[sb/4,sb/4,3.*sb/4],[0,0,sb],[r,0,0]], #16
+ [[sb/4,sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #17
+ [[sb/4,sb/4,3.*sb/4],[sb/2,0,sb/2],[r,0,0]], #18
+ [[3*sb/4,sb/4,3.*sb/4],[sb/2,0,sb/2],[r,0,0]], #22
+ [[3*sb/4,sb/4,3.*sb/4],[sb,sb/2,sb/2],[r,0,0]], #23
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb/2,sb,sb/2],[r,0,0]], #26
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb,sb/2,sb/2],[r,0,0]], #27
+ [[sb/4,3.*sb/4,3*sb/4],[0,sb,sb],[r,0,0]], #28
+ [[sb/4,3.*sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #29
+ [[sb/4,3.*sb/4,3.*sb/4],[sb/2,sb,sb/2],[r,0,0]]] #30
+
+ elif name == "5":
+
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ r=Radius
+ #angles: 56º and 45º
+ #inner star points:
+
+ s56_f=(sb/2)+((math.sin(math.radians(34))/math.sin(math.radians(56)))*sb/4)
+ s56_c=(sb/2)-((math.sin(math.radians(34))/math.sin(math.radians(56)))*sb/4)
+
+ A=[[[sb,3*sb/4,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #2
+ [[sb,3*sb/4,sb],[3*sb/4,sb/2,sb/2],[r,0,0]], #3
+ [[sb,sb/4,sb],[s56_f,sb/4,sb/2],[r,0,0]], #6
+ [[sb,sb/4,sb],[3*sb/4,sb/2,sb/2],[r,0,0]], #7
+ [[0,3*sb/4,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #8
+ [[0,3*sb/4,0],[sb/4,sb/2,sb/2],[r,0,0]], #9
+ [[0,sb/4,0],[s56_c,sb/4,sb/2],[r,0,0]], #12
+ [[0,sb/4,0],[sb/4,sb/2,sb/2],[r,0,0]], #13
+ [[sb/2,0,0],[s56_c,sb/4,sb/2],[r,0,0]], #16
+ [[sb/2,0,0],[s56_f,sb/4,sb/2],[r,0,0]], #17
+ [[sb/2,sb,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #22
+ [[sb/2,sb,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #23
+ [[sb/2,sb/2,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #24
+ [[sb/2,sb/2,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #26
+ [[sb/2,sb/2,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #28
+ [[sb/2,sb/2,sb],[s56_f,sb/4,sb/2],[r,0,0]], #31
+ [[sb/2,sb/2,sb],[s56_c,sb/4,sb/2],[r,0,0]], #33
+ [[sb/2,sb/2,sb],[sb/4,sb/2,sb/2],[r,0,0]], #35
+
+ [[sb,3*sb/4,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #0
+ [[sb,3*sb/4,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #1
+ [[sb,sb/4,0],[s56_f,sb/4,sb/2],[r,0,0]], #4
+ [[sb,sb/4,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #5
+ [[0,3*sb/4,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #10
+ [[0,3*sb/4,sb],[sb/4,sb/2,sb/2],[r,0,0]], #11
+ [[0,sb/4,sb],[s56_c,sb/4,sb/2],[r,0,0]], #14
+ [[0,sb/4,sb],[sb/4,sb/2,sb/2],[r,0,0]], #15
+ [[sb/2,0,sb],[s56_c,sb/4,sb/2],[r,0,0]], #18
+ [[sb/2,0,sb],[s56_f,sb/4,sb/2],[r,0,0]], #19
+ [[sb/2,sb,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #20
+ [[sb/2,sb,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #21
+ [[sb/2,sb/2,0],[s56_f,sb/4,sb/2],[r,0,0]], #25
+ [[sb/2,sb/2,0],[s56_c,sb/4,sb/2],[r,0,0]], #27
+ [[sb/2,sb/2,0],[sb/4,sb/2,sb/2],[r,0,0]], #29
+ [[sb/2,sb/2,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #30
+ [[sb/2,sb/2,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #32
+ [[sb/2,sb/2,sb],[3*sb/4,sb/2,sb/2],[r,0,0]]] #34
+
+ elif name == "6":
+
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ r=Radius
+ A=[[[0,0,0],[sb,0,sb],[r,0,0]],
+ [[0,0,sb],[sb,0,0],[r,0,0]],
+ [[sb,0,sb],[sb,sb,0],[r,0,0]],
+ [[sb,0,0],[sb,sb,sb],[r,0,0]],
+ [[sb,sb,sb],[0,sb,0],[r,0,0]],
+ [[0,sb,sb],[sb,sb,0],[r,0,0]],
+ [[0,sb,sb],[0,0,0],[r,0,0]],
+ [[0,0,sb],[0,sb,0],[r,0,0]],
+ [[0,0,0],[sb,sb,0],[r,0,0]],
+ [[sb,0,0],[0,sb,0],[r,0,0]],
+ [[0,0,sb],[sb,sb,sb],[r,0,0]],
+ [[sb,0,sb],[0,sb,sb],[r,0,0]],
+ [[sb/2,0,sb/2],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]],
+ [[sb/2,0,sb/2],[sb/2,sb/2,0],[r,0,0]],
+ [[sb/2,0,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[sb/2,sb/2,0],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb/2,0],[0,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb/2,0],[sb/2,sb,sb/2],[r,0,0]],
+ [[0,sb/2,sb/2],[sb/2,sb,sb/2],[r,0,0]],
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[0,sb/2,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[sb,sb/2,sb/2],[sb/2,sb/2,sb],[r,0,0]]]
+
+ elif name== "7":
+
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ #hexagon side lenght: l
+ #hexagon side to side lenght: sb
+ #lattice x lenght: l2
+
+ l=(sb/2)/(math.sin(math.radians(60)))
+ x_cell=3*l
+
+ if vfr=="b":
+ r=Radius
+
+ A=[[[l,sb,0],[2*l,sb,0],[r,0,0]], #0
+ [[l,0,0],[2*l,0,0],[r,0,0]], #3
+ [[2.5*l,sb/2,0],[3*l,sb/2,0],[r,0,0]], #7
+ [[2.5*l,sb/2,sb],[3*l,sb/2,sb],[r,0,0]], #23
+ [[l,sb,0],[l/2,sb/2,sb/2],[r,0,0]], #24
+ [[2.5*l,sb/2,0],[2*l,sb,sb/2],[r,0,0]], #29
+ [[l,0,sb],[l/2,sb/2,sb/2],[r,0,0]], #32
+ [[2.5*l,sb/2,sb],[2*l,0,sb/2],[r,0,0]], #34
+ [[l/2,sb/2,sb/2],[-l/2,sb/2,0],[r,0,0]], #38
+ [[l/2,sb/2,sb/2],[-l/2,sb/2,sb],[r,0,0]], #39
+
+ [[l,sb,0],[l/2,sb/2,0],[r,0,0]], #1
+ [[l/2,sb/2,0],[l,0,0],[r,0,0]], #2
+ [[2*l,0,0],[2.5*l,sb/2,0],[r,0,0]], #4
+ [[2.5*l,sb/2,0],[2*l,sb,0],[r,0,0]], #5
+
+ [[l,sb,sb/2],[l/2,sb/2,sb/2],[r,0,0]], #9
+ [[l/2,sb/2,sb/2],[l,0,sb/2],[r,0,0]], #10
+ [[2*l,0,sb/2],[2.5*l,sb/2,sb/2],[r,0,0]], #12
+ [[2.5*l,sb/2,sb/2],[2*l,sb,sb/2],[r,0,0]], #13
+ [[l,sb,sb],[l/2,sb/2,sb],[r,0,0]], #17
+ [[l/2,sb/2,sb],[l,0,sb],[r,0,0]], #18
+ [[2*l,0,sb],[2.5*l,sb/2,sb],[r,0,0]], #20
+ [[2.5*l,sb/2,sb],[2*l,sb,sb],[r,0,0]], #21
+
+
+
+ [[l/2,sb/2,0],[0,sb/2,0],[r,0,0]], #6
+ [[l,sb,sb/2],[2*l,sb,sb/2],[r,0,0]], #8
+ [[l,0,sb/2],[2*l,0,sb/2],[r,0,0]], #11
+ [[l/2,sb/2,sb/2],[0,sb/2,sb/2],[r,0,0]], #14
+ [[3*l,sb/2,sb/2],[2.5*l,sb/2,sb/2],[r,0,0]], #15 origin & end coordinates changed: meshing problems
+ [[l,sb,sb],[2*l,sb,sb],[r,0,0]], #16
+ [[l,0,sb],[2*l,0,sb],[r,0,0]], #19
+ [[l/2,sb/2,sb],[0,sb/2,sb],[r,0,0]], #22
+
+ [[l,sb,0],[2*l,sb,sb/2],[r,0,0]], #25
+ [[l,0,0],[l/2,sb/2,sb/2],[r,0,0]], #26
+ [[l,0,0],[2*l,0,sb/2],[r,0,0]], #27
+ [[2.5*l,sb/2,0],[2*l,0,sb/2],[r,0,0]], #28
+ [[l,sb,sb],[l/2,sb/2,sb/2],[r,0,0]], #30
+ [[l,sb,sb],[2*l,sb,sb/2],[r,0,0]], #31
+ [[l,0,sb],[2*l,0,sb/2],[r,0,0]], #33
+ [[2.5*l,sb/2,sb],[2*l,sb,sb/2],[r,0,0]], #35
+ [[2.5*l,sb/2,0],[3.5*l,sb/2,sb/2],[r,0,0]], #36
+ [[2.5*l,sb/2,sb],[3.5*l,sb/2,sb/2],[r,0,0]]] #37
+
+
+ elif name== "8":
+
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ r=Radius
+
+ #bow height: bow_h
+ #bow points: bow_p
+ #middle point:(connection) tie_l
+ #tie points:(connection) tie_p
+ #angle of bow corner: bow_angle
+ #f and c correspond to close and far from axis x,y,z
+
+ #inputs:
+ #bow_h=float(input())
+ bow_h=0.75
+ print ("bow heinght between 0 and size of unit cell, chosen:", bow_h)
+ #bow_angle=float(input())
+ bow_angle=75
+ print ("bow angle between 60º and 90º, chosen: "+ str(bow_angle)+"º")
+
+ #-------calculations
+
+ tie_l=((math.sin(math.radians(90-bow_angle))/math.sin(math.radians(bow_angle)))*sb/2)
+ bow_pc=(sb-bow_h)/2
+ bow_pf=bow_h+bow_pc
+ tie_pc=bow_pc+tie_l
+ tie_pf=bow_pf-tie_l
+
+ # only for B
+ m=(bow_pc-tie_pc)/((sb)-(sb/2))
+ pcortec=(sb/4)*m+bow_pc
+ pcortef=(-sb/4)*m+bow_pf
+ centrec=(-sb/10)*m+tie_pc
+ centref=(sb/10)*m+tie_pf
+ #-----------
+
+ # A for vox and phase
+ # B for mesh
+
+ A=[[[tie_pc,sb/2,0],[0,sb/2,0],[r,0,0]],
+ [[tie_pc,sb/2,0],[bow_pc,0,0],[r,0,0]],
+ [[tie_pc,sb/2,0],[bow_pc,sb,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[sb,sb/2,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[bow_pf,0,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[bow_pf,sb,0],[r,0,0]],
+ [[bow_pc,0,0],[bow_pf,0,0],[r,0,0]],
+ [[bow_pc,sb,0],[bow_pf,sb,0],[r,0,0]],
+ [[tie_pc,sb/2,sb],[0,sb/2,sb],[r,0,0]],
+ [[tie_pc,sb/2,sb],[bow_pc,0,sb],[r,0,0]],
+ [[tie_pc,sb/2,sb],[bow_pc,sb,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[sb,sb/2,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[bow_pf,0,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[bow_pc,0,sb],[bow_pf,0,sb],[r,0,0]],
+ [[bow_pc,sb,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb/2,0,0],[r,0,0]],
+ [[sb/2,0,tie_pc],[0,0,bow_pc],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb,0,bow_pc],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb/2,0,sb],[r,0,0]],
+ [[sb/2,0,tie_pf],[0,0,bow_pf],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb,0,bow_pf],[r,0,0]],
+ [[0,0,bow_pc],[0,0,bow_pf],[r,0,0]],
+ [[sb,0,bow_pc],[sb,0,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb/2,sb,0],[r,0,0]],
+ [[sb/2,sb,tie_pc],[0,sb,bow_pc],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb,sb,bow_pc],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb/2,sb,sb],[r,0,0]],
+ [[sb/2,sb,tie_pf],[0,sb,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb,sb,bow_pf],[r,0,0]],
+ [[0,sb,bow_pc],[0,sb,bow_pf],[r,0,0]],
+ [[sb,sb,bow_pc],[sb,sb,bow_pf],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,0,sb/2],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,bow_pc,0],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,bow_pc,sb],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,sb,sb/2],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,bow_pf,0],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,bow_pf,sb],[r,0,0]],
+ [[0,bow_pc,0],[0,bow_pf,0],[r,0,0]],
+ [[0,bow_pc,sb],[0,bow_pf,sb],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,0,sb/2],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,bow_pc,0],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,bow_pc,sb],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,sb,sb/2],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,bow_pf,0],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,bow_pf,sb],[r,0,0]],
+ [[sb,bow_pc,0],[sb,bow_pf,0],[r,0,0]],
+ [[sb,bow_pc,sb],[sb,bow_pf,sb],[r,0,0]]]
+
+
+#-----------
+ B=[[[tie_pc,sb/2,0],[0,sb/2,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[sb,sb/2,0],[r,0,0]],
+ [[bow_pc,0,0],[bow_pf,0,0],[r,0,0]],
+ [[bow_pc,sb,0],[bow_pf,sb,0],[r,0,0]],
+ [[tie_pc,sb/2,sb],[0,sb/2,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[sb,sb/2,sb],[r,0,0]],
+ [[bow_pc,0,sb],[bow_pf,0,sb],[r,0,0]],
+ [[bow_pc,sb,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb/2,0,0],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb/2,0,sb],[r,0,0]],
+ [[0,0,bow_pc],[0,0,bow_pf],[r,0,0]],
+ [[sb,0,bow_pc],[sb,0,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb/2,sb,0],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb/2,sb,sb],[r,0,0]],
+ [[0,sb,bow_pc],[0,sb,bow_pf],[r,0,0]],
+ [[sb,sb,bow_pc],[sb,sb,bow_pf],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,0,sb/2],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,sb,sb/2],[r,0,0]],
+ [[0,bow_pc,0],[0,bow_pf,0],[r,0,0]],
+ [[0,bow_pc,sb],[0,bow_pf,sb],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,0,sb/2],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,sb,sb/2],[r,0,0]],
+ [[sb,bow_pc,0],[sb,bow_pf,0],[r,0,0]],
+ [[sb,bow_pc,sb],[sb,bow_pf,sb],[r,0,0]],
+
+ [[centrec,sb/2+(sb/10),0],[pcortec,-sb/4,0],[r,0,0]],
+ [[centrec,sb/2-(sb/10),0],[pcortec,(sb+sb/4),0],[r,0,0]],
+ [[centref,sb/2+(sb/10),0],[pcortef,-sb/4,0],[r,0,0]],
+ [[centref,sb/2-(sb/10),0],[pcortef,(sb+sb/4),0],[r,0,0]],
+
+ [[centrec,sb/2+(sb/10),sb],[pcortec,-sb/4,sb],[r,0,0]],
+ [[centrec,sb/2-(sb/10),sb],[pcortec,(sb+sb/4),sb],[r,0,0]],
+ [[centref,sb/2+(sb/10),sb],[pcortef,-sb/4,sb],[r,0,0]],
+ [[centref,sb/2-(sb/10),sb],[pcortef,(sb+sb/4),sb],[r,0,0]],
+
+ [[sb/2+(sb/10),0,centrec],[-sb/4,0,pcortec],[r,0,0]],
+ [[sb/2-(sb/10),0,centrec],[(sb+sb/4),0,pcortec],[r,0,0]],
+ [[sb/2+(sb/10),0,centref],[-sb/4,0,pcortef],[r,0,0]],
+ [[sb/2-(sb/10),0,centref],[(sb+sb/4),0,pcortef],[r,0,0]],
+
+ [[sb/2+(sb/10),sb,centrec],[-sb/4,sb,pcortec],[r,0,0]],
+ [[sb/2-(sb/10),sb,centrec],[(sb+sb/4),sb,pcortec],[r,0,0]],
+ [[sb/2+(sb/10),sb,centref],[-sb/4,sb,pcortef],[r,0,0]],
+ [[sb/2-(sb/10),sb,centref],[(sb+sb/4),sb,pcortef],[r,0,0]],
+
+ [[0,centrec,sb/2+(sb/10)],[0,pcortec,-sb/4],[r,0,0]],
+ [[0,centrec,sb/2-(sb/10)],[0,pcortec,(sb+sb/4)],[r,0,0]],
+ [[0,centref,sb/2+(sb/10)],[0,pcortef,-sb/4],[r,0,0]],
+ [[0,centref,sb/2-(sb/10)],[0,pcortef,(sb+sb/4)],[r,0,0]],
+
+ [[sb,centrec,sb/2+(sb/10)],[sb,pcortec,-sb/4],[r,0,0]],
+ [[sb,centrec,sb/2-(sb/10)],[sb,pcortec,(sb+sb/4)],[r,0,0]],
+ [[sb,centref,sb/2+(sb/10)],[sb,pcortef,-sb/4],[r,0,0]],
+ [[sb,centref,sb/2-(sb/10)],[sb,pcortef,(sb+sb/4)],[r,0,0]]]
+
+
+ elif name== "9":
+
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ r=Radius
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0]], #13
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0]]] #15
+
+
+ elif name== "91":
+
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ r=Radius
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0]], #13
+ [[0,0,sb/2],[0,sb,sb/2],[r,0,0]], #16
+ [[sb,0,sb/2],[sb,sb,sb/2],[r,0,0]], #17
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0]]] #15
+
+
+ elif name=="92":
+ if vfr=="b":
+ r=Radius
+ ri=r/1.25
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0],[ri,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0],[ri,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #13
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0],[ri,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0],[ri,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]]] #15
+
+ r0=r
+
+
+ if vfr=="a":
+
+ voxelized(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B)
+ print ("7.a. Please introduce value of RADIUSf")
+ with open('Radius.dat','r') as data:
+ B_string = data.read()
+ r0=float(B_string)
+ mesh(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B, r0)
+
+ if vfr=="b":
+ mesh(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A,B,r0)
+ voxelized(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B)
+
+ GeoFile = name+"_"+rho+".geo"
+ cl_min = str(r0*0.3)
+ cl_max = str(r0*0.7)
+
+ HeadLine="""Mesh.CharacteristicLengthMin =""" + cl_min + """;
+Mesh.CharacteristicLengthMax =""" + cl_max + """;"""
+
+ EndLine = """
+Mytopx() = Surface In BoundingBox{Rx-e,-e,-e, Rx+e,Ry+e,Rz+e};
+Mytopy() = Surface In BoundingBox{-e,Ry-e,-e, Rx+e,Ry+e,Rz+e};
+Mytopz() = Surface In BoundingBox{-e,-e,Rz-e, Rx+e,Ry+e,Rz+e};
+
+AllVol() = Volume{:};
+
+Physical Surface(100) = {bottomx()};
+Physical Surface(110) = {bottomy()};
+Physical Surface(120) = {bottomz()};
+
+Physical Surface(101) = {Mytopx()};
+Physical Surface(111) = {Mytopy()};
+Physical Surface(121) = {Mytopz()};
+
+Physical Volume(51) ={AllVol()};"""
+#Physical Volume(51) ={"""+str(NumVolume)+""",1};"""
+
+ f = open(GeoFile, "r")
+ contents = f.readlines()
+ f.close()
+ contents.insert(1, HeadLine)
+ contents.append(EndLine)
+
+ f = open(GeoFile, "w")
+ contents = "".join(contents)
+ f.write(contents)
+ f.close()
+ os.system("rm vox_"+rho+".geo")
+ os.system("rm phase_map_"+rho+".geo")
+
+ return GeoFile, r0
+
+
diff --git a/Lattice_Cell_script/MicroSample.py b/Lattice_Cell_script/MicroSample.py
new file mode 100755
index 0000000000000000000000000000000000000000..fdacd210513e81e448b8cb9e0c50bf155925bd43
--- /dev/null
+++ b/Lattice_Cell_script/MicroSample.py
@@ -0,0 +1,59 @@
+#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")
diff --git a/Lattice_Cell_script/Readme.txt b/Lattice_Cell_script/Readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c600d1a71a471233a1df360814882d83b34dd1f6
--- /dev/null
+++ b/Lattice_Cell_script/Readme.txt
@@ -0,0 +1,6 @@
+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)
diff --git a/Lattice_Cell_script/compile.sh b/Lattice_Cell_script/compile.sh
new file mode 100755
index 0000000000000000000000000000000000000000..45c674a063f38966fd7fe85d30424a921fb13dd0
--- /dev/null
+++ b/Lattice_Cell_script/compile.sh
@@ -0,0 +1,5 @@
+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
+
diff --git a/Lattice_Cell_script/geometry_parametric.py b/Lattice_Cell_script/geometry_parametric.py
new file mode 100644
index 0000000000000000000000000000000000000000..1daa62bfd3eb089a8068cfb584dc015819d734bd
--- /dev/null
+++ b/Lattice_Cell_script/geometry_parametric.py
@@ -0,0 +1,672 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 14 10:39:21 2020
+
+@author: Lucia
+"""
+
+import numpy as np
+import math as m
+
+import vmf_test
+import vmf_hollow_OD
+import vmf_hollow_ID
+import vmf_bisection
+
+from utils_geo import *
+
+# f2py -c name.f95 -m name
+
+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
+
+from Auto_geometry_parametric import *
+
+file = open('Radius.dat','a+')
+file.close()
+file = open('Density.dat','a+')
+file.close()
+
+
+initial_data="""
+Geometry number
+--------------- -------
+2-Triangles 2
+Diagonals 3
+Rombo dodecahedron 4
+6-Star 5
+Octet 6
+Hexagon 7
+Auxetic 8
+Octahedron v1 9
+Octahedron v2 91
+Octahedron v3 (hollow) 92
+
+"""
+print (initial_data)
+print("1. Specify geometry number: ")
+name= input()
+print ("2. Number of voxels, 8,32,64,128...")
+#rho=input()
+rho=str(128)
+print (rho)
+
+print(name, rho)
+
+vf_radius="""
+3. What will you introduce a or b?:
+ a) volume fraction
+ b) radius
+Insert a or b
+"""
+print (vf_radius)
+vfr= input()
+if vfr=="a":
+ if name =="92":
+ print ("NOTE: geo 92 does not work with option a")
+ sys.exit()
+ print ("you chose to introduce volume fraction")
+
+if vfr=="b":
+ print ("you chose to introduce radius")
+
+#dimension of the lattice/ NUMBER OF UNIT CELL IN EACH DIRECTION
+print ("4. Number of unit cells in each direction")
+print ("length")
+length = int(input())
+print ("width")
+width = int(input())
+print("height")
+height = int(input())
+
+A=[]
+B=[]
+x_cell=0
+y_cell=0
+z_cell=0
+r=0
+print("5. Set size of unit cell")
+sb=float(input())
+if vfr=="b":
+ if name!="92":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+ elif name=="92":
+ print("HOLLOW geometry, inner radius dimension: 0.8 x outer radius")
+ print ("6.b. Specify outer radius of hollow struts, from 0 to "+ str(sb))
+ r=float(input())
+
+
+if vfr=="a":
+ if name!="8":
+ print( "6.a. Set volume fraction from 0 to 1")
+ vfraction=float(input())
+ if name=="8":
+ print( "6.a. Set volume fraction from 0 to 0.15")
+ vfraction=float(input())
+
+if vfr=="b":
+ vfraction=0.
+
+Radius=r
+VmFr=vfraction
+cell_name=name
+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)
+'''if name =="2":
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+ A=A+[[[sb/2,0,sb/2],[0,sb/2,0],[r,0,0]], #0
+ [[sb/2,0,sb/2],[sb,sb/2,sb],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[0,1.5*sb,0],[r,0,0]], #15
+ [[sb/2,sb,sb/2],[sb,1.5*sb,sb],[r,0,0]], #16
+ [[sb,sb/2,0],[0,sb/2,sb],[r,0,0]], #17
+
+ [[sb/2,0,sb/2],[sb,-sb/2,0],[r,0,0]], #2
+ [[sb/2,0,sb/2],[0,-sb/2,sb],[r,0,0]], #3
+ [[0,sb/2,0],[sb,sb/2,sb],[r,0,0]], #4
+ [[sb/2,sb,sb/2],[sb,sb/2,0],[r,0,0]], #13
+ [[sb/2,sb,sb/2],[0,sb/2,sb],[r,0,0]], #14
+
+ [[0,sb/2,0],[-sb,sb/2,sb],[r,0,0]], #5
+ [[0,sb/2,0],[-sb/2,sb,sb/2],[r,0,0]], #6
+ [[sb,sb/2,sb],[0,sb/2,2*sb],[r,0,0]], #10
+ [[sb,sb/2,sb],[1.5*sb,sb,sb/2],[r,0,0]], #11
+ [[sb,sb/2,sb],[2*sb,sb/2,0],[r,0,0]], #12
+ [[sb,sb/2,0],[sb/2,0,-sb/2],[r,0,0]], #20
+ [[sb,sb/2,0],[0,sb/2,-sb],[r,0,0]], #21
+ [[0,sb/2,sb],[sb/2,0,1.5*sb],[r,0,0]], #22
+
+ [[0,sb/2,0],[sb,sb/2,-sb],[r,0,0]], #7
+ [[0,sb/2,0],[sb/2,sb,-sb/2],[r,0,0]], #8
+ [[sb,sb/2,sb],[sb/2,sb,1.5*sb],[r,0,0]], #9
+ [[sb,sb/2,0],[1.5*sb,0,sb/2],[r,0,0]], #18
+ [[sb,sb/2,0],[2*sb,sb/2,sb],[r,0,0]], #19
+ [[0,sb/2,sb],[sb,sb/2,2*sb],[r,0,0]], #23
+ [[0,sb/2,sb],[-sb/2,0,sb/2],[r,0,0]], #24
+ [[0,sb/2,sb],[-sb,sb/2,0],[r,0,0]]] #25
+
+elif name =="3":
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+ A=A+[[[0,0,0],[sb/2,sb/2,sb/2],[r,0,0]], #0
+ [[0,sb,0],[sb/2,sb/2,sb/2],[r,0,0]], #2
+ [[sb/2,sb/2,sb/2],[0,sb,sb],[r,0,0]], #5
+ [[sb/2,sb/2,sb/2],[0,0,sb],[r,0,0]], #7
+
+ [[sb,0,0],[sb/2,sb/2,sb/2],[r,0,0]], #1
+ [[sb,sb,0],[sb/2,sb/2,sb/2],[r,0,0]], #3
+ [[sb/2,sb/2,sb/2],[sb,sb,sb],[r,0,0]], #4
+ [[sb/2,sb/2,sb/2],[sb,0,sb],[r,0,0]]] #6
+
+
+elif name== "4":
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to"+ str(sb))
+ r=float(input())
+
+ A=A+[[[sb/4,sb/4,sb/4],[0,0,0],[r,0,0]], #0
+ [[sb/4,sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #1
+ [[3.*sb/4,sb/4,sb/4],[sb,sb/2,sb/2],[r,0,0]], #7
+ [[3.*sb/4,3*sb/4,sb/4],[sb,sb/2,sb/2],[r,0,0]], #11
+ [[sb/4,3.*sb/4,sb/4],[0,sb,0],[r,0,0]], #12
+ [[sb/4,3.*sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #13
+ [[sb/4,sb/4,3.*sb/4],[0,sb/2,sb/2],[r,0,0]], #19
+ [[3*sb/4,sb/4,3.*sb/4],[sb,0,sb],[r,0,0]], #20
+ [[3*sb/4,sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #21
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb,sb,sb],[r,0,0]], #24
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #25
+ [[sb/4,3.*sb/4,3.*sb/4],[0,sb/2,sb/2],[r,0,0]], #31
+
+ [[sb/4,sb/4,sb/4],[sb/2,0,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/4],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3.*sb/4,sb/4,sb/4],[sb,0,0],[r,0,0]], #4
+ [[3.*sb/4,sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #5
+ [[3.*sb/4,sb/4,sb/4],[sb/2,0,sb/2],[r,0,0]], #6
+ [[3.*sb/4,3*sb/4,sb/4],[sb,sb,0],[r,0,0]], #8
+ [[3.*sb/4,3*sb/4,sb/4],[sb/2,sb/2,0],[r,0,0]], #9
+ [[3.*sb/4,3*sb/4,sb/4],[sb/2,sb,sb/2],[r,0,0]], #10
+ [[sb/4,3.*sb/4,sb/4],[sb/2,sb,sb/2],[r,0,0]], #14
+ [[sb/4,3.*sb/4,sb/4],[0,sb/2,sb/2],[r,0,0]], #15
+ [[sb/4,sb/4,3.*sb/4],[0,0,sb],[r,0,0]], #16
+ [[sb/4,sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #17
+ [[sb/4,sb/4,3.*sb/4],[sb/2,0,sb/2],[r,0,0]], #18
+ [[3*sb/4,sb/4,3.*sb/4],[sb/2,0,sb/2],[r,0,0]], #22
+ [[3*sb/4,sb/4,3.*sb/4],[sb,sb/2,sb/2],[r,0,0]], #23
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb/2,sb,sb/2],[r,0,0]], #26
+ [[3*sb/4,3.*sb/4,3.*sb/4],[sb,sb/2,sb/2],[r,0,0]], #27
+ [[sb/4,3.*sb/4,3*sb/4],[0,sb,sb],[r,0,0]], #28
+ [[sb/4,3.*sb/4,3.*sb/4],[sb/2,sb/2,sb],[r,0,0]], #29
+ [[sb/4,3.*sb/4,3.*sb/4],[sb/2,sb,sb/2],[r,0,0]]] #30
+
+
+
+elif name == "5":
+
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+ #angles: 56º and 45º
+ #inner star points:
+
+ s56_f=(sb/2)+((m.sin(m.radians(34))/m.sin(m.radians(56)))*sb/4)
+ s56_c=(sb/2)-((m.sin(m.radians(34))/m.sin(m.radians(56)))*sb/4)
+
+ A=[[[sb,3*sb/4,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #2
+ [[sb,3*sb/4,sb],[3*sb/4,sb/2,sb/2],[r,0,0]], #3
+ [[sb,sb/4,sb],[s56_f,sb/4,sb/2],[r,0,0]], #6
+ [[sb,sb/4,sb],[3*sb/4,sb/2,sb/2],[r,0,0]], #7
+ [[0,3*sb/4,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #8
+ [[0,3*sb/4,0],[sb/4,sb/2,sb/2],[r,0,0]], #9
+ [[0,sb/4,0],[s56_c,sb/4,sb/2],[r,0,0]], #12
+ [[0,sb/4,0],[sb/4,sb/2,sb/2],[r,0,0]], #13
+ [[sb/2,0,0],[s56_c,sb/4,sb/2],[r,0,0]], #16
+ [[sb/2,0,0],[s56_f,sb/4,sb/2],[r,0,0]], #17
+ [[sb/2,sb,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #22
+ [[sb/2,sb,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #23
+ [[sb/2,sb/2,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #24
+ [[sb/2,sb/2,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #26
+ [[sb/2,sb/2,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #28
+ [[sb/2,sb/2,sb],[s56_f,sb/4,sb/2],[r,0,0]], #31
+ [[sb/2,sb/2,sb],[s56_c,sb/4,sb/2],[r,0,0]], #33
+ [[sb/2,sb/2,sb],[sb/4,sb/2,sb/2],[r,0,0]], #35
+
+ [[sb,3*sb/4,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #0
+ [[sb,3*sb/4,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #1
+ [[sb,sb/4,0],[s56_f,sb/4,sb/2],[r,0,0]], #4
+ [[sb,sb/4,0],[3*sb/4,sb/2,sb/2],[r,0,0]], #5
+ [[0,3*sb/4,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #10
+ [[0,3*sb/4,sb],[sb/4,sb/2,sb/2],[r,0,0]], #11
+ [[0,sb/4,sb],[s56_c,sb/4,sb/2],[r,0,0]], #14
+ [[0,sb/4,sb],[sb/4,sb/2,sb/2],[r,0,0]], #15
+ [[sb/2,0,sb],[s56_c,sb/4,sb/2],[r,0,0]], #18
+ [[sb/2,0,sb],[s56_f,sb/4,sb/2],[r,0,0]], #19
+ [[sb/2,sb,0],[s56_c,3*sb/4,sb/2],[r,0,0]], #20
+ [[sb/2,sb,0],[s56_f,3*sb/4,sb/2],[r,0,0]], #21
+ [[sb/2,sb/2,0],[s56_f,sb/4,sb/2],[r,0,0]], #25
+ [[sb/2,sb/2,0],[s56_c,sb/4,sb/2],[r,0,0]], #27
+ [[sb/2,sb/2,0],[sb/4,sb/2,sb/2],[r,0,0]], #29
+ [[sb/2,sb/2,sb],[s56_f,3*sb/4,sb/2],[r,0,0]], #30
+ [[sb/2,sb/2,sb],[s56_c,3*sb/4,sb/2],[r,0,0]], #32
+ [[sb/2,sb/2,sb],[3*sb/4,sb/2,sb/2],[r,0,0]]] #34
+
+
+
+elif name == "6":
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+ A=[[[0,0,0],[sb,0,sb],[r,0,0]],
+ [[0,0,sb],[sb,0,0],[r,0,0]],
+ [[sb,0,sb],[sb,sb,0],[r,0,0]],
+ [[sb,0,0],[sb,sb,sb],[r,0,0]],
+ [[sb,sb,sb],[0,sb,0],[r,0,0]],
+ [[0,sb,sb],[sb,sb,0],[r,0,0]],
+ [[0,sb,sb],[0,0,0],[r,0,0]],
+ [[0,0,sb],[0,sb,0],[r,0,0]],
+ [[0,0,0],[sb,sb,0],[r,0,0]],
+ [[sb,0,0],[0,sb,0],[r,0,0]],
+ [[0,0,sb],[sb,sb,sb],[r,0,0]],
+ [[sb,0,sb],[0,sb,sb],[r,0,0]],
+ [[sb/2,0,sb/2],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]],
+ [[sb/2,0,sb/2],[sb/2,sb/2,0],[r,0,0]],
+ [[sb/2,0,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[sb/2,sb/2,0],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb/2,0],[0,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb/2,0],[sb/2,sb,sb/2],[r,0,0]],
+ [[0,sb/2,sb/2],[sb/2,sb,sb/2],[r,0,0]],
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]],
+ [[sb/2,sb,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[0,sb/2,sb/2],[sb/2,sb/2,sb],[r,0,0]],
+ [[sb,sb/2,sb/2],[sb/2,sb/2,sb],[r,0,0]]]
+
+
+elif name== "7":
+
+ print("5. Not cube, set size of unit cell of y and z direction")
+ sb=float(input())
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ #hexagon side lenght: l
+ #hexagon side to side lenght: sb
+ #lattice x lenght: l2
+
+ l=(sb/2)/(m.sin(m.radians(60)))
+ x_cell=3*l
+
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+
+ A=[[[l,sb,0],[2*l,sb,0],[r,0,0]], #0
+ [[l,0,0],[2*l,0,0],[r,0,0]], #3
+ [[2.5*l,sb/2,0],[3*l,sb/2,0],[r,0,0]], #7
+ [[2.5*l,sb/2,sb],[3*l,sb/2,sb],[r,0,0]], #23
+ [[l,sb,0],[l/2,sb/2,sb/2],[r,0,0]], #24
+ [[2.5*l,sb/2,0],[2*l,sb,sb/2],[r,0,0]], #29
+ [[l,0,sb],[l/2,sb/2,sb/2],[r,0,0]], #32
+ [[2.5*l,sb/2,sb],[2*l,0,sb/2],[r,0,0]], #34
+ [[l/2,sb/2,sb/2],[-l/2,sb/2,0],[r,0,0]], #38
+ [[l/2,sb/2,sb/2],[-l/2,sb/2,sb],[r,0,0]], #39
+
+ [[l,sb,0],[l/2,sb/2,0],[r,0,0]], #1
+ [[l/2,sb/2,0],[l,0,0],[r,0,0]], #2
+ [[2*l,0,0],[2.5*l,sb/2,0],[r,0,0]], #4
+ [[2.5*l,sb/2,0],[2*l,sb,0],[r,0,0]], #5
+
+ [[l,sb,sb/2],[l/2,sb/2,sb/2],[r,0,0]], #9
+ [[l/2,sb/2,sb/2],[l,0,sb/2],[r,0,0]], #10
+ [[2*l,0,sb/2],[2.5*l,sb/2,sb/2],[r,0,0]], #12
+ [[2.5*l,sb/2,sb/2],[2*l,sb,sb/2],[r,0,0]], #13
+ [[l,sb,sb],[l/2,sb/2,sb],[r,0,0]], #17
+ [[l/2,sb/2,sb],[l,0,sb],[r,0,0]], #18
+ [[2*l,0,sb],[2.5*l,sb/2,sb],[r,0,0]], #20
+ [[2.5*l,sb/2,sb],[2*l,sb,sb],[r,0,0]], #21
+
+
+
+ [[l/2,sb/2,0],[0,sb/2,0],[r,0,0]], #6
+ [[l,sb,sb/2],[2*l,sb,sb/2],[r,0,0]], #8
+ [[l,0,sb/2],[2*l,0,sb/2],[r,0,0]], #11
+ [[l/2,sb/2,sb/2],[0,sb/2,sb/2],[r,0,0]], #14
+ [[3*l,sb/2,sb/2],[2.5*l,sb/2,sb/2],[r,0,0]], #15 origin & end coordinates changed: meshing problems
+ [[l,sb,sb],[2*l,sb,sb],[r,0,0]], #16
+ [[l,0,sb],[2*l,0,sb],[r,0,0]], #19
+ [[l/2,sb/2,sb],[0,sb/2,sb],[r,0,0]], #22
+
+ [[l,sb,0],[2*l,sb,sb/2],[r,0,0]], #25
+ [[l,0,0],[l/2,sb/2,sb/2],[r,0,0]], #26
+ [[l,0,0],[2*l,0,sb/2],[r,0,0]], #27
+ [[2.5*l,sb/2,0],[2*l,0,sb/2],[r,0,0]], #28
+ [[l,sb,sb],[l/2,sb/2,sb/2],[r,0,0]], #30
+ [[l,sb,sb],[2*l,sb,sb/2],[r,0,0]], #31
+ [[l,0,sb],[2*l,0,sb/2],[r,0,0]], #33
+ [[2.5*l,sb/2,sb],[2*l,sb,sb/2],[r,0,0]], #35
+ [[2.5*l,sb/2,0],[3.5*l,sb/2,sb/2],[r,0,0]], #36
+ [[2.5*l,sb/2,sb],[3.5*l,sb/2,sb/2],[r,0,0]]] #37
+
+
+
+
+
+elif name== "8":
+
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to"+ str(sb))
+ r=float(input())
+
+ #bow height: bow_h
+ #bow points: bow_p
+ #middle point:(connection) tie_l
+ #tie points:(connection) tie_p
+ #angle of bow corner: bow_angle
+ #f and c correspond to close and far from axis x,y,z
+
+ #inputs:
+ #bow_h=float(input())
+ bow_h=0.75
+ print ("bow heinght between 0 and size of unit cell, chosen:", bow_h)
+ #bow_angle=float(input())
+ bow_angle=75
+ print ("bow angle between 60º and 90º, chosen: "+ str(bow_angle)+"º")
+
+ #-------calculations
+
+ tie_l=((m.sin(m.radians(90-bow_angle))/m.sin(m.radians(bow_angle)))*sb/2)
+ bow_pc=(sb-bow_h)/2
+ bow_pf=bow_h+bow_pc
+ tie_pc=bow_pc+tie_l
+ tie_pf=bow_pf-tie_l
+
+ # only for B
+ m=(bow_pc-tie_pc)/((sb)-(sb/2))
+ pcortec=(sb/4)*m+bow_pc
+ pcortef=(-sb/4)*m+bow_pf
+ centrec=(-sb/10)*m+tie_pc
+ centref=(sb/10)*m+tie_pf
+ #-----------
+
+ # A for vox and phase
+ # B for mesh
+
+ A=[[[tie_pc,sb/2,0],[0,sb/2,0],[r,0,0]],
+ [[tie_pc,sb/2,0],[bow_pc,0,0],[r,0,0]],
+ [[tie_pc,sb/2,0],[bow_pc,sb,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[sb,sb/2,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[bow_pf,0,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[bow_pf,sb,0],[r,0,0]],
+ [[bow_pc,0,0],[bow_pf,0,0],[r,0,0]],
+ [[bow_pc,sb,0],[bow_pf,sb,0],[r,0,0]],
+ [[tie_pc,sb/2,sb],[0,sb/2,sb],[r,0,0]],
+ [[tie_pc,sb/2,sb],[bow_pc,0,sb],[r,0,0]],
+ [[tie_pc,sb/2,sb],[bow_pc,sb,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[sb,sb/2,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[bow_pf,0,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[bow_pc,0,sb],[bow_pf,0,sb],[r,0,0]],
+ [[bow_pc,sb,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb/2,0,0],[r,0,0]],
+ [[sb/2,0,tie_pc],[0,0,bow_pc],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb,0,bow_pc],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb/2,0,sb],[r,0,0]],
+ [[sb/2,0,tie_pf],[0,0,bow_pf],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb,0,bow_pf],[r,0,0]],
+ [[0,0,bow_pc],[0,0,bow_pf],[r,0,0]],
+ [[sb,0,bow_pc],[sb,0,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb/2,sb,0],[r,0,0]],
+ [[sb/2,sb,tie_pc],[0,sb,bow_pc],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb,sb,bow_pc],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb/2,sb,sb],[r,0,0]],
+ [[sb/2,sb,tie_pf],[0,sb,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb,sb,bow_pf],[r,0,0]],
+ [[0,sb,bow_pc],[0,sb,bow_pf],[r,0,0]],
+ [[sb,sb,bow_pc],[sb,sb,bow_pf],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,0,sb/2],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,bow_pc,0],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,bow_pc,sb],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,sb,sb/2],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,bow_pf,0],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,bow_pf,sb],[r,0,0]],
+ [[0,bow_pc,0],[0,bow_pf,0],[r,0,0]],
+ [[0,bow_pc,sb],[0,bow_pf,sb],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,0,sb/2],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,bow_pc,0],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,bow_pc,sb],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,sb,sb/2],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,bow_pf,0],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,bow_pf,sb],[r,0,0]],
+ [[sb,bow_pc,0],[sb,bow_pf,0],[r,0,0]],
+ [[sb,bow_pc,sb],[sb,bow_pf,sb],[r,0,0]]]
+
+
+#-----------
+ B=[[[tie_pc,sb/2,0],[0,sb/2,0],[r,0,0]],
+ [[tie_pf,sb/2,0],[sb,sb/2,0],[r,0,0]],
+ [[bow_pc,0,0],[bow_pf,0,0],[r,0,0]],
+ [[bow_pc,sb,0],[bow_pf,sb,0],[r,0,0]],
+ [[tie_pc,sb/2,sb],[0,sb/2,sb],[r,0,0]],
+ [[tie_pf,sb/2,sb],[sb,sb/2,sb],[r,0,0]],
+ [[bow_pc,0,sb],[bow_pf,0,sb],[r,0,0]],
+ [[bow_pc,sb,sb],[bow_pf,sb,sb],[r,0,0]],
+ [[sb/2,0,tie_pc],[sb/2,0,0],[r,0,0]],
+ [[sb/2,0,tie_pf],[sb/2,0,sb],[r,0,0]],
+ [[0,0,bow_pc],[0,0,bow_pf],[r,0,0]],
+ [[sb,0,bow_pc],[sb,0,bow_pf],[r,0,0]],
+ [[sb/2,sb,tie_pc],[sb/2,sb,0],[r,0,0]],
+ [[sb/2,sb,tie_pf],[sb/2,sb,sb],[r,0,0]],
+ [[0,sb,bow_pc],[0,sb,bow_pf],[r,0,0]],
+ [[sb,sb,bow_pc],[sb,sb,bow_pf],[r,0,0]],
+ [[0,tie_pc,sb/2],[0,0,sb/2],[r,0,0]],
+ [[0,tie_pf,sb/2],[0,sb,sb/2],[r,0,0]],
+ [[0,bow_pc,0],[0,bow_pf,0],[r,0,0]],
+ [[0,bow_pc,sb],[0,bow_pf,sb],[r,0,0]],
+ [[sb,tie_pc,sb/2],[sb,0,sb/2],[r,0,0]],
+ [[sb,tie_pf,sb/2],[sb,sb,sb/2],[r,0,0]],
+ [[sb,bow_pc,0],[sb,bow_pf,0],[r,0,0]],
+ [[sb,bow_pc,sb],[sb,bow_pf,sb],[r,0,0]],
+
+ [[centrec,sb/2+(sb/10),0],[pcortec,-sb/4,0],[r,0,0]],
+ [[centrec,sb/2-(sb/10),0],[pcortec,(sb+sb/4),0],[r,0,0]],
+ [[centref,sb/2+(sb/10),0],[pcortef,-sb/4,0],[r,0,0]],
+ [[centref,sb/2-(sb/10),0],[pcortef,(sb+sb/4),0],[r,0,0]],
+
+ [[centrec,sb/2+(sb/10),sb],[pcortec,-sb/4,sb],[r,0,0]],
+ [[centrec,sb/2-(sb/10),sb],[pcortec,(sb+sb/4),sb],[r,0,0]],
+ [[centref,sb/2+(sb/10),sb],[pcortef,-sb/4,sb],[r,0,0]],
+ [[centref,sb/2-(sb/10),sb],[pcortef,(sb+sb/4),sb],[r,0,0]],
+
+ [[sb/2+(sb/10),0,centrec],[-sb/4,0,pcortec],[r,0,0]],
+ [[sb/2-(sb/10),0,centrec],[(sb+sb/4),0,pcortec],[r,0,0]],
+ [[sb/2+(sb/10),0,centref],[-sb/4,0,pcortef],[r,0,0]],
+ [[sb/2-(sb/10),0,centref],[(sb+sb/4),0,pcortef],[r,0,0]],
+
+ [[sb/2+(sb/10),sb,centrec],[-sb/4,sb,pcortec],[r,0,0]],
+ [[sb/2-(sb/10),sb,centrec],[(sb+sb/4),sb,pcortec],[r,0,0]],
+ [[sb/2+(sb/10),sb,centref],[-sb/4,sb,pcortef],[r,0,0]],
+ [[sb/2-(sb/10),sb,centref],[(sb+sb/4),sb,pcortef],[r,0,0]],
+
+ [[0,centrec,sb/2+(sb/10)],[0,pcortec,-sb/4],[r,0,0]],
+ [[0,centrec,sb/2-(sb/10)],[0,pcortec,(sb+sb/4)],[r,0,0]],
+ [[0,centref,sb/2+(sb/10)],[0,pcortef,-sb/4],[r,0,0]],
+ [[0,centref,sb/2-(sb/10)],[0,pcortef,(sb+sb/4)],[r,0,0]],
+
+ [[sb,centrec,sb/2+(sb/10)],[sb,pcortec,-sb/4],[r,0,0]],
+ [[sb,centrec,sb/2-(sb/10)],[sb,pcortec,(sb+sb/4)],[r,0,0]],
+ [[sb,centref,sb/2+(sb/10)],[sb,pcortef,-sb/4],[r,0,0]],
+ [[sb,centref,sb/2-(sb/10)],[sb,pcortef,(sb+sb/4)],[r,0,0]]]
+
+
+elif name== "9":
+
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0]], #13
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0]]] #15
+
+
+
+elif name== "91":
+
+ print("5. Set size of unit cell")
+ sb=float(input())
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ if vfr=="b":
+ print ("6.b. Specify radius of struts, from 0 to "+ str(sb))
+ r=float(input())
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0]], #13
+ [[0,0,sb/2],[0,sb,sb/2],[r,0,0]], #16
+ [[sb,0,sb/2],[sb,sb,sb/2],[r,0,0]], #17
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0]]] #15
+
+
+
+
+elif name=="92":
+
+ print("5. Set size of unit cell")
+ sb=float(input())
+ if vfr=="b":
+ print("HOLLOW geometry, inner radius dimension: 0.8 x outer radius")
+ print ("6.b. Specify outer radius of hollow struts, from 0 to "+ str(sb))
+ r=float(input())
+ ri=r/1.25
+ x_cell=x_cell+sb
+ y_cell=y_cell+sb
+ z_cell=z_cell+sb
+
+ A=[[[sb/2,0,sb/2],[0,sb/2,sb/2],[r,0,0],[ri,0,0]], #1
+ [[sb/2,sb,sb/2],[sb,sb/2,sb/2],[r,0,0],[ri,0,0]], #2
+ [[sb/4,sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #4
+ [[sb/4,sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #5
+ [[sb/2,sb/2,0],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #10
+ [[sb/2,sb/2,0],[3*sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #11
+ [[sb/2,sb/2,sb],[sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #12
+ [[sb/2,sb/2,sb],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #13
+
+ [[sb/2,0,sb/2],[sb, sb/2,sb/2],[r,0,0],[ri,0,0]], #0
+ [[sb/2,sb,sb/2],[0,sb/2,sb/2],[r,0,0],[ri,0,0]], #3
+ [[3*sb/4,3*sb/4,sb/2],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #6
+ [[3*sb/4,3*sb/4,sb/2],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #7
+ [[sb/2,sb/2,0],[sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #8
+ [[sb/2,sb/2,0],[sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]], #9
+ [[sb/2,sb/2,sb],[3*sb/4,sb/4,sb/2],[r,0,0],[ri,0,0]], #14
+ [[sb/2,sb/2,sb],[3*sb/4,3*sb/4,sb/2],[r,0,0],[ri,0,0]]] #15
+
+
+r0=r
+
+
+
+
+if vfr=="a":
+ if name!="8":
+ print( "6.a. Set volume fraction from 0 to 1")
+ vfraction=float(input())
+ if name=="8":
+ print( "6.a. Set volume fraction from 0 to 0.15")
+ vfraction=float(input())
+ voxelized(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B)
+ print ("7.a. Please introduce value of RADIUSf")
+ #r0=input()
+ with open('Radius.dat','r') as data:
+ B_string = data.read()
+ r0=float(B_string)
+ mesh(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B, r0)
+
+if vfr=="b":
+ vfraction=0.
+ mesh(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B, r0)
+ voxelized(name, vfr, vfraction, rho, sb, x_cell,y_cell,z_cell,length, width,height, A, B)
+
+print(name, rho)'''
+ #plt.show()
+
diff --git a/Lattice_Cell_script/utils_geo.py b/Lattice_Cell_script/utils_geo.py
new file mode 100644
index 0000000000000000000000000000000000000000..dadb16397315cb9c1296c8c5b9031f6d456273eb
--- /dev/null
+++ b/Lattice_Cell_script/utils_geo.py
@@ -0,0 +1,342 @@
+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()
+
diff --git a/Lattice_Cell_script/vmf_bisection.f95 b/Lattice_Cell_script/vmf_bisection.f95
new file mode 100644
index 0000000000000000000000000000000000000000..eb53463f3d3720a1c1fd647e0ebbbaf58fdb0dde
--- /dev/null
+++ b/Lattice_Cell_script/vmf_bisection.f95
@@ -0,0 +1,116 @@
+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
diff --git a/Lattice_Cell_script/vmf_hollow_ID.f95 b/Lattice_Cell_script/vmf_hollow_ID.f95
new file mode 100644
index 0000000000000000000000000000000000000000..665768d95bfe69717be413c66618cb0e57349622
--- /dev/null
+++ b/Lattice_Cell_script/vmf_hollow_ID.f95
@@ -0,0 +1,84 @@
+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
diff --git a/Lattice_Cell_script/vmf_hollow_OD.f95 b/Lattice_Cell_script/vmf_hollow_OD.f95
new file mode 100644
index 0000000000000000000000000000000000000000..f8400a210ab3765f527abaa19d4355a1a97bf952
--- /dev/null
+++ b/Lattice_Cell_script/vmf_hollow_OD.f95
@@ -0,0 +1,83 @@
+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
diff --git a/Lattice_Cell_script/vmf_test.f95 b/Lattice_Cell_script/vmf_test.f95
new file mode 100644
index 0000000000000000000000000000000000000000..49d7323ec2f402434929f66d60beca007b4b6d20
--- /dev/null
+++ b/Lattice_Cell_script/vmf_test.f95
@@ -0,0 +1,96 @@
+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