import gmsh
import math
import random
from subprocess import run

gmsh.initialize()

factory = gmsh.model.occ

gmsh.option.setNumber("Mesh.Tetrahedra", 1)
gmsh.option.setNumber("Mesh.Triangles", 1)

gmsh.option.setNumber("General.Terminal", 1)
gmsh.option.setNumber("General.Verbosity", 5)
gmsh.option.setNumber("Mesh.SaveAll", 0)
gmsh.option.setNumber("Mesh.MshFileVersion", 4.1)


fac = 1

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 2 * fac)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2 * fac)

# geometry parameters
voiden = 0.01
Vf = 0.7
#  Coefficent of variation of the fibre diameter
cvf = 0.17
Vlratio = 20
Cube = 75 * fac
frmean = 5.0 * fac
frmin = frmean - (frmean * cvf)
frmax = frmean + (frmean * cvf)
vrmin = 2.5 * fac
vrmax = 5.0 * fac
R = Cube
nSvoid = math.floor(
    (voiden * 2 * R * 2 * R * 2 * R) / ((4 / 3) * math.pi * vrmax * vrmax * vrmax)
)
nFibre = math.floor((Vf * 2 * R * 2 * R * 2 * R) / (math.pi * frmax * frmax * 2 * R))
dims = math.floor(math.sqrt(nFibre))

box_1 = factory.addBox(-R, -R, -R, 2 * R, 2 * R, 2 * R)
spheres = list()
for s in range(2, nSvoid + 1):
    if s < nSvoid:
        r = vrmin + random.uniform(0.0, vrmax - vrmin)
        x = -R + random.uniform(0.0, 2 * R)
        y = -R + random.uniform(0.0, 2 * R)
        z = -R + random.uniform(0.0, 2 * R)
        sphere = factory.addSphere(x, y, z, r)
        spheres.append(sphere)

    elif s == nSvoid:
        r = vrmin + random.uniform(0.0, vrmax - vrmin)
        x = -R + random.uniform(0.0, 2 * R)
        y = -R + random.uniform(0.0, 2 * R)
        z = -R + random.uniform(0.0, 2 * R)
        sphere = factory.addSphere(x, y, z, r)
        factory.dilate([(3, sphere)], x, y, z, Vlratio, 1, 0.5)
        spheres.append(sphere)

phy_gp = gmsh.model.addPhysicalGroup(3, [sphere])

difference = factory.cut(
    [(3, box_1)], [(3, tag) for tag in spheres], removeObject=True, removeTool=True
)
# get the tag :
difference = difference[0][0][1]

dx = (2 * R) / dims
count = 0
# parameters for  cylinder axis
dxF, dyF, dzF = 0, 0, 2 * R

cylinders = list()

for i in range(0, dims):
    for j in range(0, dims):
        # Fibre diameter
        Fr = random.uniform(frmin, frmax)
        # Center of grid box in the RVE
        cgx = (-R + 0.5 * dx) + j * dx
        cgy = (-R + 0.5 * dx) + i * dx
        # Center of cylinder base
        xF = (cgx - 0.5 * dx + Fr) + random.uniform(0.0, dx - 2 * Fr)
        yF = (cgy - 0.5 * dx + Fr) + random.uniform(0.0, dx - 2 * Fr)
        zF = -R
        cylinder = factory.addCylinder(xF, yF, zF, dxF, dyF, dzF, Fr)
        cylinders.append(cylinder)

phy_cylinders = list()
for cyl in cylinders:
    phy_cyl = gmsh.model.addPhysicalGroup(3, [cyl])
    phy_cylinders.append(phy_cyl)

matrix = factory.cut(
    [(3, difference)], [(3, t) for t in cylinders], removeObject=True, removeTool=False
)
matrix = matrix[0][0][1]
phy_matrix = gmsh.model.addPhysicalGroup(3, [matrix])
gmsh.model.setPhysicalName(3, phy_matrix, "Matrix")
phy_fibres = gmsh.model.addPhysicalGroup(3, cylinders)
gmsh.model.setPhysicalName(3, phy_fibres, "Fibres")

all_entities = gmsh.model.getEntities(3)
# remove the matrix and the fibres :
remove_entities = [e for e in all_entities if e[1] not in cylinders + [matrix]]

gmsh.model.removeEntities(remove_entities, recursive=True)
# factory.fragment([(3, matrix)], [(3, t) for t in cylinders])
factory.synchronize()

gmsh.model.mesh.generate(3)
gmsh.model.mesh.removeDuplicateNodes()
gmsh.model.mesh.renumberNodes()
gmsh.model.mesh.renumberElements()

gmsh.write("rve.msh")
try:
    run(f"meshio-convert rve.msh rve.xdmf", shell=True, check=True)
except:
    pass

gmsh.fltk.run()