# ------------------------------------------------------------------------------
#
#  Gmsh Mesh of a cylindrical pipe
#
#  Dimensions:
#
# ------------------------------------------------------------------------------

# The Python API is entirely defined in the `gmsh.py' module (which contains the
# full documentation of all the functions in the API):
import gmsh
import sys
import numpy as np
from math import isclose
# Before using any functions in the Python API, Gmsh must be initialized:
gmsh.initialize()

# Next we add a new model named "pipe1" (if gmsh.model.add() is not called a new
# unnamed model will be created on the fly, if necessary):
gmsh.model.add("pipev1")

# We start by defining some points and some lines. To make the code shorter we
# can redefine a namespace:
fG = gmsh.model.geo

# Pipe Dimensions
Ri=7.0 # in
Ro=8.0 # in 
ft2in=12.0 
Lz=4.0*ft2in
lin=0.5 # Target elements
lout=0.5
nz=50
###############################################################################
# Adapted from t1.py
# The first type of `elementary entity' in Gmsh is a `Point'. To create a point
# with the built-in CAD kernel, the Python API function is
# gmsh.model.geo.addPoint():
# - the first 3 arguments are the point coordinates (x, y, z)
# - the next (optional) argument is the target mesh size close to the point
# - the last (optional) argument is the point tag (a stricly positive integer
#   that uniquely identifies the point)
# 5 points to define the circle origin and circumference
rad=Ri
cin_pts=[]
cin_pts.append(fG.addPoint(0.0,0.0,0.0,lin))
cin_pts.append(fG.addPoint(rad,0.0,0.0,lin))
cin_pts.append(fG.addPoint(-rad,0.0,0.0,lin)) 
cin_pts.append(fG.addPoint(0.0,rad,0.0,lin)) 
cin_pts.append(fG.addPoint(0.0,-rad,0.0,lin)) 

# 5 points to define the outer circle 
rad=Ro
cout_pts=[]
cout_pts.append(fG.addPoint(0.0,0.0,0.0,lout))
cout_pts.append(fG.addPoint(rad,0.0,0.0,lout))
cout_pts.append(fG.addPoint(-rad,0.0,0.0,lout)) 
cout_pts.append(fG.addPoint(0.0,rad,0.0,lout)) 
cout_pts.append(fG.addPoint(0.0,-rad,0.0,lout))

# Curves are Gmsh's second type of elementery entities, and, amongst curves,
# straight lines are the simplest. The API to create straight line segments with
# the built-in kernel follows the same conventions: the first 2 arguments are
# point tags (the start and end points of the line), and the last (optional one)
# is the line tag.
# Note that curve tags are separate from point tags - hence we can reuse tag `1'
# for our first curve. And as a general rule, elementary entity tags in Gmsh
# have to be unique per geometrical dimension.

# Define 4 arcs of the outer circle
cout_arcs=[]
cout_arcs.append(fG.addCircleArc(cout_pts[3],cout_pts[0] , cout_pts[1]))
cout_arcs.append(fG.addCircleArc(cout_pts[1],cout_pts[0] , cout_pts[4]))
cout_arcs.append(fG.addCircleArc(cout_pts[4],cout_pts[0] , cout_pts[2]))
cout_arcs.append(fG.addCircleArc(cout_pts[2],cout_pts[0] , cout_pts[3]))
#print(cout_arcs)
# Define 4 arcs of the inner circle
cin_arcs=[]
cin_arcs.append(fG.addCircleArc(cin_pts[3],cin_pts[0] , cin_pts[1]))
cin_arcs.append(fG.addCircleArc(cin_pts[1],cin_pts[0] , cin_pts[4]))
cin_arcs.append(fG.addCircleArc(cin_pts[4],cin_pts[0] , cin_pts[2]))
cin_arcs.append(fG.addCircleArc(cin_pts[2],cin_pts[0] , cin_pts[3]))
#print(cin_arcs)
# Define 4 lines b/w circles
clines=[fG.addLine(cout_pts[index], cin_pts[index]) for index in range(1,5)]
#print(clines)

# The third elementary entity is the surface. In order to define a simple
# circular surface from the four curves defined above, a curve loop has first
# to be defined. A curve loop is defined by an ordered list of connected curves,
# a sign being associated with each curve (depending on the orientation of the
# curve to form a loop). The API function to create curve loops takes a list
# of integers as first argument, and the curve loop tag (which must be unique
# amongst curve loops) as the second (optional) argument:

quads=[]
quads.append(fG.addCurveLoop([cout_arcs[0],clines[0],-cin_arcs[0],-clines[2] ]))
quads.append(fG.addCurveLoop([cout_arcs[1],clines[3],-cin_arcs[1],-clines[0] ]) )
quads.append(fG.addCurveLoop([cout_arcs[2],clines[1],-cin_arcs[2],-clines[3] ]))
quads.append(fG.addCurveLoop([cout_arcs[3],clines[2],-cin_arcs[3],-clines[1] ]))

surfs=[]
extsurfs=[]
vols=[]
extvec=[0.0,0.0,Lz]
for quad in quads:
    circquad=fG.addPlaneSurface([quad])
    cylquad=fG.extrude([(2,circquad)],extvec[0],extvec[1],extvec[2],[nz],[1],recombine=True)
    surfs.append(circquad)
    extsurfs.append(cylquad[0])
    vols.append(cylquad[1])
fG.synchronize()    


def avgdist(taglist,origin):
    dist=0.0
    for e in taglist:
        nodeCoords = gmsh.model.getValue(e[0], e[1],[])
        dist=dist+np.linalg.norm(nodeCoords-origin)
    return dist/len(taglist)
coutext_arcs=[]
cinext_arcs=[]
for surf in extsurfs:
    topsurfcurves=gmsh.model.getBoundary([surf])
    for c in topsurfcurves:
        bdypts=gmsh.model.getBoundary([c])
        dist=avgdist(bdypts,extvec)
        if isclose(dist,Ro):
            coutext_arcs.append(c[1])
        elif isclose(dist,Ri):
            cinext_arcs.append(c[1])
print(coutext_arcs)
print(cinext_arcs)
#entities = gmsh.model.getEntities()
#for e in entities:
#    print("Entity " + str(e) + " of type " + gmsh.model.getType(e[0], e[1]))


# All the extrusion functions return a list of (dim,tag) tuples of extruded entities: the
# "top" of the extruded surface (in `[0]'), the newly created volume (in
# `[1]') and the lateral surfaces (in `[2]', `[3]', ...).
#print(surfs)
#print(vquads)

# Before they can be meshed (and, more generally, before they can be used by API
# functions outside of the built-in CAD kernel functions), the CAD entities must
# be synchronized with the Gmsh model, which will create the relevant Gmsh data
# structures. This is achieved by the gmsh.model.geo.synchronize() API call for
# the built-in CAD kernel. Synchronizations can be called at any time, but they
# involve a non trivial amount of processing; so while you could synchronize the
# internal CAD data after every CAD command, it is usually better to minimize
# the number of synchronization points.

fG.synchronize()


# At this level, Gmsh knows everything to display the pipe and
# to mesh it. An optional step is needed if we want to group elementary
# geometrical entities into more meaningful groups, e.g. to define some
# mathematical ("domain", "boundary"), functional ("left wing", "fuselage") or
# material ("steel", "carbon") properties.
#
# Such groups are called "Physical Groups" in Gmsh. By default, if physical
# groups are defined, Gmsh will export in output files only mesh elements that
# belong to at least one physical group. (To force Gmsh to save all elements,
# whether they belong to physical groups or not, set the `Mesh.SaveAll' option
# to 1.) Physical groups are also identified by tags, i.e. stricly positive
# integers, that should be unique per dimension (0D, 1D, 2D or 3D). Physical
# groups can also be given names.
#
# Here we define a physical curve that groups the arcs of inner and outer
# circles into individual groups :
loop_cin=fG.addCurveLoop(cin_arcs)    # Inner Circle
loop_cout=fG.addCurveLoop(cout_arcs)  # Outer Circle
# Physical curve groups are started from 1000
# Physical surfaces are started from 2000
# Physical volumes are started from 3000
gmsh.model.addPhysicalGroup(1, cout_arcs,1001)
gmsh.model.addPhysicalGroup(1, cin_arcs,1002)
gmsh.model.addPhysicalGroup(1, coutext_arcs,1003)
gmsh.model.addPhysicalGroup(1, cinext_arcs,1004)

surftags=surfs
extsurftags=[t[1] for t in extsurfs]
voltags=[v[1] for v in vols]

gmsh.model.addPhysicalGroup(2, surftags,2001)
gmsh.model.addPhysicalGroup(2, extsurftags,2002)
gmsh.model.addPhysicalGroup(3, voltags,3001)

#ps = gmsh.model.addPhysicalGroup(2, [1])
#gmsh.model.setPhysicalName(2, ps, "My surface")


#################################################################################
# Adapted from t11.py
# To generate quadrangles instead of triangles, we can simply add
#gmsh.model.mesh.setRecombine(2, pl)
# If we'd had several surfaces, we could have used the global option
gmsh.option.setNumber("Mesh.RecombineAll", 1)

# The default recombination algorithm is called "Blossom": it uses a minimum
# cost perfect matching algorithm to generate fully quadrilateral meshes from
# triangulations. More details about the algorithm can be found in the
# following paper: J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise,
# A. Johnen and C. Geuzaine, "Blossom-Quad: a non-uniform quadrilateral mesh
# generator using a minimum cost perfect matching algorithm", International
# Journal for Numerical Methods in Engineering 89, pp. 1102-1119, 2012.

# For even better 2D (planar) quadrilateral meshes, you can try the
# experimental "Frontal-Delaunay for quads" meshing algorithm, which is a
# triangulation algorithm that enables to create right triangles almost
# everywhere: J.-F. Remacle, F. Henrotte, T. Carrier-Baudouin, E. Bechet,
# E. Marchandise, C. Geuzaine and T. Mouton. A frontal Delaunay quad mesh
# generator using the L^inf norm. International Journal for Numerical Methods
# in Engineering, 94, pp. 494-512, 2013. Uncomment the following line to try
# the Frontal-Delaunay algorithms for quads:
#
gmsh.option.setNumber("Mesh.MshFileVersion", 2)
gmsh.option.setNumber("Mesh.Algorithm", 8)

# You can also set the subdivision step alone, with
# 0=none, 1=all quadrangles, 2=all hexahedra
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2)
gmsh.model.mesh.generate(3)

# Note that you could also apply the recombination algorithm and/or the
# subdivision step explicitly after meshing, as follows:
#
# gmsh.model.mesh.generate(2)
# gmsh.model.mesh.recombine()
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
# gmsh.model.mesh.refine()

# ... and save it to disk
gmsh.write("pipev1.msh")
########################################################################

# To visualize the model we can run the graphical user interface with
# `gmsh.fltk.run()'. Here we run it only if the "-nopopup" is not provided in
# the command line arguments:
if '-nopopup' not in sys.argv:
    gmsh.fltk.run()

gmsh.finalize()

#ToDo
# Check t6.py and use setTransfiniteCurve options to customize mesh 
# Look into Mesh Smoothing option in t6.py