import salome
salome.salome_init()

theStudy = salome.myStudy

import time
import sys
sys.path.insert( 0, r'/home/mattkoch/Desktop/CodeSaturne-Turbomachinery')


from math import pi
from math import sqrt
from math import cos, acos
from math import sin, asin

# === Science & Technology Consultants (SciTeX) ===== Tue-07-17-2018 ===
# Geometry Functions
# ======================================================================

# --- Functions --------------------------------------------------------

def RotateCoordinates(PosX,PosY,Ang):
  PosXAng =  PosX*cos(Ang) + PosY*sin(Ang)
  PosYAng = -PosX*sin(Ang) + PosY*cos(Ang)

  return PosXAng, PosYAng


def MakeHelix(Hlx):
  Spl = []

  for Idx in range(Hlx.Cnt):
    PosX = Hlx.Rds
    PosY = 0.0

    Pos0X,Pos0Y = RotateCoordinates(PosX,PosY,Hlx.Rot + Hlx.Ang * Idx/(Hlx.Cnt-1))
    PosZ =                                    Hlx.Off + Hlx.Len * Idx/(Hlx.Cnt-1)

    Pnt = geompy.MakeVertex(Pos0X, Pos0Y, PosZ)
    Spl.append(Pnt)

  Hlx.Obj = geompy.MakeInterpol(Spl)

  return Hlx


def MakeShaft(Shf,Bld):
  Vtr = geompy.MakeVectorDXDYDZ(0.0, 0.0, 1.0)

  Pnt = geompy.MakeVertex(0.0, 0.0, 0.0)
  Shf.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia, Bld.Off - Bld.Thk))

  Pnt = geompy.MakeVertex(0.0, 0.0, Bld.Off - Bld.Thk)
  Shf.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia, Bld.Len + 2.0*Bld.Thk))

  Pnt = geompy.MakeVertex(0.0, 0.0, Bld.Off + Bld.Len + Bld.Thk)
  Shf.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia, Shf.Len - Bld.Off - Bld.Len - Bld.Thk))

  return Shf


def MakeBlade(Bld,Hlx):
  Ang = asin(Bld.Thk/Shf.Dia)

#  PosXMin =  0.5*Shf.Dia*cos(Ang)
  PosXMin =  0.5*Shf.Dia - Bld.Thk
  PosXMax =  0.5*Shf.Dia + Bld.Hgt
  PosYMin = -0.5*Bld.Thk
  PosYMax =  0.5*Bld.Thk
  PosZ = Bld.Off

  Pos0X,Pos0Y = RotateCoordinates(PosXMin,PosYMin,Hlx.Rot)
  Pos1X,Pos1Y = RotateCoordinates(PosXMax,PosYMin,Hlx.Rot)
  Pos2X,Pos2Y = RotateCoordinates(PosXMax,PosYMax,Hlx.Rot)
  Pos3X,Pos3Y = RotateCoordinates(PosXMin,PosYMax,Hlx.Rot)

  Pnt0 = geompy.MakeVertex(Pos0X, Pos0Y, PosZ)
  Pnt1 = geompy.MakeVertex(Pos1X, Pos1Y, PosZ)
  Pnt2 = geompy.MakeVertex(Pos2X, Pos2Y, PosZ)
  Pnt3 = geompy.MakeVertex(Pos3X, Pos3Y, PosZ)

  Edg0 = geompy.MakeEdge(Pnt0, Pnt1)
  Edg1 = geompy.MakeEdge(Pnt1, Pnt2)
  Edg2 = geompy.MakeEdge(Pnt2, Pnt3)
  Edg3 = geompy.MakeEdge(Pnt3, Pnt0)

  Wre = geompy.MakeWire([Edg0, Edg1, Edg2, Edg3])
  Fce = geompy.MakeFace(Wre, True)
  Vtr = geompy.MakeVectorDXDYDZ(0.0, 0.0, 1.0)
  Bld.Obj.append(geompy.MakePipeBiNormalAlongVector(Fce, Hlx.Obj, Vtr))

  return Bld


def MakeFluid(Shf,Bld):
  Vtr = geompy.MakeVectorDXDYDZ(0.0, 0.0, 1.0)

  Pnt = geompy.MakeVertex(0.0, 0.0, 0.0)
  Fld.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia + Bld.Hgt + Bld.Thk, Bld.Off - Bld.Thk))

  Pnt = geompy.MakeVertex(0.0, 0.0, Bld.Off - Bld.Thk)
  Fld.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia + Bld.Hgt + Bld.Thk, Bld.Len + 2.0*Bld.Thk))

  Pnt = geompy.MakeVertex(0.0, 0.0, Bld.Off + Bld.Len + Bld.Thk)
  Fld.Obj.append(geompy.MakeCylinder(Pnt, Vtr, 0.5*Shf.Dia + Bld.Hgt + Bld.Thk, Shf.Len - Bld.Off - Bld.Len - Bld.Thk))

  return Fld

# =============================================== mattkoch@scitex.us ===


# === Science & Technology Consultants (SciTeX) ===== Tue-07-17-2018 ===
# Geometry Generation
# ======================================================================

import GEOM
from salome.geom import geomBuilder
geompy = geomBuilder.New(salome.myStudy)

# --- Classes ----------------------------------------------------------

class Conversions:
  in_m =  25.4e-03; # [m/in]
  ft_m = 304.8e-03; # [m/ft]

class Shaft:
  def __init__(self):
    self.Obj = []

class Helix:
  def __init__(self):
    self.Obj = []

class Blade:
  def __init__(self):
    self.Obj = []

class Fluid:
  def __init__(self):
    self.Obj = []

# --- Instances --------------------------------------------------------

Cvn = Conversions()

Shf = Shaft()

Shf.Dia = 1.500*Cvn.in_m # [in]
Shf.Len = 3.000*Cvn.ft_m # [ft]

Bld = Blade()

Bld.Thk = 0.125*Cvn.in_m # [in]
Bld.Hgt = 3.000*Cvn.in_m # [in]
Bld.Len = 1.000*Cvn.ft_m # [ft]
Bld.Trn = 1 # [1]
Bld.Cnt = 3 # [1]
Bld.Off = 0.5*(Shf.Len - Bld.Len)

Hlx = Helix()

Hlx.Cnt = 25
Hlx.Rds = 0.5*Shf.Dia
Hlx.Len =     Bld.Len
Hlx.Ang =     Bld.Trn*2*pi
Hlx.Off =     Bld.Off

Fld = Fluid()

# --- Calculations -----------------------------------------------------

print("Start Geometry = ",time.time())

Shf = MakeShaft(Shf, Bld)

geompy.addToStudy(Shf.Obj[0], "Shaft-0")
geompy.addToStudy(Shf.Obj[1], "Shaft-1")
geompy.addToStudy(Shf.Obj[2], "Shaft-2")

for Idx in range(Bld.Cnt):
  Hlx.Rot = float(Idx)/float(Bld.Cnt)*2*pi

  Hlx = MakeHelix(Hlx)
  Bld = MakeBlade(Bld,Hlx)

  geompy.addToStudy(Hlx.Obj, "Helix-"+str(Idx))
  geompy.addToStudy(Bld.Obj[Idx], "Blade-"+str(Idx))

Fld = MakeFluid(Shf, Bld)

geompy.addToStudy(Fld.Obj[0], "Fluid-0")
geompy.addToStudy(Fld.Obj[1], "Fluid-1")
geompy.addToStudy(Fld.Obj[2], "Fluid-2")

print("Stop Geometry = ",time.time())

# --- Volumes ----------------------------------------------------------

Blade = geompy.MakeFuseList([Shf.Obj[1], Bld.Obj[0], Bld.Obj[1], Bld.Obj[2]], True, True)

Inlet = geompy.MakeCutList(Fld.Obj[0], [Shf.Obj[0]], True)

Turbine = geompy.MakeCutList(Fld.Obj[1], [Blade], True)

Outlet = geompy.MakeCutList(Fld.Obj[2], [Shf.Obj[2]], True)

Domain = geompy.MakePartition([Inlet, Turbine, Outlet], [], [], [], geompy.ShapeType["SOLID"], 0, [], 0)

# --- Faces ------------------------------------------------------------

Inlet_Inlet = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet_Inlet, [11])
Inlet_Shaft = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet_Shaft, [21])
Inlet_Wall = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet_Wall, [4])
Inlet_Turbine = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet_Turbine, [16])

Outlet_Turbine = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet_Turbine, [31])
Outlet_Shaft = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet_Shaft, [163])
Outlet_Wall = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet_Wall, [153])
Outlet_Outlet = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet_Outlet, [158])

Turbine_Shaft = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Turbine_Shaft, [36, 65])
Turbine_Wall = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Turbine_Wall, [26])
Turbine_Blade = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Turbine_Blade, [131, 97, 114, 109, 92, 126, 145, 149, 147, 136, 121, 142, 85, 104, 139])

# --- Solids -----------------------------------------------------------

Domain_Inlet = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(Domain_Inlet, [2])
Domain_Turbine = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(Domain_Turbine, [24])
Domain_Outlet = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(Domain_Outlet, [151])

# --- Add to Study -----------------------------------------------------

geompy.addToStudy(Inlet, "Inlet")
geompy.addToStudy(Outlet, "Outlet")
geompy.addToStudy(Blade, "Blade")
geompy.addToStudy(Turbine, "Turbine")
geompy.addToStudy(Domain, "Domain")

geompy.addToStudyInFather(Domain, Inlet_Inlet, 'Inlet_Inlet')
geompy.addToStudyInFather(Domain, Inlet_Wall, 'Inlet_Wall')
geompy.addToStudyInFather(Domain, Inlet_Shaft, 'Inlet_Shaft')
geompy.addToStudyInFather(Domain, Inlet_Turbine, 'Inlet_Turbine')

geompy.addToStudyInFather(Domain, Outlet_Turbine, 'Outlet_Turbine')
geompy.addToStudyInFather(Domain, Outlet_Wall, 'Outlet_Wall')
geompy.addToStudyInFather(Domain, Outlet_Shaft, 'Outlet_Shaft')
geompy.addToStudyInFather(Domain, Outlet_Outlet, 'Outlet_Outlet')

geompy.addToStudyInFather(Domain, Turbine_Shaft, 'Turbine_Shaft')
geompy.addToStudyInFather(Domain, Turbine_Wall, 'Turbine_Wall')
geompy.addToStudyInFather(Domain, Turbine_Blade, 'Turbine_Blade')

# =============================================== mattkoch@scitex.us ===


# === Science & Technology Consultants (SciTeX) ===== Tue-07-17-2018 ===
# Mesh Generation
# ======================================================================

import  SMESH, SALOMEDS
from salome.smesh import smeshBuilder

smesh = smeshBuilder.New(theStudy)
Domain_Mesh = smesh.Mesh(Domain)

NETGEN_1D_2D_3D = Domain_Mesh.Tetrahedron(algo=smeshBuilder.NETGEN_1D2D3D)
NETGEN_3D_Settings = NETGEN_1D_2D_3D.Parameters()
NETGEN_3D_Settings.SetMaxSize( 0.01 )
NETGEN_3D_Settings.SetSecondOrder( 0 )
NETGEN_3D_Settings.SetOptimize( 1 )
NETGEN_3D_Settings.SetFineness( 3 )
NETGEN_3D_Settings.SetMinSize( 0.005 )
NETGEN_3D_Settings.SetUseSurfaceCurvature( 1 )
NETGEN_3D_Settings.SetFuseEdges( 1 )
NETGEN_3D_Settings.SetQuadAllowed( 0 )

print("Start Mesh = ",time.time())

isDone = Domain_Mesh.Compute()

print("Stop Mesh = ",time.time())

# --- Volumes ----------------------------------------------------------

Domain_Inlet_Volume = Domain_Mesh.GroupOnGeom(Domain_Inlet,'Domain_Inlet',SMESH.VOLUME)
Domain_Turbine_Volume = Domain_Mesh.GroupOnGeom(Domain_Turbine,'Domain_Turbine',SMESH.VOLUME)
Domain_Outlet_Volume = Domain_Mesh.GroupOnGeom(Domain_Outlet,'Domain_Outlet',SMESH.VOLUME)

# --- Faces ------------------------------------------------------------

Inlet_Inlet_Faces = Domain_Mesh.GroupOnGeom(Inlet_Inlet,'Inlet_Inlet',SMESH.FACE)
Inlet_Shaft_Faces = Domain_Mesh.GroupOnGeom(Inlet_Shaft,'Inlet_Shaft',SMESH.FACE)
Inlet_Wall_Faces = Domain_Mesh.GroupOnGeom(Inlet_Wall,'Inlet_Wall',SMESH.FACE)
Inlet_Turbine_Faces = Domain_Mesh.GroupOnGeom(Inlet_Turbine,'Inlet_Turbine',SMESH.FACE)

Outlet_Turbine_Faces = Domain_Mesh.GroupOnGeom(Outlet_Turbine,'Outlet_Turbine',SMESH.FACE)
Outlet_Shaft_Faces = Domain_Mesh.GroupOnGeom(Outlet_Shaft,'Outlet_Shaft',SMESH.FACE)
Outlet_Wall_Faces = Domain_Mesh.GroupOnGeom(Outlet_Wall,'Outlet_Wall',SMESH.FACE)
Outlet_Outlet_Faces = Domain_Mesh.GroupOnGeom(Outlet_Outlet,'Outlet_Outlet',SMESH.FACE)

Turbine_Shaft_Faces = Domain_Mesh.GroupOnGeom(Turbine_Shaft,'Turbine_Shaft',SMESH.FACE)
Turbine_Wall_Faces = Domain_Mesh.GroupOnGeom(Turbine_Wall,'Turbine_Wall',SMESH.FACE)
Turbine_Blade_Faces = Domain_Mesh.GroupOnGeom(Turbine_Blade,'Turbine_Blade',SMESH.FACE)

# --- Nodes ------------------------------------------------------------

Domain_Inlet_Nodes = Domain_Mesh.GroupOnGeom(Domain_Inlet,'Domain_Inlet',SMESH.NODE)
Domain_Turbine_Nodes = Domain_Mesh.GroupOnGeom(Domain_Turbine,'Domain_Turbine',SMESH.NODE)
Domain_Outlet_Nodes = Domain_Mesh.GroupOnGeom(Domain_Outlet,'Domain_Outlet',SMESH.NODE)

Inlet_Inlet_Nodes = Domain_Mesh.GroupOnGeom(Inlet_Inlet,'Inlet_Inlet',SMESH.NODE)
Inlet_Shaft_Nodes = Domain_Mesh.GroupOnGeom(Inlet_Shaft,'Inlet_Shaft',SMESH.NODE)
Inlet_Wall_Nodes = Domain_Mesh.GroupOnGeom(Inlet_Wall,'Inlet_Wall',SMESH.NODE)
Inlet_Turbine_Nodes = Domain_Mesh.GroupOnGeom(Inlet_Turbine,'Inlet_Turbine',SMESH.NODE)

Outlet_Turbine_Nodes = Domain_Mesh.GroupOnGeom(Outlet_Turbine,'Outlet_Turbine',SMESH.NODE)
Outlet_Shaft_Nodes = Domain_Mesh.GroupOnGeom(Outlet_Shaft,'Outlet_Shaft',SMESH.NODE)
Outlet_Wall_Nodes = Domain_Mesh.GroupOnGeom(Outlet_Wall,'Outlet_Wall',SMESH.NODE)
Outlet_Outlet_Nodes = Domain_Mesh.GroupOnGeom(Outlet_Outlet,'Outlet_Outlet',SMESH.NODE)

Turbine_Shaft_Nodes = Domain_Mesh.GroupOnGeom(Turbine_Shaft,'Turbine_Shaft',SMESH.NODE)
Turbine_Wall_Nodes = Domain_Mesh.GroupOnGeom(Turbine_Wall,'Turbine_Wall',SMESH.NODE)
Turbine_Blade_Nodes = Domain_Mesh.GroupOnGeom(Turbine_Blade,'Turbine_Blade',SMESH.NODE)

# =============================================== mattkoch@scitex.us ===