#!/usr/bin/env python
# coding: utf-8

import gmsh
import sys
import math
import json

ExportPATH="/home/banerjee/Salome/forma11/"

gmsh.initialize()


# Choose OpenCascade kernel
model = gmsh.model
occ = model.occ
mesh = model.mesh
model.add("forma11c_gmsh")

# Default size
el_size = 0.7

# Read the profile coordinates
file_id = open( r''+ExportPATH+'forma11c_profile.json'+'', 'r')
coords = json.load(file_id)
file_id.close()

# Create profile points
v_profile = []
for coord in coords:
    v = occ.addPoint(coord[0], coord[1], coord[2], el_size)
    v_profile.append(v)

# Create spline going through profile points
l1 = occ.addBSpline(v_profile)

# Create copies and rotate
l2 = occ.copy([(1,l1)])
l3 = occ.copy([(1,l1)])
l4 = occ.copy([(1,l1)])

# Rotate the copies
occ.rotate(l2, 0, 0, 0, 0, 0, 1, math.pi/2)
occ.rotate(l3, 0, 0, 0, 0, 0, 1, math.pi)
occ.rotate(l4, 0, 0, 0, 0, 0, 1, 3*math.pi/2)

# Sweep the lines
surf1 = occ.revolve([(1, l1)], 0, 0, 0, 0, 0, 1, math.pi/2)
surf2 = occ.revolve(l2, 0, 0, 0, 0, 0, 1, math.pi/2)
surf3 = occ.revolve(l3, 0, 0, 0, 0, 0, 1, math.pi/2)
surf4 = occ.revolve(l4, 0, 0, 0, 0, 0, 1, math.pi/2)

# Join the surfaces
surf5 = occ.fragment(surf1, surf2)
surf6 = occ.fragment(surf3, surf4)
surf7 = occ.fragment(surf5[0], surf6[0])

# Update
occ.synchronize()


# Get all geometric entities
faces = [face[1] for face in occ.getEntities(2)]
edges = [curve[1] for curve in occ.getEntities(1)]         
vertices = [point[1] for point in occ.getEntities(0)]

# Add physical groups for the mesh
chimney_faces = model.addPhysicalGroup(2, faces)
model.setPhysicalName(2, chimney_faces, "chimney_faces")
chimney_edges = model.addPhysicalGroup(1, edges)
model.setPhysicalName(1, chimney_edges, "chimney_edges")
chimney_verts = model.addPhysicalGroup(0, vertices)
model.setPhysicalName(0, chimney_verts, "chimney_verts")

# Add physical group for the base edges
bb1 = model.getBoundingBox(2, 1)
bb2 = model.getBoundingBox(2, 2)
bb3 = model.getBoundingBox(2, 3)
bb4 = model.getBoundingBox(2, 4)
xmax = bb1[3]
xmin = bb2[0]
ymax = xmax
ymin = xmin
zmin = bb1[2]
zmax = zmin

eps = 1.0e-3
base_edges = model.getEntitiesInBoundingBox(xmin-eps, ymin-eps, zmin-eps, xmax+eps, ymax+eps, zmax+eps, dim=1)
edges = [curve[1] for curve in base_edges]
chimney_base = model.addPhysicalGroup(1, edges)
model.setPhysicalName(1, chimney_base, "chimney_base")

# Assign transfinite interpolation to enetities
num_nodes_circ = 15
for curve in occ.getEntities(1):
    mesh.setTransfiniteCurve(curve[1], num_nodes_circ)

num_nodes_vert = 32
vertical_curves = [7, 10, 13, 17]
for curve in vertical_curves:
    mesh.setTransfiniteCurve(curve, num_nodes_vert)

for surf in occ.getEntities(2):
    mesh.setTransfiniteSurface(surf[1])

# Set up meshing options
gmsh.option.setNumber('Mesh.RecombineAll', 1)
gmsh.option.setNumber('Mesh.RecombinationAlgorithm', 1)
gmsh.option.setNumber('Mesh.Recombine3DLevel', 2)
gmsh.option.setNumber('Mesh.ElementOrder', 1)
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)

# Generate mesh
mesh.generate(2)

# Save mesh
gmsh.write(r''+ExportPATH+'forma11c_gmsh.msh'+'')

# Check in GUI
gmsh.fltk.run()

# Clean up
gmsh.finalize()