1.9. 3D Orthogonal Grid Generation with multiple block IDs

This is very similar to the previous 3D Orthogonal Grid tutorial.

However, we assign several block IDs.

To run the code, simply type: jupyter nbconvert --to python --execute <basename>.ipynb.

To convert it to a python file (named <basename>.py), simply type: jupyter nbconvert --to python <basename>.ipynb

To run the python file from the terminal, using N processes, simply type: mpiexec -n <N> python <basename>.py

import os
import sys

sys.path.append("../../..")

from pyopensn.mesh import OrthogonalMeshGenerator, PETScGraphPartitioner
from pyopensn.context import UseColor, Finalize

UseColor(False)

1.9.1. List of Nodes, Mesh GGeneration and Partition

We create a list of geometry nodes that are evenly spaced from 0 to +L. These nodes are employed to define block IDs interfaces.

We create a list of mesh nodes that are evenly spaced from 0 to +L. These nodes are the standard ones, used to define the vertices of our computational grid.

length = 10.
n_cells_geo = 10
dx_geo = length / n_cells_geo
nodes_geo = [i * dx_geo for i in range(n_cells_geo + 1)]

n_refinement = 5
n_cells = n_cells_geo * n_refinement
dx = length / n_cells
nodes_msh = [i * dx for i in range(n_cells + 1)]

# Setup mesh
meshgen = OrthogonalMeshGenerator(
    node_sets=[nodes_msh, nodes_msh, nodes_msh],
    partitioner=PETScGraphPartitioner(type='parmetis'),
)
grid = meshgen.Execute()

1.9.2. Material IDs

We set all block IDs to a uniform initial value 0 for each cell in the spatial domain.

We create a function factory to create a material ID function for a given (i, j, k) triplet.

• It returns a two-argument function (mat_id) that checks whether a point pt lies

  • between nodes[i] and nodes[i+1] in the x-direction,

  • between nodes[j] and nodes[j+1] in the y-direction,

  • and between nodes[k] and nodes[k+1] in the z-direction.

• If the condition is met, it returns the specified material_id; otherwise, it leaves the current block ID unchanged.

grid.SetUniformBlockID(0)

def make_mat_id_function(i, j, k, material_id):
    """
    Returns a function that sets the block to 'material_id' if the point 'pt' lies within:
    x between nodes[i] and nodes[i+1],
    y between nodes[j] and nodes[j+1],
    z between nodes[k] and nodes[k+1].
    """
    def mat_id(pt, cur_id):
        if nodes_geo[i] < pt.x < nodes_geo[i+1] and \
                nodes_geo[j] < pt.y < nodes_geo[j+1] and \
                nodes_geo[k] < pt.z < nodes_geo[k+1]:
            return material_id
        return cur_id
    return mat_id

1.9.2.1. Example: List of target blocks, where each entry is a tuple: ((i, j, k), material_id)

target_blocks = [
    ((4, 4, 7), 9),  # Block defined by nodes[4] to nodes[5] in all dimensions gets material_id 9.
    ((4, 7, 7), 55)   # Another block for demonstration.
    ]

# Create a list of material ID functions for each target block.
mat_id_functions = [make_mat_id_function(i, j, k, material_id)
                       for (i, j, k), material_id in target_blocks]

# Apply each material ID function to the grid.
for func in mat_id_functions:
    grid.SetBlockIDFromFunction(func)

1.9.3. Export the mesh

We export to vtu format. The resulting mesh partition is shown below

grid.ExportToPVTU("ortho_3D_BlockIDs")

1.9.4. Finalize (for Jupyter Notebook only)

In Python script mode, PyOpenSn automatically handles environment termination. However, this automatic finalization does not occur when running in a Jupyter notebook, so explicit finalization of the environment at the end of the notebook is required. Do not call the finalization in Python script mode, or in console mode.

Note that PyOpenSn’s finalization must be called before MPI’s finalization.

from IPython import get_ipython

def finalize_env():
    Finalize()
    MPI.Finalize()

ipython_instance = get_ipython()
if ipython_instance is not None:
    ipython_instance.events.register("post_execute", finalize_env)