import salome

salome.salome_init()

import salome_notebook
import GEOM
from salome.geom import geomBuilder
from salome.smesh import smeshBuilder
import SALOMEDS
import SMESH

import numpy as np


notebook = salome_notebook.NoteBook()

#=======================================================================================================
#                               Défnition des étapes à effectuer
#=======================================================================================================
generation_parametres_maillage = True
generation_maillage = True

#=======================================================================================================
#                               Génération des aspects géométriques 
#=======================================================================================================
geompy = geomBuilder.New()

# Création des objets de base
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)


# Définittion du domaine
x_min, x_max = -100, 100
y_min, y_max = -100, 100
z_min, z_max = -100, 100

Box = geompy.MakeBoxDXDYDZ(x_max-x_min, y_max-y_min, z_max-z_min)
Box = geompy.MakeTranslation(Box, -(x_max-x_min)/2, -(y_max-y_min)/2, -(z_max-z_min)/2)

Box_compound = geompy.MakeCompound(Box)
geompy.addToStudy(Box_compound, "Box")

maSphere = geompy.MakeSphere(0, 0, 0, 20, "maSphere")

Domaine = geompy.MakeCut(Box_compound, maSphere, True) 

# Create group for finer meshing
Sphere_surface = geompy.CreateGroup(Domaine, geompy.ShapeType["FACE"])           #Création d'un groupe de faces vide pour l'instant appelé Torus_surface 
geompy.UnionIDs(Sphere_surface, [3])                                             #Ajout à ce groupe des faces par leur IDs

Domaine_compound = geompy.MakeCompound([Domaine, Sphere_surface])                 #Création d'un Assemblage comprenant le domaine et les faces de la surface de la géométrie groupées dans un groupe

Box_inside_Domaine = geompy.ExtractShapes(Domaine_compound, geompy.ShapeType["SOLID"])

Sphere_faces_inside_Domaine = geompy.GetSubShape(Domaine_compound, [36], "Sphere_faces_inside_Domaine")

geompy.addToStudy(Domaine, "Domaine")

#=======================================================================================================
#                               Génération des paramètres de maillage 
#=======================================================================================================

if generation_parametres_maillage :

    smesh = smeshBuilder.New()

    print("Creating Parallel Mesh")
    par_mesh = smesh.ParallelMesh(Domaine, name="par_mesh")

    #Création de groupes de maillage : permet de pouvoir identifier les faces après l'export
    mesh_Sphere_surface = par_mesh.GroupOnGeom(Sphere_surface, "Sphere_surface", SMESH.FACE)

    # Assign finer meshing to the surface group
    print("Adding refined hypothesis to Sphere_surface")
    Algo2D_surface = par_mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D, geom=Sphere_faces_inside_Domaine)
    Algo2D_surface_Parameters=Algo2D_surface.Parameters()
    Algo2D_surface_Parameters.SetMaxSize( 2 )
    Algo2D_surface_Parameters.SetMinSize( 0.05 )
    Algo2D_surface_Parameters.SetSecondOrder( 0 )

    print("Adding global 2D hypothesis")
    Algo2D_general = par_mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D)
    Algo2D_general_Parameters=Algo2D_general.Parameters()
    Algo2D_general_Parameters.SetMaxSize( 8 )
    Algo2D_general_Parameters.SetMinSize( 0.05 )
    Algo2D_general_Parameters.SetSecondOrder( 0 )

    print("Adding global hypothesis")
    Algo3D = par_mesh.Tetrahedron(algo=smeshBuilder.NETGEN_3D_Remote) # <------- Change NETGEN_3D_Remote by NETGEN_3D to see viscous layers appear
    Algo3D_Parameters=Algo3D.Parameters()
    Algo3D_Parameters.SetMaxSize( 8 )
    Algo3D_Parameters.SetMinSize( 0.05 )
    Algo3D_Parameters.SetSecondOrder( 0 )


    Viscous_Layers= Algo3D.ViscousLayers(
        1.5,                            #Thickness
        8,                              #Nb of layer
        1.2,                            #stretchFactor
        Sphere_faces_inside_Domaine,    #faces
        False                           #isFaceToIgnore
        #"SURF_OFFSET_SMOOTH"           #extrMethod (SURF_OFFSET_SMOOTH | FACE_SMOOTH | NODE_OFFSET)
    )

    print("Setting parallelism method")
    par_mesh.SetParallelismMethod(smeshBuilder.MULTITHREAD)

    print("Setting parallelism options")
    param = par_mesh.GetParallelismSettings()
    param.SetNbThreads(4)

#=======================================================================================================
#                               Génération du maillage et affichage de ses propriétés
#=======================================================================================================

if generation_maillage:

    print("Starting parallel compute")
    is_done = par_mesh.Compute()
    if not is_done:
        raise Exception("Error when computing Mesh")

    print("  Tetrahedron: ",  par_mesh.NbTetras())
    print("  Triangle: ", par_mesh.NbTriangles())
    print("  Edge: ", par_mesh.NbEdges())

    assert par_mesh.NbTetras() > 0
    assert par_mesh.NbTriangles() > 0
    assert par_mesh.NbEdges() > 0

salome.sg.updateObjBrowser()