Cells in a cylindrical geometry
This demo shows cell arangement in a cylindrical geometry. The representation deals with cells in the x-y plane at a given z slice.
Import modules and load parfis lib
[1]:
import parfis as pfs
from parfis import Parfis
import matplotlib.pyplot as plt
import matplotlib.patches as patches
Parfis.load_lib()
Successfully loaded lib file: libparfis64.so
Create simulation data
[2]:
id = Parfis.newParfis()
Parfis.loadCfgData(id)
Parfis.loadSimData(id)
success = Parfis.runCommandChain(id, "create")
print("success") if success == 0 else print("fail")
success
Get pointers to configuration and simulation data
[3]:
Parfis.setPyCfgData(id)
ptrCfgData = Parfis.getPyCfgData(id)
Parfis.setPySimData(id)
ptrSimData = Parfis.getPySimData(id)
Get geometry and cell measures
[4]:
cellSize =ptrCfgData.cellSize[0]
geoSize =ptrCfgData.geometrySize[0]
cellCount =ptrCfgData.cellCount[0]
print(f"cellSize = {cellSize} meters")
print(f"geoSize = {geoSize} meters")
print(f"cellCount = {cellCount}")
cellSize = {x: 0.001, y: 0.001, z: 0.001} meters
geoSize = {x: 0.02, y: 0.02, z: 0.4} meters
cellCount = {x: 20, y: 20, z: 400}
Get all cells in the x-y plane at z=10
Every absolute id by value is Const.noCellId or a certain id that corresponds to the position in the vector of existing cells. We gather two vectors, one with absolute cell ids and one with real cell ids.
Absolute id values \(A_{id}\) are calculated according to the following formula:
where \(C_x\), \(C_y\) and \(C_z\) are number of cells in the given direction and \(x\), \(y\) and \(z\) are the components of the cell position vector. If geometry is cylindrical, then there are cells that lie outside the geometry. Cell outside of the geometry have real id that is equal to Const.noCellIc (the maximum uint32 number).
[5]:
cellPos = pfs.Vec3DClass(pfs.Type.cellPos_t)(z=10)
absCellId = [] # Absolute cell ids
realCellId = [] # Real cell ids
boundCellId = [] # Cell ids that lie on the boundary
print(f"Cell pos Abs. id Rel. id Cell pos from data")
print(f"-----------------------------------------------------------------")
for i in range(cellCount.x):
cellPos.x = i
for j in range(cellCount.y):
cellPos.y = j
absCellId.append(pfs.getAbsoluteCellId(cellCount, cellPos))
cellId =ptrSimData.cellIdVec.ptr[absCellId[-1]]
cellPosFromData = None
if cellId != pfs.Const.noCellId:
realCellId.append(cellId)
cellPosFromData =ptrSimData.cellVec.ptr[cellId].pos
nodeFlag =ptrSimData.nodeFlagVec.ptr[cellId]
if nodeFlag != 0xFF:
boundCellId.append(cellId)
if i==0 and j < 10:
print(f"{cellPos} {absCellId[-1]:7} {cellId:10} {cellPosFromData}")
Cell pos Abs. id Rel. id Cell pos from data
-----------------------------------------------------------------
{x: 0, y: 0, z: 10} 10 4294967295 None
{x: 0, y: 1, z: 10} 410 4294967295 None
{x: 0, y: 2, z: 10} 810 4294967295 None
{x: 0, y: 3, z: 10} 1210 4294967295 None
{x: 0, y: 4, z: 10} 1610 4294967295 None
{x: 0, y: 5, z: 10} 2010 10 {x: 0, y: 5, z: 10}
{x: 0, y: 6, z: 10} 2410 410 {x: 0, y: 6, z: 10}
{x: 0, y: 7, z: 10} 2810 810 {x: 0, y: 7, z: 10}
{x: 0, y: 8, z: 10} 3210 1210 {x: 0, y: 8, z: 10}
{x: 0, y: 9, z: 10} 3610 1610 {x: 0, y: 9, z: 10}
Draw all cells that exist in the simulation.
Draw only cells that have relative id, because they exist in the cellVec vector.
[6]:
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_title("Existing cells in the simulation")
ax.plot()
for cellId in realCellId:
pos =ptrSimData.cellVec.ptr[cellId].pos
ax.add_patch(
patches.Rectangle(
(cellSize.x*pos.x, cellSize.y*pos.y),
cellSize.x, cellSize.y,
edgecolor = 'blue',
alpha = 0.5,
fill=False
)
)
# Draw the geometry bound
ax.add_patch(
patches.Circle(
(geoSize.x*0.5, geoSize.y*0.5),
radius=geoSize.x*0.5,
lw = 2,
edgecolor = 'red',
alpha = 0.5,
fill=False
)
)
ax.set_xlim(0 - cellSize.x, geoSize.x + cellSize.x)
ax.set_ylim(0 - cellSize.y, geoSize.y + cellSize.y)
ax.set_aspect('equal')
plt.show()
Represent cells that are not fully inside the geometry
These cells have nodeFlag that is different than uint8 max, which means that the cell has some nodes outside the boundary. Cells on the boundary are important for the interaction of particles and walls since for particles in these cells that interaction can occur.
[7]:
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_title("Cells on the boundary")
ax.plot()
for cellId in boundCellId:
pos =ptrSimData.cellVec.ptr[cellId].pos
ax.add_patch(
patches.Rectangle(
(cellSize.x*pos.x, cellSize.y*pos.y),
cellSize.x, cellSize.y,
edgecolor = 'blue',
alpha = 0.5,
fill=False
)
)
# Draw the geometry bound
ax.add_patch(
patches.Circle(
(geoSize.x*0.5, geoSize.y*0.5),
radius=geoSize.x*0.5,
lw = 2,
edgecolor = 'red',
alpha = 0.5,
fill=False
)
)
ax.set_xlim(0 - cellSize.x, geoSize.x + cellSize.x)
ax.set_ylim(0 - cellSize.y, geoSize.y + cellSize.y)
ax.set_aspect('equal')
plt.show()
Represent cells that have neighbours that are not fully inside the geometry
If a particle finds itself in one of these cells then it can reach the boundary in one timestep. These cells are important since for them additional conditional must be checked (did the particle from the cell crossed a boundary). Since these cells are important, ids of the cells are stored in a separate vector, the “B” vector of cell id.
[8]:
fig, ax = plt.subplots(figsize=(10, 10))
ax.set_title("Cells on the x-y boundary")
ax.plot()
for i in range(ptrSimData.cellIdBVec.size):
cellId = ptrSimData.cellIdBVec.ptr[i]
pos = ptrSimData.cellVec.ptr[cellId].pos
if pos.z != 10:
continue
ax.add_patch(
patches.Rectangle(
(cellSize.x*pos.x, cellSize.y*pos.y),
cellSize.x, cellSize.y,
edgecolor = 'blue',
alpha = 0.5,
fill=False
)
)
# Draw the geometry bound
ax.add_patch(
patches.Circle(
(geoSize.x*0.5, geoSize.y*0.5),
radius=geoSize.x*0.5,
lw = 2,
edgecolor = 'red',
alpha = 0.5,
fill=False
)
)
ax.set_xlim(0 - cellSize.x, geoSize.x + cellSize.x)
ax.set_ylim(0 - cellSize.y, geoSize.y + cellSize.y)
ax.set_aspect('equal')
plt.show()
Let’s check the x-slice and the cells in the vicinity of the z-boundary
[9]:
fig, ax = plt.subplots(1, 2, figsize=(20, 10))
ax[0].set_title("Cells on the -z boundary")
ax[1].set_title("Cells on the +z boundary")
ax[0].plot()
ax[1].plot()
for i in range(ptrSimData.cellIdBVec.size):
cellId = ptrSimData.cellIdBVec.ptr[i]
pos = ptrSimData.cellVec.ptr[cellId].pos
if pos.x != cellCount.x//2:
continue
if pos.z < 5:
ax[0].add_patch(
patches.Rectangle(
(cellSize.z*pos.z, cellSize.y*pos.y),
cellSize.y, cellSize.z,
edgecolor = 'blue',
alpha = 0.5,
fill=False
)
)
elif pos.z > cellCount.z-6:
ax[1].add_patch(
patches.Rectangle(
(cellSize.z*pos.z, cellSize.y*pos.y),
cellSize.y, cellSize.z,
edgecolor = 'blue',
alpha = 0.5,
fill=False
)
)
# Draw the geometry bound
ax[0].add_patch(
patches.Rectangle(
(0.0, 0.0),
geoSize.z, geoSize.y,
lw = 2,
edgecolor = 'red',
alpha = 0.5,
fill=False
)
)
# Draw the geometry bound
ax[1].add_patch(
patches.Rectangle(
(0.0, 0.0),
geoSize.z, geoSize.y,
lw = 2,
edgecolor = 'red',
alpha = 0.5,
fill=False
)
)
ax[0].set_xlim(0 - cellSize.z, cellSize.z*5)
ax[0].set_ylim(0 - cellSize.y, geoSize.y + cellSize.y)
ax[0].set_aspect('equal')
ax[1].set_xlim((cellCount.z - 5)*cellSize.z, (cellCount.z + 1)*cellSize.z)
ax[1].set_ylim(0 - cellSize.y, geoSize.y + cellSize.y)
ax[1].set_aspect('equal')
plt.show()
