Split Blocks

When solving a plot3D block it is often useful to break it into smaller blocks of a certain size. This will improve the speed by splitting the blocks and allowing each CPU core to solve a part of the mesh. BUT we also need to maintain something called multi-grid.

Multi-grid concept

Mulit-grid is a concept where you take a gird say 4x4 and you solve it as a 2x2 then interpolate the results on to the larger grid. The idea of solving a coarse grid and translating the solution onto a finer grid allows you to reach a converged solution much faster. So that’s the benefits, what are the requirements?

Splitting Blocks Example

In this example we will use the file PahtCascade-ASCII We will go from reading -> splitting -> connectivity -> periodicity -> viewing it in paraview

Reading the file and split the blocks

When you specify the number of cells per block to split into, the code will try to match multigrid so the number of cells may change slightly

from plot3d import split_blocks, Direction, read_plot3D,write_plot3D

blocks = read_plot3D('PahtCascade-ASCII.xyz',binary=False)  # Reading plot3D
blocks_split = split_blocks(blocks,300000, direction=Direction.i)   # Here we split into chunks of size 300000 cells per block
write_plot3D('PahtCascade-Split.xyz',blocks_split,binary=True)

Finding Connectivity

face_matches, outer_faces_formatted = connectivity(blocks_split)
with open('connectivity-block-split.pickle','wb') as f:
    pickle.dump({"face_matches":face_matches, "outer_faces":outer_faces_formatted},f)

Finding Periodicity

This will find the periodic faces of the split block

with open('connectivity-block-split.pickle','rb') as f:
    data = pickle.load(f)
    face_matches = data['face_matches']
    outer_faces = data['outer_faces']

blocks = read_plot3D('PahtCascade-Split.xyz', binary = True, big_endian=True)
periodic_surfaces, outer_faces_to_keep,periodic_faces,outer_faces = periodicity(blocks,outer_faces,face_matches,periodic_direction='k',rotation_axis='x',nblades=55)
with open('connectivity-block-split_v02.pickle','wb') as f:
    [m.pop('match',None) for m in face_matches] # Remove the dataframe
    pickle.dump({"face_matches":face_matches, "outer_faces":outer_faces_to_keep, "periodic_surfaces":periodic_surfaces},f)

# Append periodic surfaces to face_matches
face_matches.extend(periodic_surfaces)

Plotting Split Blocks with Paraview

Library file that you will need in the same directory

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 

    """
    extension = os.path.basename(filename).split('.')[1]
    plot3D = PLOT3DReader(registrationName=os.path.splitext(filename)[0],
            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.BlockIndices = 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 = '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

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

    extractSubset1Display.SetRepresentationType('Surface')

    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 = 'Surface'
    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 script that will plot the mesh and periodicity + connectivity + outer faces

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

# Extracts the mesh block
# Block indicies should be an array format
def ExtractBlocks(source,View,BlockIndicies:List[int]):
    """Extracts block(s) from a mesh in paraview 

    Args:
        source ([unknown]): This is a paraview object that describes a source 
        View ([unknown]): Thisi s a paraview object that describes a view 
        BlockIndicies (List[int]): List of blocks you wish to extract to look at. Index should start at 1

    Returns:
        (tuple): tuple containing:

            - **extractBlock1** (paraview source): source object representing the block.
            - **extractBlock1Display** (paraview display): display settings for the source.
            - **extractBlock1Display** (paraview display): display settings for the source.
    """
    extractBlock1 = ExtractBlock(Input=source)
    extractBlock1.BlockIndices = 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 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-block-split_v02.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))
    blocks_to_extract.sort()
    print(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 = 'PahtCascade-Split.xyz'
    plot3D_source,plot3D_Display,View,LUT = Load(plot3d_binary_filename)
    
    '''
    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+1])
        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']['IMIN'], f['block1']['IMAX'], f['block1']['JMIN'], f['block1']['JMAX'],f['block1']['KMIN'], f['block1']['KMAX']]
                else:
                    voi = [f['block2']['IMIN'], f['block2']['IMAX'], f['block2']['JMIN'], f['block2']['JMAX'],f['block2']['KMIN'], f['block2']['KMAX']]
                CreateSubset(block_source, voi, name='match '+str(match_indx))
        
        # Plot the outer faces  
        for surface_indx, o in enumerate(outer_faces):
            # Add Plots for Outer Faces
            if o['block_index'] == b:
                voi = [o['IMIN'], o['IMAX'], o['JMIN'], o['JMAX'],o['KMIN'], o['KMAX']]
                CreateSubset(block_source, voi, name='outer_face '+str(surface_indx+1),opacity=0.2) 
        
        # Plot the outer faces  
        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']['IMIN'], p['block1']['IMAX'], p['block1']['JMIN'], p['block1']['JMAX'],p['block1']['KMIN'], p['block1']['KMAX']]
                CreateSubset(block_source, voi, name='periodic '+str(periodic_indx))
                voi = [p['block2']['IMIN'], p['block2']['IMAX'], p['block2']['JMIN'], p['block2']['JMAX'],p['block2']['KMIN'], p['block2']['KMAX']]
                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']['IMIN'], p['block1']['IMAX'], p['block1']['JMIN'], p['block1']['JMAX'],p['block1']['KMIN'], p['block1']['KMAX']]
                else:
                    voi = [p['block2']['IMIN'], p['block2']['IMAX'], p['block2']['JMIN'], p['block2']['JMAX'],p['block2']['KMIN'], p['block2']['KMAX']]
                CreateSubset(block_source, voi, name='periodic '+str(periodic_indx))