# -*- coding: utf-8 -*-

###
### This file is generated automatically by SALOME v8.2.0 with dump python functionality
###

import sys
import salome

salome.salome_init()
theStudy = salome.myStudy

import salome_notebook
notebook = salome_notebook.NoteBook(theStudy)
#sys.path.insert( 0, r'/home/apetras/Documents/radiofrequencyablation-hpc/mesh')

###
### GEOM component
###

import GEOM
from salome.geom import geomBuilder
import math
import os


###
### SMESH component
###

import  SMESH
from salome.smesh import smeshBuilder

# ----------------------- Salome study -----------------------------
geompy = geomBuilder.New(theStudy)

# ---------------------- Global usage ------------------------------
# Cartesian axis
O = geompy.MakeVertex(0, 0, 0)
OX = geompy.MakeVectorDXDYDZ(1, 0, 0)
OY = geompy.MakeVectorDXDYDZ(0, 1, 0)
OZ = geompy.MakeVectorDXDYDZ(0, 0, 1)

# ----------------------- Geometry constants --------------
Lx = 80e-3 # Length of box in x-axis (m)
Ly = 80e-3 # Length of box in y-axis (m)
Lz = 80e-3 # Length of box in z-axis (m)
Lz_tissue = 40e-3 # Length of tissue box in z-axis (m)
Lz_plastic = 20e-3 # Length of tissue box in z-axis (m)
R_catheter = 1.165e-3 # Radius of cylindrical catheter (m)
R_thermistor = 2*R_catheter/3 # Radius of cylindrical thermistor (m)
L_electrode = 3.5e-3 # Length of electrode (m)
R_tube_holes = 0.25e-3 # Radius of electrode holes (m)
R_tube = math.sqrt(6*0.15e-3**2) # Radius of cylindrical tube (m)
L_thermistor = 3e-3 # Length of thermistor (m)
R_sph_cyl = L_electrode/13 # Radius of rounded corners (m)
L_holes = 0.8e-3 # The position of the holes (m)

# Add user defined force and Young's modulus of elasticity
try:
    inp = sys.argv[1]
    inp = map(int, inp.strip('[]').split(','))
    F = int(inp[0])
    E = int(inp[1])
    if F >0 and F <100 :
        F = F*1e-3
        if E >= 50.0:
            E = E*1e3
        else:
            raise """The Young's modulus should be larger or equal to 50kPa"""
    else:
        raise 'The contact force should be between 0 and 100'
except:
    raise 'Provide an input of the form [F in kg, E in kPa]'
#F = 40e-3
#E = 50e3 # Young's modulus of elasticity (Pa)
print('F = '+str(F*1e3)+'g and E = '+str(E*1e-3)+'kPa')

# ----------------------- Calculation of insertion depth (using elasticity) --------------
# This code is a model for a spherical indenter from the paper by Sneddon.
import scipy.optimize as opt 
import scipy.integrate as integr
from math import pi
import numpy as np

def curveShift(x,u,a,u_a,omega_max,f,tol,sl_tol):
    u1 = u[x>a]; x1 = x[x>a];
    sl = (u1-u_a)/(x1-a)
    dist = np.sqrt(np.power(x1-a,2) + np.power(u1 - u_a,2))
    
    for j in range(len(dist)):
        if dist[j] > tol:
            if sl[j]>sl_tol:
                print 'Curve shifting is required. Current slope = '+str(sl[j])
                a = np.min(R_catheter,x1[j])
                u_a = max(-omega_max+f(a,1),u1[j])
                d1 = a - x1[j]
                x1 = x1 + d1
                # Tolerance to avoid nearly flat triangles around tube holes
                if abs(omega_max + u_a - L_holes)<=tol: 
                    u_a = -omega_max+L_holes
                j = np.argmax(u1>u_a)
                x1 = np.concatenate(([a],x1[j:]),axis=0);
                u1 = np.concatenate(([u_a],u1[j:]),axis=0);
                u1 = u1[x1>a]; x1 = x1[x1>a]
                sl = (u1-u_a)/(x1-a)
                dist = np.sqrt(np.power(x1-a,2) + np.power(u1 - u_a,2))
                print 'Curve shifted, new a = '+str(a)+', and new slope = '+str(sl[np.argmax(dist>tol)])
                break
            print 'No need for curve shifting'
            break
    u1 = np.concatenate(([u_a],u1[x1>a]),axis=0)
    x1 = np.concatenate(([a],x1[x1>a]),axis=0)
    return x1,u1,a,u_a

# Parameters
nu = 0.499 # Poisson's ratio
mu = E/(2*(1+nu)) # Shear modulus of elasticity
g = 9.81 # Gravitational acceleration constant (m/s^2)

R1 = R_sph_cyl; R2 = R_catheter - R_sph_cyl

f = lambda a,X: (R1-np.sqrt(R1**2-(a*X-R2)**2) if abs(a*X-R2>=0) else 0)
fprime = lambda a,X: a*(a*X-R2)/np.sqrt(R1**2 - (a*X-R2)**2) if abs(a*X-R2>=0) else 0
    
d = lambda a: integr.quad(lambda X: fprime(a,X)/np.sqrt(1 - X**2),0,1)[0]
P = lambda a: 4*mu*a/(1-nu)*integr.quad(lambda X: X**2*fprime(a,X)/np.sqrt(1 - X**2),0,1)[0]

a = opt.fminbound(lambda X: abs(P(X)-F*g), R2, R_catheter,xtol = 1e-15)
omega_max = d(a)

uz = lambda X: 2.0/pi*(d(a)*integr.quad(lambda t: 1/np.sqrt(X**2-t**2),0,1)[0] -\
    integr.dblquad(lambda t2,t1: fprime(a,t2)/(np.sqrt(t1**2 - t2**2)*np.sqrt(X**2-t1**2)),\
    0,1,lambda y: 0,lambda y: y)[0])

x = np.linspace(a-1e-15,2*R_catheter,num = 100)/a
x = np.append(x,np.linspace(2.2*R_catheter,0.03,num = 50)/a); x = np.append(x,0.1/a)
u = np.zeros(len(x))
u_a = -d(a)+f(a,1)

for j in range(len(x)):
    if x[j] > 1:
        u[j] = -uz(x[j])
    else:
        u[j] = -d(a)+f(a,x[j]) 

slope_tol = 6; tol = 5e-5
xNew,uNew,a,u_a = curveShift(x*a,u,a,u_a,omega_max,f,tol,slope_tol)
while True:
    dist = np.sqrt(np.power(xNew-a,2) + np.power(uNew - u_a,2))
    sl = (uNew[np.argmax(dist>tol)]-u_a)/(xNew[np.argmax(dist>tol)]-a)
    if sl<=slope_tol:
        break
    else:
        xNew,uNew,a,u_a = curveShift(xNew,uNew,a,u_a,omega_max,f,tol,slope_tol)
            
# Points for elastic plane deformation with adaptive tolerance
tol = 1e-4
x = np.concatenate(([a],xNew[xNew>a]),axis=0);
y = np.concatenate(([u_a],uNew[xNew>a]),axis=0);
X = np.array([[x[0]],[y[0]]]);
for j in np.arange(1,len(x)):
    if np.sqrt(sum(np.power(X[:,-1] - [x[j],y[j]],2))) > tol:
        X = np.append(X,[[x[j]],[y[j]]],axis = 1)
        if x[j] > 4*R_catheter:
            tol = 2*tol
X = np.append(X,[[x[j]],[y[j]]],axis = 1)

# Construction of the geometry
Box = geompy.MakeBoxDXDYDZ(Lx, Ly, Lz)
interface_1 = geompy.MakeFaceHW(3*Lx, 3*Ly, 1)
geompy.TranslateDXDYDZ(interface_1, 0, 0, Lz_plastic)

# Catheter construction
CatheterOuterCylinder = geompy.MakeCylinderRH(R_catheter, 2*Lz)
geompy.TranslateDXDYDZ(CatheterOuterCylinder, 0, 0, R_sph_cyl)
CatheterInnerCylinder = geompy.MakeCylinderRH(R_catheter - R_sph_cyl, 2*Lz)
Torus_1 = geompy.MakeTorusRR(R_catheter - R_sph_cyl, R_sph_cyl)
geompy.TranslateDXDYDZ(Torus_1, 0, 0, R_sph_cyl)
CatheterBody = geompy.MakeFuseList([CatheterOuterCylinder, CatheterInnerCylinder, Torus_1], True, True)
geompy.TranslateDXDYDZ(CatheterBody, Lx/2, Ly/2, Lz_tissue - omega_max)
Cut1 = geompy.MakeFaceHW(Lx/2, Ly/2, 1)
geompy.TranslateDXDYDZ(Cut1, Lx/2, Ly/2, Lz_tissue - omega_max + L_electrode)
Cylinder_1 = geompy.MakeCylinderRH(R_thermistor, 2*Lz)
geompy.TranslateDXDYDZ(Cylinder_1, Lx/2, Ly/2, Lz_tissue - omega_max + L_electrode - L_thermistor)
Channel = geompy.MakeCylinderRH(R_tube, Lz_tissue)
Tube1 = geompy.MakeCylinderRH(R_tube_holes, L_electrode)
geompy.TranslateDXDYDZ(Tube1, 0, 0, -L_electrode/2)
geompy.Rotate(Tube1, OY, 90*math.pi/180.0)
geompy.Rotate(Tube1, OZ, 60*math.pi/180.0)
geompy.Rotate(Tube1, OZ, 30*math.pi/180.0)
Tube2 = geompy.MakeRotation(Tube1, OZ, 60*math.pi/180.0)
Tube3 = geompy.MakeRotation(Tube2, OZ, 60*math.pi/180.0)
geompy.TranslateDXDYDZ(Tube1, 0, 0, R_tube_holes)
geompy.TranslateDXDYDZ(Tube2, 0, 0, R_tube_holes)
geompy.TranslateDXDYDZ(Tube3, 0, 0, R_tube_holes)
ChannelwHoles = geompy.MakeFuseList([Channel, Tube1, Tube2, Tube3], True, True)
geompy.TranslateDXDYDZ(ChannelwHoles, Lx/2, Ly/2, Lz_tissue - omega_max + L_holes)
geompy.TranslateDXDYDZ(Channel, Lx/2, Ly/2, Lz_tissue - omega_max)
Partition_1 = geompy.MakePartition([CatheterBody], [Cut1], [], [], geompy.ShapeType["SOLID"], 0, [], 1)
[Tip] = geompy.SubShapes(Partition_1, [15])
[Hollow] = geompy.SubShapes(Partition_1, [2])
PartitionedTip = geompy.MakePartition([Tip], [ChannelwHoles], [], [], geompy.ShapeType["SOLID"], 0, [], 1)
[Inner_Channel] = geompy.SubShapes(PartitionedTip, [130])
Cut_1 = geompy.MakeCutList(Box, [Hollow], True)

# Line in the middle of box for revolutions
Vertex_1 = geompy.MakeVertex(Lx/2, Ly/2, Lz_tissue/2)
Vertex_2 = geompy.MakeVertex(Lx/2, Ly/2, Lz_tissue*2)
Line_1 = geompy.MakeLineTwoPnt(Vertex_2, Vertex_1)

# Elastic Tissue
ptArray = [];
for n in range(len(X[0,:])):
    pt = geompy.MakeVertex(X[0,n], 0.0, X[1,n])
    ptArray.append(pt)
    geompy.addToStudy(pt, "Pt_%s"%(n))
Curve_1 = geompy.MakePolyline(ptArray, False)
geompy.TranslateDXDYDZ(Curve_1, Lx/2, Ly/2, Lz_tissue)
Extrusion_1 = geompy.MakePrismDXDYDZ(Curve_1, 0, 0, -Lz_tissue/2)
Revolution_1 = geompy.MakeRevolution(Extrusion_1, Line_1, 360*math.pi/180.0)
Tissue_top = geompy.MakeCommonList([Box,Revolution_1], True)
tissue_box = geompy.MakeBoxDXDYDZ(Lx, Ly, Lz_tissue+u_a)
Cut_tissue = geompy.MakeFuseList([Tissue_top, tissue_box], True, True)
tissue_full = geompy.MakeCutList(Cut_tissue, [PartitionedTip], True)

# Final geometry
Partition_5 = geompy.MakePartition([Cut_1], [interface_1, tissue_full, Cylinder_1, PartitionedTip], [], [], geompy.ShapeType["SOLID"], 0, [], 1)
Domain = geompy.MakeCutList(Partition_5, [Inner_Channel], True)

# Add everything to salome gui
geompy.addToStudy( O, 'O' )
geompy.addToStudy( OX, 'OX' )
geompy.addToStudy( OY, 'OY' )
geompy.addToStudy( OZ, 'OZ' )
geompy.addToStudy( Box, 'Box' )
geompy.addToStudy( interface_1, 'interface_1' )
geompy.addToStudy( CatheterOuterCylinder, 'CatheterOuterCylinder' )
geompy.addToStudy( CatheterInnerCylinder, 'CatheterInnerCylinder' )
geompy.addToStudy( Torus_1, 'Torus_1' )
geompy.addToStudy( CatheterBody, 'CatheterBody' )
geompy.addToStudy( Cut1, 'Cut1' )
geompy.addToStudy( Cylinder_1, 'Cylinder_1' )
geompy.addToStudy( Channel, 'Channel' )
geompy.addToStudy( Partition_1, 'Partition_1' )
geompy.addToStudyInFather( Partition_1, Tip, 'Tip' )
geompy.addToStudyInFather( Partition_1, Hollow, 'Hollow' )
geompy.addToStudy( Tube1, 'Tube1' )
geompy.addToStudy( Tube2, 'Tube2' )
geompy.addToStudy( Tube3, 'Tube3' )
geompy.addToStudy( ChannelwHoles, 'ChannelwHoles' )
geompy.addToStudy( Partition_1, 'Partition_1' )
geompy.addToStudyInFather( Partition_1, Tip, 'Tip' )
geompy.addToStudyInFather( Partition_1, Hollow, 'Hollow' )
geompy.addToStudy( PartitionedTip, 'PartitionedTip' )
geompy.addToStudy( Inner_Channel, 'Inner_Channel' )
geompy.addToStudy( Cut_1, 'Cut_1' )
geompy.addToStudy( Curve_1, 'Curve_1' )
geompy.addToStudy( Extrusion_1, 'Extrusion_1' )
geompy.addToStudy( Revolution_1, 'Revolution_1' )
geompy.addToStudy( Tissue_top, 'Tissue_top' )
geompy.addToStudy( tissue_box, 'tissue_box' )
geompy.addToStudy( tissue_full, 'tissue_full' )
geompy.addToStudy( Partition_5, 'Partition_5' )
geompy.addToStudy( Domain, 'Domain' )

##  Define geometrical groups
geompy.SubShapeAll(Domain,geompy.ShapeType["SOLID"])
geompy.SubShapeAll(Domain,geompy.ShapeType["FACE"])

## Solids
# Board
Box_board = geompy.MakeBoxDXDYDZ(Lx, Ly, Lz_plastic)
BoardDomain = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(BoardDomain, geompy.GetShapesOnShapeIDs(Box_board, Domain, geompy.ShapeType["SOLID"], GEOM.ST_ONIN))
geompy.addToStudyInFather( Domain, BoardDomain, 'BoardDomain' )
# Tissue
Box_tissue = geompy.MakeBoxDXDYDZ(Lx, Ly, Lz_plastic)
geompy.TranslateDXDYDZ(Box_tissue,0,0,Lz_plastic)
SolidDomain = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(SolidDomain, geompy.GetShapesOnShapeIDs(Box_tissue, Domain, geompy.ShapeType["SOLID"], GEOM.ST_ONIN))
geompy.addToStudyInFather( Domain, SolidDomain, 'SolidDomain' )
# Blood 
Box_blood = geompy.MakeBoxDXDYDZ(Lx, Ly, Lz)
geompy.TranslateDXDYDZ(Box_blood,0,0,Lz_tissue-Lz_plastic/2)
Cut_blood = geompy.MakeCutList(Box_blood, [CatheterBody,SolidDomain], True)
FluidDom = geompy.MakeCommonList([Domain,Cut_blood], True)
FluidDomain = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(FluidDomain, geompy.GetSameIDs(Domain, FluidDom))
geompy.addToStudyInFather( Domain, FluidDomain, 'FluidDomain' )
# Thermistor
Thermistor = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(Thermistor, geompy.GetShapesOnShapeIDs(Cylinder_1, Domain, geompy.ShapeType["SOLID"], GEOM.ST_ONIN))
geompy.addToStudyInFather( Domain, Thermistor, 'Thermistor' )
# Electrode
Cut_elec = geompy.MakeCutList(Domain, [FluidDomain,Thermistor,SolidDomain,BoardDomain], True)
Electrode = geompy.CreateGroup(Domain, geompy.ShapeType["SOLID"])
geompy.UnionIDs(Electrode, geompy.GetShapesOnShapeIDs(Cut_elec, Domain, geompy.ShapeType["SOLID"], GEOM.ST_ONIN))
geompy.addToStudyInFather( Domain, Electrode, 'Electrode' )

## Faces
# RightFluidWall
PlaneXZ = geompy.MakePlane(O,OY,0.5)
FluidRight = geompy.MakeCommonList([PlaneXZ,FluidDomain], True)
RightFluidWall = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(RightFluidWall, geompy.GetSameIDs(Domain, FluidRight))
geompy.addToStudyInFather( Domain, RightFluidWall, 'RightFluidWall' )
# LeftFluidWall
PlaneXZ2 = geompy.MakeTranslation(PlaneXZ,0,Ly,0)
FluidLeft = geompy.MakeCommonList([PlaneXZ2,FluidDomain], True)
LeftFluidWall = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(LeftFluidWall, geompy.GetSameIDs(Domain, FluidLeft))
geompy.addToStudyInFather( Domain, LeftFluidWall, 'LeftFluidWall' )
# Inlet
PlaneYZ = geompy.MakePlane(O,OX,0.5)
FluidInlet = geompy.MakeCommonList([PlaneYZ,FluidDomain], True)
Inlet = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Inlet, geompy.GetSameIDs(Domain, FluidInlet))
geompy.addToStudyInFather( Domain, Inlet, 'Inlet' )
# OutletRevolution_1
PlaneYZ2 = geompy.MakeTranslation(PlaneYZ,Lx,0,0)
FluidOutlet = geompy.MakeCommonList([PlaneYZ2,FluidDomain], True)
Outlet = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Outlet, geompy.GetSameIDs(Domain, FluidOutlet))
geompy.addToStudyInFather( Domain, Outlet, 'Outlet' )
# WallsDownstream
Downstream = geompy.MakeCommonList([PlaneYZ2,Domain], True)
WallsDownstream = geompy.GetInPlace(Domain, Downstream)
geompy.addToStudyInFather( Domain, WallsDownstream, 'WallsDownstream' )
# Bottom
PlaneXY = geompy.MakePlane(O,OZ,0.5)
SolidBottom = geompy.MakeCommonList([PlaneXY,Domain], True)
Bottom = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Bottom, geompy.GetSameIDs(Domain, SolidBottom))
geompy.addToStudyInFather( Domain, Bottom, 'Bottom' )
# Top
PlaneXY2 = geompy.MakeTranslation(PlaneXY,0,0,Lz)
FluidTop = geompy.MakeCommonList([PlaneXY2,Domain], True)
Top = geompy.CreateGroup(Domain, geompy.ShapeType["FACE"])
geompy.UnionIDs(Top, geompy.GetSameIDs(Domain, FluidTop))
geompy.addToStudyInFather( Domain, Top, 'Top' )
# SolidWalls
OpenBox = geompy.MakeFuseList([PlaneXZ,PlaneXZ2,PlaneYZ,PlaneYZ2],True,True)
Walls1 = geompy.MakeCommonList([OpenBox,SolidDomain], True)
Walls2 = geompy.MakeCommonList([OpenBox,BoardDomain], True)
SolidWalls = geompy.GetInPlace(Domain, geompy.MakeCompound([Walls1,Walls2]))
geompy.addToStudyInFather( Domain, SolidWalls, 'SolidWalls' )
# CatheterBody
Catheter = geompy.MakeTranslation(CatheterOuterCylinder, Lx/2, Ly/2, Lz_tissue - omega_max + L_electrode - R_sph_cyl)
CathBody = geompy.GetShapesOnShapeAsCompound(Catheter, FluidDomain, geompy.ShapeType["FACE"], GEOM.ST_ON)
CatheterBody = geompy.GetInPlace(Domain, CathBody, isNewImplementation=True)
geompy.addToStudyInFather( Domain, CatheterBody, 'CatheterBody' )
# Interface
Interf = geompy.GetShapesOnShapeAsCompound(FluidDomain, SolidDomain, geompy.ShapeType["FACE"], GEOM.ST_ON)
Interface = geompy.GetInPlace(Domain, Interf, isNewImplementation=True)
geompy.addToStudyInFather( Domain, Interface, 'Interface' )
# SalineInlet
DomainCut = geompy.MakeCutList(Domain, [Electrode,Thermistor], True)
Holes = geompy.GetShapesOnShapeAsCompound(Inner_Channel, DomainCut, geompy.ShapeType["FACE"], GEOM.ST_ON)
SalineInlet = geompy.GetInPlace(Domain, Holes, isNewImplementation=True)
geompy.addToStudyInFather( Domain, SalineInlet, 'SalineInlet' )
# ThermistorUpperInterface
ThermUp = geompy.GetShapesOnShapeAsCompound(Catheter, Thermistor, geompy.ShapeType["FACE"], GEOM.ST_ON)
ThermistorUpperInterface = geompy.GetInPlace(Domain, ThermUp, isNewImplementation=True)
geompy.addToStudyInFather( Domain, ThermistorUpperInterface, 'ThermistorUpperInterface' )
# ElectrodeUpperInterface
ElecUp = geompy.GetShapesOnShapeAsCompound(Catheter, Electrode, geompy.ShapeType["FACE"], GEOM.ST_ON)
ElectrodeUpperInterface = geompy.GetInPlace(Domain, ElecUp, isNewImplementation=True)
geompy.addToStudyInFather( Domain, ElectrodeUpperInterface, 'ElectrodeUpperInterface' )
# InnerChannel
ChannelTotal = geompy.GetShapesOnShapeAsCompound(Inner_Channel, Domain, geompy.ShapeType["FACE"], GEOM.ST_ON)
InnerChannel = geompy.GetInPlace(Domain, geompy.MakeCutList(ChannelTotal, [SalineInlet], True), isNewImplementation=True)
geompy.addToStudyInFather( Domain, InnerChannel, 'InnerChannel' )
# Elec_tissue
Shared_Face_1 = geompy.GetSharedShapesMulti([SolidDomain, Electrode], geompy.ShapeType["FACE"], False)
Elec_tissue = geompy.GetInPlace(Domain, geompy.MakeCompound(Shared_Face_1), isNewImplementation=True)
geompy.addToStudyInFather( Domain, Elec_tissue, 'Elec_tissue' )


# Mesh objects
smesh = smeshBuilder.New(theStudy)
Mesh = smesh.Mesh(Domain)
NETGEN_2D = Mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D)
GeneralParameters2D = NETGEN_2D.Parameters()
GeneralParameters2D.SetMaxSize( 0.002 )
GeneralParameters2D.SetSecondOrder( 0 )
GeneralParameters2D.SetOptimize( 1 )
GeneralParameters2D.SetFineness( 3 )
GeneralParameters2D.SetMinSize( 0 )
GeneralParameters2D.SetUseSurfaceCurvature( 1 )
GeneralParameters2D.SetFuseEdges( 1 )
GeneralParameters2D.SetQuadAllowed( 0 )
NETGEN_3D = Mesh.Tetrahedron()
GeneralParameters3D = NETGEN_3D.Parameters()
GeneralParameters3D.SetMaxSize( 0.001 )
GeneralParameters3D.SetOptimize( 1 )
GeneralParameters3D.SetFineness( 3 )
GeneralParameters3D.SetMinSize( 0 )
GeneralParameters3D.SetUseSurfaceCurvature( 0 )
GeneralParameters3D.SetFuseEdges( 0 )
GeneralParameters3D.SetSecondOrder( 0 )
GeneralParameters3D.SetQuadAllowed( 0 )
NETGEN_2D_1 = Mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D,geom=SalineInlet)
SalineInflowMesh = NETGEN_2D_1.GetSubMesh()
SalineInflowParameters = NETGEN_2D_1.Parameters()
SalineInflowParameters.SetMaxSize( 0.0001 )
SalineInflowParameters.SetSecondOrder( 0 )
SalineInflowParameters.SetOptimize( 1 )
SalineInflowParameters.SetFineness( 2 )
SalineInflowParameters.SetMinSize( 0 )
SalineInflowParameters.SetUseSurfaceCurvature( 1 )
SalineInflowParameters.SetFuseEdges( 1 )
SalineInflowParameters.SetQuadAllowed( 0 )
NETGEN_2D_2 = Mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D,geom=FluidDomain)
FluidMesh = NETGEN_2D_2.GetSubMesh()
FluidParameters2D = NETGEN_2D_2.Parameters()
FluidParameters2D.SetMaxSize( 0.001 )
FluidParameters2D.SetSecondOrder( 0 )
FluidParameters2D.SetOptimize( 1 )
FluidParameters2D.SetFineness( 3 )
FluidParameters2D.SetMinSize( 0 )
FluidParameters2D.SetUseSurfaceCurvature( 1 )
FluidParameters2D.SetFuseEdges( 1 )
FluidParameters2D.SetQuadAllowed( 0 )
NETGEN_3D_1 = Mesh.Tetrahedron(geom=FluidDomain)
FluidParameters3D = NETGEN_3D_1.Parameters()
FluidParameters3D.SetMaxSize( 0.00075 )
FluidParameters3D.SetSecondOrder( 0 )
FluidParameters3D.SetOptimize( 1 )
FluidParameters3D.SetFineness( 3 )
FluidParameters3D.SetMinSize( 0 )
FluidParameters3D.SetQuadAllowed( 0 )
FluidParameters3D.SetUseSurfaceCurvature( 63 )
FluidParameters3D.SetFuseEdges( 176 )
NETGEN_2D_3 = Mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D,geom=Electrode)
ElectrodeMesh = NETGEN_2D_3.GetSubMesh()
status = Mesh.AddHypothesis(FluidParameters2D,Electrode)
NETGEN_3D_2 = Mesh.Tetrahedron(geom=Electrode)
status = Mesh.AddHypothesis(GeneralParameters3D,Electrode)
NETGEN_2D_4 = Mesh.Triangle(algo=smeshBuilder.NETGEN_1D2D,geom=Elec_tissue)
Elec_tissueMesh = NETGEN_2D_4.GetSubMesh()
status = Mesh.AddHypothesis(SalineInflowParameters,Elec_tissue)
isDone = Mesh.SetMeshOrder( [ [ Elec_tissueMesh, SalineInflowMesh, ElectrodeMesh, FluidMesh ] ])
RightFluidWall_1 = Mesh.GroupOnGeom(RightFluidWall,'RightFluidWall',SMESH.FACE)
Outlet_1 = Mesh.GroupOnGeom(Outlet,'Outlet',SMESH.FACE)
Inlet_1 = Mesh.GroupOnGeom(Inlet,'Inlet',SMESH.FACE)
Top_1 = Mesh.GroupOnGeom(Top,'Top',SMESH.FACE)
LeftFluidWall_1 = Mesh.GroupOnGeom(LeftFluidWall,'LeftFluidWall',SMESH.FACE)
Bottom_1 = Mesh.GroupOnGeom(Bottom,'Bottom',SMESH.FACE)
SolidWalls_1 = Mesh.GroupOnGeom(SolidWalls,'SolidWalls',SMESH.FACE)
CatheterBody_1 = Mesh.GroupOnGeom(CatheterBody,'CatheterBody',SMESH.FACE)
Interface_1 = Mesh.GroupOnGeom(Interface,'Interface',SMESH.FACE)
SalineInlet_1 = Mesh.GroupOnGeom(SalineInlet,'SalineInlet',SMESH.FACE)
WallsDownstream_1 = Mesh.GroupOnGeom(WallsDownstream,'WallsDownstream',SMESH.FACE)
ThermistorUpperInterface_1 = Mesh.GroupOnGeom(ThermistorUpperInterface,'ThermistorUpperInterface',SMESH.FACE)
ElectrodeUpperInterface_1 = Mesh.GroupOnGeom(ElectrodeUpperInterface,'ElectrodeUpperInterface',SMESH.FACE)
InnerChannel_1 = Mesh.GroupOnGeom(InnerChannel,'InnerChannel',SMESH.FACE)
FluidDomain_1 = Mesh.GroupOnGeom(FluidDomain,'FluidDomain',SMESH.VOLUME)
SolidDomain_1 = Mesh.GroupOnGeom(SolidDomain,'SolidDomain',SMESH.VOLUME)
SolidDomain_1 = Mesh.GroupOnGeom(BoardDomain,'BoardDomain',SMESH.VOLUME)
Electrode_1 = Mesh.GroupOnGeom(Electrode,'Electrode',SMESH.VOLUME)
Thermistor_1 = Mesh.GroupOnGeom(Thermistor,'Thermistor',SMESH.VOLUME)
isDone = Mesh.Compute()
[ RightFluidWall_1, Outlet_1, Inlet_1, Top_1, LeftFluidWall_1, Bottom_1, \
SolidWalls_1, CatheterBody_1, Interface_1, SalineInlet_1, WallsDownstream_1, \
ThermistorUpperInterface_1, ElectrodeUpperInterface_1, InnerChannel_1, FluidDomain_1, \
SolidDomain_1, BoardDomain_1, Electrode_1, Thermistor_1 ] = Mesh.GetGroups()


# Set names of Mesh objects
smesh.SetName(NETGEN_2D.GetAlgorithm(), 'NETGEN_2D')
smesh.SetName(CatheterBody_1, 'CatheterBody')
smesh.SetName(Interface_1, 'Interface')
smesh.SetName(NETGEN_3D.GetAlgorithm(), 'NETGEN_3D')
smesh.SetName(GeneralParameters2D, 'GeneralParameters2D')
smesh.SetName(SalineInflowParameters, 'SalineInflowParameters')
smesh.SetName(GeneralParameters3D, 'GeneralParameters3D')
smesh.SetName(FluidParameters2D, 'FluidParameters2D')
smesh.SetName(RightFluidWall_1, 'RightFluidWall')
smesh.SetName(FluidParameters3D, 'FluidParameters3D')
smesh.SetName(Outlet_1, 'Outlet')
smesh.SetName(Inlet_1, 'Inlet')
smesh.SetName(Top_1, 'Top')
smesh.SetName(LeftFluidWall_1, 'LeftFluidWall')
smesh.SetName(SalineInflowMesh, 'SalineInflowMesh')
smesh.SetName(Bottom_1, 'Bottom')
smesh.SetName(SolidWalls_1, 'SolidWalls')
smesh.SetName(Mesh.GetMesh(), 'Mesh')
smesh.SetName(WallsDownstream_1, 'WallsDownstream')
smesh.SetName(SalineInlet_1, 'SalineInlet')
smesh.SetName(ElectrodeUpperInterface_1, 'ElectrodeUpperInterface')
smesh.SetName(Thermistor_1, 'Thermistor')
smesh.SetName(ThermistorUpperInterface_1, 'ThermistorUpperInterface')
smesh.SetName(Electrode_1, 'Electrode')
smesh.SetName(SolidDomain_1, 'SolidDomain')
smesh.SetName(InnerChannel_1, 'InnerChannel')
smesh.SetName(FluidDomain_1, 'FluidDomain')
smesh.SetName(BoardDomain_1, 'BoardDomain')
smesh.SetName(FluidMesh, 'FluidMesh')
smesh.SetName(ElectrodeMesh, 'ElectrodeMesh')
smesh.SetName(Elec_tissueMesh, 'Elec_tissueMesh')


dest = os.getcwd()
try:
  Mesh.ExportUNV( dest+'/WallsDownstream_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  WallsDownstream_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Top_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Top_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/ThermistorUpperInterface_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  ThermistorUpperInterface_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/SolidWalls_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  SolidWalls_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/SalineInlet_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  SalineInlet_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/RightFluidWall_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  RightFluidWall_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Outlet_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Outlet_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/LeftFluidWall_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  LeftFluidWall_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/InnerChannel_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  InnerChannel_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Interface_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Interface_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/ElectrodeUpperInterface_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  ElectrodeUpperInterface_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/CatheterBody_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  CatheterBody_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Bottom_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Bottom_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Inlet_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Inlet_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Thermistor_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Thermistor_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/SolidDomain_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  SolidDomain_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/BoardDomain_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  BoardDomain_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Electrode_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  Electrode_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/FluidDomain_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv',  FluidDomain_1)
except:
  print 'ExportPartToUNV() failed. Invalid file name?'
try:
  Mesh.ExportUNV( dest+'/Domain_elastic_cyl_'+str(int(F*1e3))+'g'+str(int(E*1e-3))+'kPa.unv')
except:
  print 'ExportTounvX() failed. Invalid file name?'
