Block to Block Connectivity

Computational Grids are often divided into blocks for easier gridding. It is important to know how these blocks connect with each other in order to model the transfer of information. Computational Fluid Dynamics (CFD) requires the solver to know the “connectivity” of these blocks.

The plot3D python library has a function for determining the connectivity. This is found by taking the combination of each block with another block and scanning the nodes on the exterior faces of block 1 with the exterior face of block 2. If there is a match then the two faces are connected. The code can search for full face matches or partial matching. Full face matching is much faster and should be used if your grid has a definite match.

Below is an example of a turbine domain with two blocks: HGrid in (red) and ogrid in (white). The two share a connected face with partial matching. There’s also the o-mesh which within itself has a connected face.

To find connectivity simply call the function connectivity and pass in the blocks. This will return two things: face_matches and outer faces. Face matches are faces that are connected between the blocks. Outer faces are the remaining outer faces that aren’t connected to anything. This can take a long time to run for larger meshes.

blocks = read_plot3D('PahtCascade.xyz', binary = True, big_endian=True)
# Block 1 is the blade O-Mesh k=0
# outer_faces, _ = get_outer_faces(blocks[0]) # lets check
face_matches, outer_faces_formatted = connectivity(blocks)

This is an example of the output. The figure below shows the block 1 and 2 matching faces found in the variable face_matches.

This is the matching face within the omesh

Face Orientation: Cross-Plane Connectivity

When two connected faces share the same constant axis (e.g., both K-constant), the two varying axes are identical and only 4 direction combinations exist. These are fully encodable by the lb/ub diagonal convention.

When faces lie on different constant planes (e.g., a K-face connects to a J-face), the varying axes differ. Face A varies in (i, j) while Face B varies in (i, k). This introduces an additional degree of freedom — which axis maps to which:

• No swap: A(i,j) maps to B(i,k) — 4 direction combos (encodable by lb/ub)

• Swapped: A(i,j) maps to B(k,i) — 4 direction combos (need orientation)

4 × 2 = 8 total permutations. lb/ub only encodes 4.

Permutation Matrix System

Both the Python and Rust implementations encode all 8 orientations as a 3-bit permutation_index into an array of 2×2 signed matrices (PERMUTATION_MATRICES):

permutation_index = u_reversed | (v_reversed << 1) | (swapped << 2)

The 8 Permutation Matrices

Index	Binary	u_rev	v_rev	swap	Matrix	Effect
0	000	no	no	no	[[ 1, 0],[ 0, 1]]	identity
1	001	yes	no	no	[[-1, 0],[ 0, 1]]	flip u
2	010	no	yes	no	[[ 1, 0],[ 0,-1]]	flip v
3	011	yes	yes	no	[[-1, 0],[ 0,-1]]	flip both
4	100	no	no	yes	[[ 0, 1],[ 1, 0]]	transpose
5	101	yes	no	yes	[[ 0,-1],[ 1, 0]]	transpose + flip u
6	110	no	yes	yes	[[ 0, 1],[-1, 0]]	transpose + flip v
7	111	yes	yes	yes	[[ 0,-1],[-1, 0]]	transpose + both

In Python, the PERMUTATION_MATRICES constant is defined in connectivity.py and _orient_vec_to_permutation() converts the legacy orientation vector to a permutation_index.

In Rust, the PERMUTATION_MATRICES constant is in face_record.rs and each FaceMatch carries an Orientation { permutation_index, plane } set by verify_connectivity.

The u and v names are abstract — they map to concrete i/j/k axes depending on which axis is constant:

u/v to axis mapping

Constant axis	u (outer loop)	v (inner loop)
I-constant	j	k
J-constant	i	k
K-constant	i	j

So u_reversed: true means the outer-loop axis runs opposite direction on block2 vs block1, and v_reversed: true means the inner-loop axis runs opposite.

Legacy Orientation Vector

The older Python _compute_orientation() in connectivity.py produces an orientation vector that maps each face1 axis to a face2 axis:

# orientation = [d1_maps_to, d2_maps_to, d3_maps_to]
# 1-indexed: 1=I, 2=J, 3=K
#
# Example: orientation = [2, 1, 3]
#   face1 I-axis -> face2 J-axis
#   face1 J-axis -> face2 I-axis
#   face1 K-axis -> face2 K-axis

Direction (forward/reverse) is encoded in the lb/ub values, not in the orientation vector. The _orient_vec_to_permutation() function converts this vector to a permutation_index for the unified system. The verification modules (verify_connectivity / verify_periodicity) test all 8 permutations to confirm the match.

For a detailed analysis with diagrams, see the Root Cause Analysis document.

Directed Diagonal for GHT_CONN Export

The GlennHT connectivity file (.ght_conn) uses a directed diagonal convention where each face is specified by two corners (IMIN, JMIN, KMIN) and (IMAX, JMAX, KMAX). For cross-plane matches, the traversal direction is encoded by allowing the “max” corner to be less than the “min” corner on reversed axes. For example:

1   1   1   1   25  409   1
2 409   1   1    1    1  25

Here block 2’s i-axis runs 409 → 1 (reversed), encoding the cross-plane orientation without a separate permutation matrix.

When connectivity is computed in Python (via connectivity_fast), the lb/ub values already encode the directed diagonal from the point-match traversal order.

When connectivity is computed in Rust (via the connectivity-finder binary in grid-packed), the JSON output uses ascending lo/hi bounds with a permutation_index (0-7). Use reconstruct_directed_diagonal() to convert these to directed lb/ub before exporting:

from plot3d.connectivity import reconstruct_directed_diagonal

# face_match has ascending lb/ub + permutation_index from Rust JSON
directed_match = reconstruct_directed_diagonal(face_match)
# directed_match now has directed lb/ub and permutation_index = -1

The export_to_glennht_conn() function applies this reconstruction automatically, so callers do not need to call it explicitly.

Plotting Connectivity using Paraview

This example shows how you can take a mesh created Numeca Autogrid and plot the connecitivity using paraview. With this information, anyone should be able to comprehend the format and export it to whatever solver. This code has been tested using paraview 5.9. Modifications may need to be added for Paraview 5.10

Note

Important Note: This is a 2 step process because: Paraview has it’s own python environment called pvpython which it uses to automate plotting. You do not have full access to this library, so pip installs won’t work. This means you cannot directly call the plot3d library from paraview.

We need to split this into 2 codes.

Part 1 One code runs on python from your own environment.

Part 2 the other code runs from paraview’s python environment.

The goal if this code is to create a pickle file that is read in part 2. This pickle file contains the matching faces and outer faces.

from itertools import combinations
import os, sys
sys.path.insert(0,'../../')
sys.path.insert(1,"../Cascade")
from plot3d import write_plot3D, read_plot3D, find_periodicity
from plot3d import find_matching_blocks, get_outer_faces, connectivity
from glennht_con import export_to_glennht_conn
import pickle

# Convert to binary because of size 
if not os.path.exists('connectivity.pickle'):
    blocks = read_plot3D('[Path to your Plot3D File]', binary = False) # Paraview has issues with ASCII 
    write_plot3D('[Path to your Plot3D File]',blocks, binary = True)   # Converted to binary for easy reading 

    # Block 1 is the blade O-Mesh k=0
    face_matches, outer_faces_formatted = connectivity(blocks)
    with open('connectivity.pickle','wb') as f:
        pickle.dump({"face_matches":face_matches, "outer_faces":outer_faces_formatted},f)

This part of the code shows how you can read the pickle file in part 1 from within paraview to display all your plots.

This code was tested with Paraview version 5.9 https://www.paraview.org/download/

Paraview executable can be called with the script option to have it automatically create the plots for you.

For windows this looks like this. You have the [working folder]>”[path to paraview executable]” –script=[your_python_script.py]

Note

Has been tested on Paraview 5.10. Not backwards compatible with older paraview versions

C:\Github\Plot3D_utilities\python\test\paraview_plot > “C:\Program Files\ParaView 5.10.0-Windows-Python3.9-msvc2017-AMD64\bin\paraview.exe” –script=compressor_plot_test.py

import sys
import os
import pickle
from typing import List
sys.path.insert(0,os.getcwd()) # This allows you to select files locally
from pv_library import Load, ExtractBlocks, ExtractSurface
from paraview.simple import *
import random


def CreateSubset(block_source,voi:List[int],name:str,opacity:float=1):
    """Creates a subset within paraview to display the mesh 

    Args:
        block_source ([type]): [description]
        voi (List[int]): This is the volume of interest. When (I,J,K)MIN=(I,J,K)MAX it is a surface/face. When two quantities is the same e.g. IMIN=IMAX and JMIN=JMAX you have an edge. If none of them are equal then you have a volume 
        name (str): Name you want paraview to display in the left 
        opacity (float, optional): [description]. Defaults to 1.

    Returns:
        (tuple): tuple containing:

            - **extractSubset1** (paraview source): source object representing the subset.
            - **extractSubset1Display** (paraview display): display settings for the source.
    """
    # Plot the Face 
    Hide(block_source, View)
    extractSubset1 = ExtractSubset(registrationName=name, Input=block_source)
    extractSubset1.VOI = voi

    renderView1 = GetActiveViewOrCreate('RenderView')
    extractSubset1Display = Show(extractSubset1, renderView1, 'StructuredGridRepresentation')
    # trace defaults for the display properties.
    extractSubset1Display.Representation = 'Surface'
    extractSubset1Display.ColorArrayName = [None, '']
    extractSubset1Display.SelectTCoordArray = 'None'
    extractSubset1Display.SelectNormalArray = 'None'
    extractSubset1Display.SelectTangentArray = 'None'
    extractSubset1Display.OSPRayScaleFunction = 'PiecewiseFunction'
    extractSubset1Display.SelectOrientationVectors = 'None'
    extractSubset1Display.ScaleFactor = 7.410221099853516
    extractSubset1Display.SelectScaleArray = 'None'
    extractSubset1Display.GlyphType = 'Arrow'
    extractSubset1Display.GlyphTableIndexArray = 'None'
    extractSubset1Display.GaussianRadius = 0.3705110549926758
    extractSubset1Display.SetScaleArray = [None, '']
    extractSubset1Display.ScaleTransferFunction = 'PiecewiseFunction'
    extractSubset1Display.OpacityArray = [None, '']
    extractSubset1Display.OpacityTransferFunction = 'PiecewiseFunction'
    extractSubset1Display.DataAxesGrid = 'GridAxesRepresentation'
    extractSubset1Display.PolarAxes = 'PolarAxesRepresentation'
    extractSubset1Display.ScalarOpacityUnitDistance = 6.758095838007181
    extractSubset1Display.SetRepresentationType('Surface')
    extractSubset1Display.Opacity = opacity
    # Add in the face color and update 
    extractSubset1Display.AmbientColor = rgb_face_matches[match_indx]
    extractSubset1Display.DiffuseColor = rgb_face_matches[match_indx]
    renderView1.Update()

    return extractSubset1, extractSubset1Display

'''
Main Code
'''

if __name__=="__main__":
    '''
    Read the connectivity file
    '''
    with open('connectivity.pickle','rb') as f:
        data = pickle.load(f)
        face_matches = data['face_matches']
        outer_faces = data['outer_faces']
        if 'periodic_surfaces' in data:
            periodic_faces = data['periodic_surfaces']
        else:
            periodic_faces = []

    blocks_to_extract = [f['block1']['block_index'] for f in face_matches]
    blocks_to_extract.extend([f['block2']['block_index'] for f in face_matches])
    blocks_to_extract = list(set(blocks_to_extract))

    '''
    Generate Random Colors 
    
    This part generates random colors so that each color is associated with a face match. 
    Doesn't matter what the block is a single match is assigned the same color. 
    '''
    rgb_face_matches = list()
    for i in range(len(face_matches)):
        rgb_face_matches.append([random.randint(0,255)/255, random.randint(0,255)/255, random.randint(0,255)/255])

    rgb_outer_faces = list()
    for i in range(len(outer_faces)):
        rgb_outer_faces.append([random.randint(0,255)/255, random.randint(0,255)/255, random.randint(0,255)/255])

    # Load mesh
    plot3d_binary_filename = 'R4_testfile_binary.xyz'
    plot3D_source,plot3D_Display,View,LUT = Load(plot3d_binary_filename)
    surface_indx = 1 
    
    '''
    Loop through all the blocks and create within each block the match and the outer surfaces
    '''
    for b in blocks_to_extract: # Block indicies 
        block_source,block_display,LUT = ExtractBlocks(plot3D_source,View,[b])
        RenameSource('Block '+str(b), block_source)
        block_source = FindSource('Block '+str(b))
        
        # Plot the face matches 
        for match_indx, f in enumerate(face_matches):
            # Add Plots for Matched Faces 
            if f['block1']['block_index'] == b or f['block2']['block_index'] == b : 
                if f['block1']['block_index'] == b:
                    voi = [f['block1']['lb'][0], f['block1']['ub'][0], f['block1']['lb'][1], f['block1']['ub'][1],f['block1']['lb'][2], f['block1']['ub'][2]]
                else:
                    voi = [f['block2']['lb'][0], f['block2']['ub'][0], f['block2']['lb'][1], f['block2']['ub'][1],f['block2']['lb'][2], f['block2']['ub'][2]]
                CreateSubset(block_source, voi, name='match '+str(match_indx))
        
        # Plot the outer faces  
        for o in outer_faces:
            # Add Plots for Outer Faces
            if o['block_index'] == b:
                    voi = [o['lb'][0], o['ub'][0], o['lb'][1], o['ub'][1],o['lb'][2], o['ub'][2]]
                    CreateSubset(block_source, voi, name='outer_face '+str(surface_indx),opacity=0.2)
                    surface_indx +=1 
               
        for periodic_indx, p in enumerate(periodic_faces):
            # Add Plots for Outer Faces
            if p['block1']['block_index'] == b and p['block2']['block_index'] == b: # Periodicity within the block 
                voi = [p['block1']['lb'][0], p['block1']['ub'][0], p['block1']['lb'][1], p['block1']['ub'][1],p['block1']['lb'][2], p['block1']['ub'][2]]
                CreateSubset(block_source, voi, name='periodic '+str(periodic_indx))
                voi = [p['block2']['lb'][0], p['block2']['ub'][0], p['block2']['lb'][1], p['block2']['ub'][1],p['block2']['lb'][2], p['block2']['ub'][2]]
                CreateSubset(block_source, voi, name='periodic '+str(periodic_indx))

            elif p['block1']['block_index'] == b or p['block2']['block_index'] == b: # Periodicity from block to block
                if p['block1']['block_index'] == b:
                    voi = [p['block1']['lb'][0], p['block1']['ub'][0], p['block1']['lb'][1], p['block1']['ub'][1],p['block1']['lb'][2], p['block1']['ub'][2]]
                else:
                    voi = [p['block2']['lb'][0], p['block2']['ub'][0], p['block2']['lb'][1], p['block2']['ub'][1],p['block2']['lb'][2], p['block2']['ub'][2]]
                CreateSubset(block_source, voi, name='periodic '+str(periodic_indx))

This is the contents of pv_library.py

from typing import List
from paraview.simple import *
import os
import math
#pylint: skip-file
paraview.simple._DisableFirstRenderCameraReset()

def Load(filename:str):
    """Calls pvpython and displays the file

    Args:
        filename (str): [description]

    Returns:
        (tuple): tuple containing:

            - **Plot3D** (paraview.source): source object for plot3d 
            - **plot3D_Display** (paraview display): this is the display object for the source. Colors and things 
            - **View** (paraview.source): source object for plot3d 

    """
    plot3D = PLOT3DReader(registrationName=filename,
            QFileName='',
            FileName=filename,
            FunctionFileName='')

    # set active source
    SetActiveSource(plot3D)
    # get active view
    View = GetActiveViewOrCreate('RenderView')    
    # show data in view
    plot3D_Display = Show(plot3D, View)
    # trace defaults for the display properties.
    plot3D_Display.Representation = 'Outline'
    plot3D_Display.ColorArrayName = ['POINTS', '']
    plot3D_Display.OSPRayScaleFunction = 'PiecewiseFunction'
    plot3D_Display.SelectOrientationVectors = 'None'
    plot3D_Display.ScaleFactor = 0.5376288175582886
    plot3D_Display.SelectScaleArray = 'None'
    plot3D_Display.GlyphType = 'Arrow'
    plot3D_Display.PolarAxes = 'PolarAxesRepresentation'
    plot3D_Display.ScalarOpacityUnitDistance = 0.11673090489063437

    # reset view to fit data
    View.ResetCamera()
    
    # show color bar/color legend
    plot3D_Display.SetRepresentationType('Outline')
    ColorBy(plot3D_Display, ('FIELD', 'vtkBlockColors'))
    LUT = GetColorTransferFunction('SolidColor')
    HideScalarBarIfNotNeeded(LUT, View) # Change to solid color
    return plot3D,plot3D_Display,View,LUT

def SetCamera(View,position,focalPoint,ViewUp):  
    # current camera placement for renderView1
    View.CameraPosition = position
    View.CameraFocalPoint = focalPoint
    View.CameraViewUp = ViewUp
    View.CameraParallelScale = 3.242501630387953

# Extracts the mesh block
# Block indicies should be an array format
def ExtractBlocks(source,View,BlockIndicies):
    extractBlock1 = ExtractBlock(Input=source)
    extractBlock1.Selectors = [f"/Root/Block{i}" for i in BlockIndicies]
    extractBlock1Display = Show(extractBlock1, View)
    extractBlock1Display.Representation = 'Outline'
    extractBlock1Display.ColorArrayName = ['POINTS', '']
    extractBlock1Display.OSPRayScaleFunction = 'PiecewiseFunction'
    extractBlock1Display.SelectOrientationVectors = 'None'
    extractBlock1Display.ScaleFactor = 0.2489664673805237
    extractBlock1Display.SelectScaleArray = 'None'
    extractBlock1Display.GlyphType = 'Arrow'
    extractBlock1Display.PolarAxes = 'PolarAxesRepresentation'
    extractBlock1Display.ScalarOpacityUnitDistance = 0.07851226208722488
    Hide(source, View)
    SetActiveSource(extractBlock1)
    # set scalar coloring
    ColorBy(extractBlock1Display, ('FIELD', 'vtkBlockColors'))

    # show color bar/color legend
    extractBlock1Display.SetScalarBarVisibility(View, True)

    # show color bar/color legend
    ColorBy(extractBlock1Display, ('FIELD', 'vtkBlockColors'))
    LUT = GetColorTransferFunction('vtkBlockColors')
    extractBlock1Display.SetRepresentationType('Surface With Edges')
    ColorBy(extractBlock1Display, ('FIELD', 'Solid Color'))
    HideScalarBarIfNotNeeded(LUT, View) # Change to solid color
    return extractBlock1,extractBlock1Display,LUT

def ExtractSurface(source,name:str,VOI:List[int]):
    """[summary]

    Args:
        source ([type]): [description]
        name (string): name of surface 
        VOI (List[int]): Volume of interest [imin,imax,jmin,jmax,kmin,kmax]
    Returns:
        [type]: [description]
    """
    renderView1 = GetActiveViewOrCreate('RenderView')

    Hide(source, renderView1)
    extractSubset1 = ExtractSubset(registrationName=name, Input=source)
    extractSubset1.VOI = voi

    renderView1 = GetActiveViewOrCreate('RenderView')
    extractSubset1Display = Show(extractSubset1, renderView1, 'StructuredGridRepresentation')
    # trace defaults for the display properties.
    extractSubset1Display.Representation = 'Outline'
    extractSubset1Display.ColorArrayName = [None, '']
    extractSubset1Display.SelectTCoordArray = 'None'
    extractSubset1Display.SelectNormalArray = 'None'
    extractSubset1Display.SelectTangentArray = 'None'
    extractSubset1Display.OSPRayScaleFunction = 'PiecewiseFunction'
    extractSubset1Display.SelectOrientationVectors = 'None'
    extractSubset1Display.ScaleFactor = 7.410221099853516
    extractSubset1Display.SelectScaleArray = 'None'
    extractSubset1Display.GlyphType = 'Arrow'
    extractSubset1Display.GlyphTableIndexArray = 'None'
    extractSubset1Display.GaussianRadius = 0.3705110549926758
    extractSubset1Display.SetScaleArray = [None, '']
    extractSubset1Display.ScaleTransferFunction = 'PiecewiseFunction'
    extractSubset1Display.OpacityArray = [None, '']
    extractSubset1Display.OpacityTransferFunction = 'PiecewiseFunction'
    extractSubset1Display.DataAxesGrid = 'GridAxesRepresentation'
    extractSubset1Display.PolarAxes = 'PolarAxesRepresentation'
    extractSubset1Display.ScalarOpacityUnitDistance = 6.758095838007181

    # update the view to ensure updated data information
    renderView1.Update()

    extractSubset1Display.SetRepresentationType('Outline')

    return extractSubset1,extractSubset1Display

def ChangeRepresentationMesh(source,sourceDisplay,View,LUT):
    SetActiveSource(source) # set active source
    HideScalarBarIfNotNeeded(LUT, View) # Hide the scalar bar for this color map if no visible data is colored by it.    
    sourceDisplay.SetRepresentationType('Surface With Edges') # set active source
    ColorBy(sourceDisplay, ('FIELD', 'Solid Color'))


def Periodicity(source,nBlades,BlockIndices,View):  
    # hide data in view
    Hide(source, View)
    angularPeriodicFilter1 = AngularPeriodicFilter(Input=source)

    # Properties modified on angularPeriodicFilter1
    angularPeriodicFilter1.BlockIndices = BlockIndices
    angularPeriodicFilter1.IterationMode = 'Manual'
    angularPeriodicFilter1.NumberOfPeriods = 2
    angularPeriodicFilter1.RotationAngle = 360.0/nBlades

    SetActiveSource(angularPeriodicFilter1)
    # get color transfer function/color map for 'Density'
    LUT = GetColorTransferFunction('Density')

    # show data in view
    angularPeriodicFilter1Display = Show(angularPeriodicFilter1, View)
    # trace defaults for the display properties.
    angularPeriodicFilter1Display.Representation = 'Outline'
    angularPeriodicFilter1Display.ColorArrayName = ['POINTS', 'Density']
    angularPeriodicFilter1Display.LookupTable = LUT
    angularPeriodicFilter1Display.OSPRayScaleArray = 'Density'
    angularPeriodicFilter1Display.OSPRayScaleFunction = 'PiecewiseFunction'
    angularPeriodicFilter1Display.SelectOrientationVectors = 'Momentum'
    angularPeriodicFilter1Display.ScaleFactor = 0.6140377521514893
    angularPeriodicFilter1Display.SelectScaleArray = 'Density'
    angularPeriodicFilter1Display.GlyphType = 'Arrow'
    angularPeriodicFilter1Display.PolarAxes = 'PolarAxesRepresentation'
    angularPeriodicFilter1Display.ScalarOpacityUnitDistance = 0.35345957752629076
    # show color bar/color legend
    angularPeriodicFilter1Display.SetScalarBarVisibility(View, True)
    return angularPeriodicFilter1,angularPeriodicFilter1Display,LUT


def CreateNewLayout(source,layoutName):
    CreateLayout(layoutName)

    # Create a new 'Render View'
    renderView2 = CreateView('RenderView')
    #renderView2.ViewSize = [1233, 814]
    renderView2.AxesGrid = 'GridAxes3DActor'
    renderView2.StereoType = 0
    renderView2.Background = [0.32, 0.34, 0.43]
    layout = GetLayout()
    
    # set active source
    if (source):
        SetActiveSource(source)
        # show data in view
        sourceDisplay = Show(source, renderView2)
        # trace defaults for the display properties.
        sourceDisplay.Representation = 'Outline'
        sourceDisplay.ColorArrayName = ['POINTS', '']
        sourceDisplay.OSPRayScaleArray = 'Density'
        sourceDisplay.OSPRayScaleFunction = 'PiecewiseFunction'
        sourceDisplay.SelectOrientationVectors = 'Momentum'
        sourceDisplay.ScaleFactor = 0.5376288175582886
        sourceDisplay.SelectScaleArray = 'Density'
        sourceDisplay.GlyphType = 'Arrow'
        sourceDisplay.PolarAxes = 'PolarAxesRepresentation'
        sourceDisplay.ScalarOpacityUnitDistance = 0.11673090489063437
        # reset view to fit data
        renderView2.ResetCamera()
        return layout,sourceDisplay,renderView2
    
    return layout,renderView2

The output should look like the following