Simulating lid driven cavity flow

import math

import vtk
from ENigMA import ENigMA

# create a rendering window and renderer
ren = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)

renWin.SetSize(800, 600)

# create a renderwindowinteractor
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

# enable user interface interactor
iren.Initialize()

vertex1 = ENigMA.CGeoCoordinate(+0.00, +0.00, -0.1)
vertex2 = ENigMA.CGeoCoordinate(+1.00, +0.00, -0.1)
vertex3 = ENigMA.CGeoCoordinate(+1.00, +1.00, -0.1)
vertex4 = ENigMA.CGeoCoordinate(+0.00, +1.00, -0.1)
vertex5 = ENigMA.CGeoCoordinate(+0.00, +0.00, +0.1)
vertex6 = ENigMA.CGeoCoordinate(+1.00, +0.00, +0.1)
vertex7 = ENigMA.CGeoCoordinate(+1.00, +1.00, +0.1)
vertex8 = ENigMA.CGeoCoordinate(+0.00, +1.00, +0.1)

hexahedron = ENigMA.CGeoHexahedron()

hexahedron.addVertex(vertex1)
hexahedron.addVertex(vertex2)
hexahedron.addVertex(vertex3)
hexahedron.addVertex(vertex4)
hexahedron.addVertex(vertex5)
hexahedron.addVertex(vertex6)
hexahedron.addVertex(vertex7)
hexahedron.addVertex(vertex8)

mesher = ENigMA.CMshBasicMesher()

mesher.generate(hexahedron, 40, 40, 1, False)

mesh = mesher.mesh()

mesh.generateFaces(1E-3)

unstructured_grid = vtk.vtkUnstructuredGrid();

points = vtk.vtkPoints()

for i in range(0, mesh.nbNodes()):

    node_id = mesh.nodeId(i)

    node = mesh.node(node_id)

    points.InsertNextPoint(node.x(), node.y(), node.z())

unstructured_grid.SetPoints(points)

for i in range(0, mesh.nbElements()):

    element_id = mesh.elementId(i)

    element = mesh.element(element_id)

    if (element.elementType() == ENigMA.ET_HEXAHEDRON):

        hexahedron = vtk.vtkHexahedron()
        hexahedron.GetPointIds().SetId(0, mesh.nodeIndex(element.nodeId(0)))
        hexahedron.GetPointIds().SetId(1, mesh.nodeIndex(element.nodeId(1)))
        hexahedron.GetPointIds().SetId(2, mesh.nodeIndex(element.nodeId(2)))
        hexahedron.GetPointIds().SetId(3, mesh.nodeIndex(element.nodeId(3)))
        hexahedron.GetPointIds().SetId(4, mesh.nodeIndex(element.nodeId(4)))
        hexahedron.GetPointIds().SetId(5, mesh.nodeIndex(element.nodeId(5)))
        hexahedron.GetPointIds().SetId(6, mesh.nodeIndex(element.nodeId(6)))
        hexahedron.GetPointIds().SetId(7, mesh.nodeIndex(element.nodeId(7)))

        unstructured_grid.InsertNextCell(hexahedron.GetCellType(), hexahedron.GetPointIds())

mesh.calculateFaceCentroid()
mesh.calculateElementCentroid()

fvm_mesh = ENigMA.CFvmMesh(mesh)

U = 1.0    # lid velocity
mu = 0.001 # dynamic viscosity
rho = 1.0  # density

pisoSolver = ENigMA.CFvmPisoSolver(fvm_mesh)

pisoSolver.setGravity(0.0, 0.0, 0.0)
pisoSolver.setMaterialProperties(rho, mu)

faceIds = ENigMA.StdVectorInt()

faceIds.clear()
for i in range(0, fvm_mesh.nbFaces()):   

    faceId = fvm_mesh.faceId(i)
    face = fvm_mesh.face(faceId)

    face.calculateCentroid()

    if (face.centroid().x() == 0.0 or
        face.centroid().x() == 1.0 or
        face.centroid().y() == 0.0 or
        face.centroid().y() == 1.0):
        faceIds.push_back(faceId)

pisoSolver.setBoundaryVelocity(faceIds, ENigMA.BT_WALL_NO_SLIP, 0.0, 0.0, 0.0)
pisoSolver.setBoundaryPressure(faceIds, ENigMA.BT_WALL_NO_SLIP, 0.0)

faceIds.clear()
for i in range(0, fvm_mesh.nbFaces()):

    faceId = fvm_mesh.faceId(i)
    face = fvm_mesh.face(faceId)

    face.calculateCentroid()

    if (face.centroid().y() == 1.0):
        faceIds.push_back(faceId)

pisoSolver.setBoundaryVelocity(faceIds, ENigMA.BT_WALL_NO_SLIP, U, 0.0, 0.0)

scalars = vtk.vtkFloatArray()
scalars.SetNumberOfTuples(fvm_mesh.nbControlVolumes())
geometry_filter = vtk.vtkGeometryFilter()
cell_data_to_point_data_filter1 = vtk.vtkCellDataToPointData()
bounded_filter = vtk.vtkBandedPolyDataContourFilter()
banded_mapper = vtk.vtkPolyDataMapper()

lookup_table = vtk.vtkLookupTable()
lookup_table.SetNumberOfColors(256)
lookup_table.SetHueRange(0.667, 0.0)
lookup_table.Build()

# add actor to the renderer
mesh_actor = vtk.vtkActor()
ren.AddActor(mesh_actor)

scalar_bar_actor = vtk.vtkScalarBarActor()
ren.AddActor(scalar_bar_actor)

# Courant < 1
dt = U / 40			

iter = 1000

# Flow in a rectangle
for it in range(0, iter):
    print('Iter = ' + str(it + 1) + ' of ' + str(iter))
    pisoSolver.iterate(dt)
    
    for i in range(0, fvm_mesh.nbControlVolumes()):

        control_volume_id = fvm_mesh.controlVolumeId(i)

        u = pisoSolver.u(control_volume_id)
        v = pisoSolver.v(control_volume_id)

        vel = math.sqrt(u * u + v * v)

        scalars.SetTuple1(control_volume_id, vel)

    scalars.Modified()

    unstructured_grid.GetCellData().SetScalars(scalars)
    geometry_filter.SetInputData(unstructured_grid)
    geometry_filter.Update()

    cell_data_to_point_data_filter1.SetInputConnection(geometry_filter.GetOutputPort())
    cell_data_to_point_data_filter1.Update()

    bounded_filter.SetInputData(cell_data_to_point_data_filter1.GetOutput())
    bounded_filter.GenerateValues(24, scalars.GetRange()[0], scalars.GetRange()[1])

    banded_mapper.SetInputConnection(bounded_filter.GetOutputPort())
    banded_mapper.SetScalarModeToUsePointData()
    banded_mapper.SetScalarRange(scalars.GetRange()[0], scalars.GetRange()[1])
    banded_mapper.SetLookupTable(lookup_table)

    mesh_actor.SetMapper(banded_mapper)

    scalar_bar_actor.SetLookupTable(lookup_table)
    scalar_bar_actor.SetTitle("Velocity")
    scalar_bar_actor.SetNumberOfLabels(6)

    renWin.Render()

    if (it == 0):
        ren.ResetCamera()

iren.Start()