Repository: DARMA-tasking/LB-analysis-framework

Issue:
ParaView scripting for post-hoc visualization
The goal of this issue is to add in-situ 2D visualization to LBAF, in order to give the user the option to entirely skip ParaView-based post-processing.

Discussion:

Patch:
diff --git a/src/lbaf/Applications/AnimationViewer.py b/src/lbaf/Applications/AnimationViewer.py
deleted file mode 100644
index 58b562c..0000000
--- a/src/lbaf/Applications/AnimationViewer.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import sys
-
-import paraview.simple as pv
-
-from lbaf.Applications.ParaviewViewer import ParaviewViewer
-from lbaf.Applications.ParaviewViewerBase import ViewerParameters, ParaviewViewerBase
-from lbaf.Utils.logger import logger
-
-
-class AnimationViewer(ParaviewViewer):
-    """ A concrete class providing an Animation Viewer
-    """
-
-    def __init__(self, exodus=None, file_name=None, viewer_type=None):
-        # Call superclass init
-        super().__init__(exodus, file_name, viewer_type)
-
-        # Starting logger
-        self.__logger = logger()
-
-    def saveView(self, reader):
-        """ Save animation
-        """
-        # Get animation scene
-        animationScene = pv.GetAnimationScene()
-        animationScene.PlayMode = "Snap To TimeSteps"
-
-        # Save animation images
-        for t in reader.TimestepValues.GetData()[:]:
-            animationScene.AnimationTime = t
-
-        # Save animation movie
-        self.__logger.info("###  Generating AVI animation...")
-        pv.AssignViewToLayout()
-        pv.WriteAnimation(f"{self.file_name}.avi", Magnification=1, Quality=2, FrameRate=1.0, Compression=True)
-        self.__logger.info(f"### AVI animation generated.")
-
-
-if __name__ == '__main__':
-    # Assign logger to variable
-    lgr = logger()
-
-    # Print startup information
-    sv = sys.version_info
-    lgr.info(f"### Started with Python {sv.major}.{sv.minor}.{sv.micro}")
-
-    # Instantiate parameters and set values from command line arguments
-    lgr.info("Parsing command line arguments")
-
-    params = ViewerParameters()
-    if params.parse_command_line():
-        raise SystemExit(1)
-
-    # Check if arguments were correctly parsed
-    animationViewer = ParaviewViewerBase.factory(params.exodus, params.file_name, "Animation")
-
-    # Create view from AnimationViewer instance
-    reader = animationViewer.createViews()
-
-    # Save generated view
-    animationViewer.saveView(reader)
-
-    # If this point is reached everything went fine
-    lgr.info(f"{animationViewer.file_name} file views generated ###")
diff --git a/src/lbaf/Applications/LBAF_app.py b/src/lbaf/Applications/LBAF_app.py
index 1adfc8c..f947511 100644
--- a/src/lbaf/Applications/LBAF_app.py
+++ b/src/lbaf/Applications/LBAF_app.py
@@ -21,7 +21,7 @@ from lbaf.Applications.rank_object_enumerator import compute_min_max_arrangement
 from lbaf.Execution.lbsRuntime import Runtime
 from lbaf.IO.configurationValidator import ConfigurationValidator
 from lbaf.IO.lbsVTDataWriter import VTDataWriter
-from lbaf.IO.lbsMeshWriter import MeshWriter
+from lbaf.IO.lbsMeshBasedVisualizer import MeshBasedVisualizer
 import lbaf.IO.lbsStatistics as lbstats
 from lbaf.Model.lbsPhase import Phase
 from lbaf.Utils.logger import logger
@@ -204,16 +204,16 @@ class LBAFApp:
             for phase_id in self.params.phase_ids:
                 # Create a phase and populate it
                 if "file_suffix" in self.params.__dict__:
-                    phase = Phase(self.logger, 0, self.params.file_suffix)
+                    phase = Phase(self.logger, phase_id, self.params.file_suffix)
                 else:
-                    phase = Phase(self.logger, 0)
+                    phase = Phase(self.logger, phase_id)
                 phase.populate_from_log(
                     self.params.n_ranks,
                     phase_id,
                     self.params.data_stem)
                 phases.append(phase)
         else:
-            # Populate phase pseudo-randomly
+            # Populate phase pseudo-randomly a phase 0
             phase = Phase(self.logger, 0)
             phase.populate_from_samplers(
                 self.params.n_ranks,
@@ -292,33 +292,21 @@ class LBAFApp:
                 output_dir=self.params.output_dir)
             vt_writer.write()
 
-        # If prefix parsed from command line
-        if "generate_meshes" in self.params.__dict__:
-            # Instantiate phase to mesh writer if requested
-            ex_writer = MeshWriter(
-                phase_0.get_number_of_ranks(),
+        # Generate meshes and multimedia when requested
+        gen_meshes = self.params.__dict__.get("generate_meshes")
+        gen_mulmed = self.params.__dict__.get("generate_multimedia")
+        if gen_meshes or gen_mulmed:
+            # Instantiate mesh based visualizer and execute as requested
+            ex_writer = MeshBasedVisualizer(
+                self.logger,
+                phases,
                 self.params.grid_size,
                 self.params.object_jitter,
-                self.logger,
+                self.params.output_dir,
                 self.params.output_file_stem,
-                output_dir=self.params.output_dir,
-                )
-            ex_writer.write(phases, rt.distributions, rt.statistics)
-
-        # Create a viewer if paraview is available
-        file_name = self.params.output_file_stem
-        if self.params.__dict__.get("generate_multimedia") is not None \
-                and self.params.__dict__.get("generate_multimedia"):
-            from ParaviewViewerBase import ParaviewViewerBase
-            if self.params.output_dir is not None:
-                file_name = os.path.join(self.params.output_dir, file_name)
-                output_file_stem = file_name
-                viewer = ParaviewViewerBase.factory(
-                    exodus=output_file_stem,
-                    file_name=file_name,
-                    viewer_type='')
-                reader = viewer.createViews()
-            viewer.saveView(reader)
+                rt.distributions,
+                rt.statistics)
+            ex_writer.generate(gen_meshes, gen_mulmed)
 
         # Compute and print final rank load and edge volume statistics
         _, _, l_ave, _, _, _, _, _ = lbstats.print_function_statistics(
diff --git a/src/lbaf/Applications/MoveCountsViewer.py b/src/lbaf/Applications/MoveCountsViewer.py
index 0a930f2..df07981 100644
--- a/src/lbaf/Applications/MoveCountsViewer.py
+++ b/src/lbaf/Applications/MoveCountsViewer.py
@@ -13,10 +13,48 @@ import sys
 
 import vtk
 
-from lbaf.Applications.MoveCountsViewerParameters import MoveCountsViewerParameters
 from lbaf.Utils.logger import logger
 
 
+class MoveCountsViewerParameters:
+    """ A class to describe MoveCountsViewer parameters
+    """
+
+    def __init__(self, viewer):
+        # Set renderer parameters
+        self.renderer_background = [1, 1, 1]
+
+        # Set actor_vertices parameters
+        self.actor_vertices_screen_size = 50 if viewer.interactive else 5000
+        self.actor_vertices_color = [0, 0, 0]
+        self.actor_vertices_opacity = .3 if viewer.interactive else .5
+
+        # Set actor_labels parameters
+        self.actor_labels_color = [0, 0, 0]
+        self.actor_labels_font_size = 16 if viewer.interactive else 150
+        self.actor_edges_opacity = .5 if viewer.interactive else 1
+        self.actor_edges_line_width = 2 if viewer.interactive else 15
+
+        # Set actor_arrows parameters
+        self.actor_arrows_edge_glyph_position = .5
+        self.actor_arrows_source_scale = .075
+
+        # Set actor_bar parameters
+        self.actor_bar_number_of_labels = 2
+        self.actor_bar_width = .2
+        self.actor_bar_heigth = .08
+        self.actor_bar_position = [.4, .91]
+        self.actor_bar_title_color = [0, 0, 0]
+        self.actor_bar_label_color = [0, 0, 0]
+
+        # Set window parameters
+        self.window_size_x = 600
+        self.window_size_y = 600
+
+        # Set wti (WindowToImageFilter) parameters
+        self.wti_scale = 10
+
+
 class MoveCountsViewer:
     """ A class to describe MoveCountsViewer attributes
     """
@@ -48,12 +86,11 @@ class MoveCountsViewer:
         self.logger = logger()
         self.logging_level = "info"
 
-    def usage(self):
+    @staticmethod
+    def usage():
         """ Provide online help
         """
-
-        print("Usage:")
-
+        print("# Usage:")
         print("\t [-p <np>]   number of processors")
         print("\t [-f <fn>]   input file name")
         print("\t [-s]        input file format suffix")
@@ -66,7 +103,6 @@ class MoveCountsViewer:
     def parse_command_line(self):
         """ Parse command line
         """
-
         # Try to hash command line with respect to allowable flags
         try:
             opts, args = getopt.getopt(sys.argv[1:], "p:f:s:o:t:ih")
@@ -159,7 +195,6 @@ class MoveCountsViewer:
                     if src_id != i:
                         directed_moves[(src_id, i)] = directed_moves.get(
                             (src_id, i), 0) + 1
-
         # Compute undirected move counts
         undirected_moves = {
             (i, j): directed_moves.get((i, j), 0) + directed_moves.get(
@@ -308,6 +343,7 @@ class MoveCountsViewer:
         actor_edges.GetProperty().SetLineWidth(
             viewerParams.actor_edges_line_width)
         renderer.AddViewProp(actor_edges)
+
         # Reset camera to set it up based on edge actor
         renderer.ResetCamera()
 
@@ -359,17 +395,23 @@ class MoveCountsViewer:
             # Write PNG image
             writer = vtk.vtkPNGWriter()
             writer.SetInputConnection(wti.GetOutputPort())
-            writer.SetFileName(f"{self.output_file_name}.{self.output_file_suffix}")
+            writer.SetFileName(
+                f"{self.output_file_name}.{self.output_file_suffix}")
             writer.Write()
 
 
 if __name__ == "__main__":
+    # Default settings
     n_processors = 8
-    input_file_name = "data/data/lb_iter"
-    input_file_suffix = "out"
+    input_file_name = "../data/lb50-data/data"
+    input_file_suffix = "vom"
     output_file_name = "move_counts"
-    params = MoveCountsViewer(n_processors=n_processors, input_file_name=input_file_name,
-                              input_file_suffix=input_file_suffix, output_file_name=output_file_name, interactive=False)
+    params = MoveCountsViewer(
+        n_processors=n_processors,
+        input_file_name=input_file_name,
+        input_file_suffix=input_file_suffix,
+        output_file_name=output_file_name,
+        interactive=False)
 
     # Assign logger to variable
     lgr = params.logger
@@ -379,9 +421,9 @@ if __name__ == "__main__":
     lgr.info(f"### Started with Python {sv.major}.{sv.minor}.{sv.micro}")
 
     # Instantiate parameters and set values from command line arguments
-    lgr.info("Parsing command line arguments")
-
+    lgr.info("# Parsing command line arguments")
     if params.parse_command_line():
         raise SystemExit(1)
 
+    # Execute viewer
     params.computeMoveCountsViewer()
diff --git a/src/lbaf/Applications/MoveCountsViewerParameters.py b/src/lbaf/Applications/MoveCountsViewerParameters.py
index c08ad9f..e69de29 100644
--- a/src/lbaf/Applications/MoveCountsViewerParameters.py
+++ b/src/lbaf/Applications/MoveCountsViewerParameters.py
@@ -1,40 +0,0 @@
-class MoveCountsViewerParameters:
-    """ A class to describe MoveCountsViewer parameters
-    """
-
-    def __init__(self, viewer):
-
-        # Set parameters based on viewer's attribute values
-
-        # Set renderer parameters
-        self.renderer_background = [1, 1, 1]
-
-        # Set actor_vertices parameters
-        self.actor_vertices_screen_size = 50 if viewer.interactive else 5000
-        self.actor_vertices_color = [0, 0, 0]
-        self.actor_vertices_opacity = .3 if viewer.interactive else .5
-
-        # Set actor_labels parameters
-        self.actor_labels_color = [0, 0, 0]
-        self.actor_labels_font_size = 16 if viewer.interactive else 150
-        self.actor_edges_opacity = .5 if viewer.interactive else 1
-        self.actor_edges_line_width = 2 if viewer.interactive else 15
-
-        # Set actor_arrows parameters
-        self.actor_arrows_edge_glyph_position = .5
-        self.actor_arrows_source_scale = .075
-
-        # Set actor_bar parameters
-        self.actor_bar_number_of_labels = 2
-        self.actor_bar_width = .2
-        self.actor_bar_heigth = .08
-        self.actor_bar_position = [.4, .91]
-        self.actor_bar_title_color = [0, 0, 0]
-        self.actor_bar_label_color = [0, 0, 0]
-
-        # Set window parameters
-        self.window_size_x = 600
-        self.window_size_y = 600
-
-        # Set wti (WindowToImageFilter) parameters
-        self.wti_scale = 10
diff --git a/src/lbaf/Applications/PNGViewer.py b/src/lbaf/Applications/PNGViewer.py
deleted file mode 100644
index f65d7a0..0000000
--- a/src/lbaf/Applications/PNGViewer.py
+++ /dev/null
@@ -1,62 +0,0 @@
-import sys
-
-import paraview.simple as pv
-
-from lbaf.Applications.ParaviewViewer import ParaviewViewer
-from lbaf.Applications.ParaviewViewerBase import ViewerParameters, ParaviewViewerBase
-from lbaf.Utils.logger import logger
-
-
-class PNGViewer(ParaviewViewer):
-    """ A concrete class providing a PNG Viewer
-    """
-
-    def __init__(self, exodus=None, file_name=None, viewer_type=None):
-
-        # Call superclass init
-        super(PNGViewer, self).__init__(exodus, file_name, viewer_type)
-
-        # Starting logger
-        self.__logger = logger()
-
-    def saveView(self, reader):
-        """ Save figure
-        """
-
-        # Get animation scene
-        animationScene = pv.GetAnimationScene()
-        animationScene.PlayMode = "Snap To TimeSteps"
-
-        # Save animation images
-        self.__logger.info("###  Generating PNG images...")
-        for t in reader.TimestepValues.GetData()[:]:
-            animationScene.AnimationTime = t
-            pv.WriteImage(f"{self.file_name}.{t:.6f}.png")
-
-        self.__logger.info("### All PNG images generated.")
-
-
-if __name__ == '__main__':
-    # Assign logger to variable
-    lgr = logger()
-
-    # Print startup information
-    sv = sys.version_info
-    lgr.info(f"### Started with Python {sv.major}.{sv.minor}.{sv.micro}")
-
-    # Instantiate parameters and set values from command line arguments
-    lgr.info("Parsing command line arguments")
-
-    params = ViewerParameters()
-    if params.parse_command_line():
-        raise SystemExit(1)
-    pngViewer = ParaviewViewerBase.factory(params.exodus, params.file_name, "PNG")
-
-    # Create view from PNGViewer instance
-    reader = pngViewer.createViews()
-
-    # Save generated view
-    pngViewer.saveView(reader)
-
-    # If this point is reached everything went fine
-    lgr.info(f"{pngViewer.file_name} file views generated ###")
diff --git a/src/lbaf/Applications/ParaviewViewer.py b/src/lbaf/Applications/ParaviewViewer.py
deleted file mode 100644
index 4c895c0..0000000
--- a/src/lbaf/Applications/ParaviewViewer.py
+++ /dev/null
@@ -1,54 +0,0 @@
-import sys
-
-from lbaf.Applications.ParaviewViewerBase import ViewerParameters, ParaviewViewerBase
-from lbaf.Utils.logger import logger
-
-
-class ParaviewViewer(ParaviewViewerBase):
-    """ A concrete class providing a Paraview Viewer
-    """
-
-    def __init__(self, exodus=None, file_name=None, viewer_type=None):
-
-        # Call superclass init
-        super(ParaviewViewer, self).__init__(exodus, file_name, viewer_type)
-
-    def saveView(self, reader):
-        """ Save figure
-        """
-        from lbaf.Applications.AnimationViewer import AnimationViewer
-        from lbaf.Applications.PNGViewer import PNGViewer
-        self.__class__ = PNGViewer
-        self.saveView(reader)
-        self.__class__ = AnimationViewer
-        self.saveView(reader)
-
-
-if __name__ == '__main__':
-    # Assign logger to variable
-    lgr = logger()
-
-    # Print startup information
-    sv = sys.version_info
-    lgr.info(f"### Started with Python {sv.major}.{sv.minor}.{sv.micro}")
-
-    # Instantiate parameters and set values from command line arguments
-    lgr.info("Parsing command line arguments")
-    params = ViewerParameters()
-    if params.parse_command_line():
-        raise SystemExit(1)
-    viewer = ParaviewViewerBase.factory(params.exodus, params.file_name, "")
-
-    # Create view from PNGViewer instance
-    reader = viewer.createViews()
-
-    from lbaf.Applications.AnimationViewer import AnimationViewer
-    from lbaf.Applications.PNGViewer import PNGViewer
-    # Save generated view
-    viewer.__class__ = PNGViewer
-    viewer.saveView(reader)
-    viewer.__class__ = AnimationViewer
-    viewer.saveView(reader)
-
-    # If this point is reached everything went fine
-    lgr.info(f"{viewer.file_name} file views generated ###")
diff --git a/src/lbaf/Applications/ParaviewViewerBase.py b/src/lbaf/Applications/ParaviewViewerBase.py
deleted file mode 100644
index ad25c33..0000000
--- a/src/lbaf/Applications/ParaviewViewerBase.py
+++ /dev/null
@@ -1,416 +0,0 @@
-import abc
-import getopt
-import sys
-
-import paraview.simple as pv
-
-from lbaf.Utils.logger import logger
-
-# Assign logger to variable
-LGR = logger()
-
-
-class ViewerParameters:
-    """ A class to describe ParaviewViewerBase parameters
-    """
-
-    def usage(self):
-        """ Provide online help
-        """
-
-        print("Usage:")
-        print("\t [-e]        ExodusII file name")
-        print("\t [-f]        visualization file name")
-        print("\t [-h]        help: print this message and exit")
-        print('')
-
-    def parse_command_line(self):
-        """ Parse command line
-        """
-        # Try to hash command line with respect to allowable flags
-        try:
-            opts, args = getopt.getopt(sys.argv[1:], "he:f:")
-
-        except getopt.GetoptError:
-            LGR.error("Incorrect command line arguments.")
-            self.usage()
-            return True
-
-        # Parse arguments and assign corresponding member variable values
-        for o, a in opts:
-            if o == "-h":
-                self.usage()
-                sys.exit(0)
-            elif o == "-e":
-                self.exodus = a
-            elif o == "-f":
-                self.file_name = a
-
-        # Ensure that exactly one ExodusII file has been provided
-        if not self.exodus:
-            LGR.error("Provide an ExodusII file")
-            return True
-
-        # Set default visualization file name prefix
-        if not self.file_name:
-            self.file_name = self.exodus
-
-        # Set viewer type
-        self.viewer_type = None
-
-        # Set material library
-        self.material_library = pv.GetMaterialLibrary()
-
-        # No line parsing error occurred
-        return False
-
-
-class ParaviewViewerBase(object):
-    __metaclass__ = abc.ABCMeta
-
-    @abc.abstractmethod
-    def __init__(self, exodus=None, file_name=None, viewer_type=None):
-
-        # ExodusII file to be displayed
-        self.exodus = f"{exodus}.e"
-
-        # visualization file name
-        self.file_name = f"{file_name}.e"
-
-        # Viewer type
-        self.viewer_type = viewer_type
-
-        # Material library
-        self.material_library = pv.GetMaterialLibrary()
-
-    @staticmethod
-    def factory(exodus, file_name, viewer_type):
-        """ Produce the necessary concrete backend instance
-        """
-        from AnimationViewer import AnimationViewer
-        from PNGViewer import PNGViewer
-        from ParaviewViewer import ParaviewViewer
-
-        # Unspecified ExodusII file name
-        if not exodus:
-            LGR.error("An ExodusII file name needs to be provided. Exiting.")
-            raise SystemExit(1)
-
-        # Unspecified visualization file name
-        if (not file_name) or file_name == "''":
-            LGR.warning("Visualization file name has not been provided. Using ExodusII file name by default.")
-            file_name = exodus
-
-        # PNG viewer
-        if viewer_type == "PNG":
-            ret_object = PNGViewer(exodus, file_name, viewer_type)
-
-        # Animation viewer
-        elif viewer_type == "Animation":
-            ret_object = AnimationViewer(exodus, file_name, viewer_type)
-
-        # Paraview viewer
-        elif viewer_type == "":
-            ret_object = ParaviewViewer(exodus, file_name)
-
-        # Unspecified viewer type
-        elif viewer_type == None:
-            LGR.error("A viewer type needs to be provided. Exiting.")
-            raise SystemExit(1)
-
-        # Unsupported viewer type
-        else:
-            LGR.error(f"{viewer_type} type viewer unsupported. Exiting.")
-            raise SystemExit(1)
-
-        # Report not instantiated
-        if not ret_object:
-            LGR.error(f"{viewer_type} viewer not instantiated. Exiting.")
-            raise SystemExit(1)
-
-        # Return instantiated object
-        ret_object.exodus = "{}.e".format(exodus)
-        ret_object.file_name = "{}.e".format(file_name)
-        ret_object.viewer_type = viewer_type
-        ret_object.material_library = pv.GetMaterialLibrary()
-        LGR.info(f"Instantiated {viewer_type} viewer.")
-        return ret_object
-
-    def get_exodus(self):
-        """ Convenience method to get ExodusII file name
-        """
-        # Return value of ExodusII file name
-        return self.exodus
-
-    def get_file_name(self):
-        """ Convenience method to get visualization file name
-        """
-        # Return value of visualization file name
-        return self.file_name
-
-    def get_viewer_type(self):
-        """ Convenience method to get viewer type
-        """
-        # Return value of viewer type
-        return self.viewer_type
-
-    def createRenderView(self, view_size=[1024, 1024]):
-        """ Create a new 'Render View'
-        """
-
-        renderView = pv.CreateView('RenderView')
-        if view_size:
-            renderView.ViewSize = view_size
-        renderView.InteractionMode = '2D'
-        renderView.AxesGrid = 'GridAxes3DActor'
-        renderView.OrientationAxesVisibility = 0
-        renderView.CenterOfRotation = [1.5, 1.5, 0.0]
-        renderView.StereoType = 0
-        renderView.CameraPosition = [1.5, 1.5, 10000.0]
-        renderView.CameraFocalPoint = [1.5, 1.5, 0.0]
-        renderView.CameraParallelScale = 2.1213203435596424
-        renderView.CameraParallelProjection = 1
-        renderView.Background = [1.0, 1.0, 1.0]
-        renderView.OSPRayMaterialLibrary = self.material_library
-
-        # init the 'GridAxes3DActor' selected for 'AxesGrid'
-        renderView.AxesGrid.XTitleFontFile = ''
-        renderView.AxesGrid.YTitleFontFile = ''
-        renderView.AxesGrid.ZTitleFontFile = ''
-        renderView.AxesGrid.XLabelFontFile = ''
-        renderView.AxesGrid.YLabelFontFile = ''
-        renderView.AxesGrid.ZLabelFontFile = ''
-
-        return renderView
-
-    def createExodusIIReader(self, elt_var, pt_var):
-        """ Create a new 'ExodusIIReader'
-        """
-
-        reader = pv.ExodusIIReader(FileName=[self.exodus])
-        reader.GenerateObjectIdCellArray = 0
-        reader.GenerateGlobalElementIdArray = 0
-        reader.ElementVariables = [elt_var]
-        reader.GenerateGlobalNodeIdArray = 0
-        reader.PointVariables = [pt_var]
-        reader.GlobalVariables = []
-        reader.ElementBlocks = ['Unnamed block ID: 3 Type: edge']
-
-        return reader
-
-    def createCalculator(self, reader, fct, var):
-        """ Create a new 'Calculator'
-        """
-
-        calculator = pv.Calculator(Input=reader)
-        calculator.ResultArrayName = "{}_{}".format(fct, var.lower())
-        calculator.Function = "{}({})".format(fct, var)
-
-        return calculator
-
-    def createGlyph(self, input, type='Box', factor=0.1, mode="All Points"):
-        """ Create a new 'Glyph'
-        """
-
-        glyph = pv.Glyph(Input=input, GlyphType=type)
-        glyph.OrientationArray = ['POINTS', 'No orientation array']
-        glyph.ScaleArray = ['POINTS', '{}'.format(input.ResultArrayName)]
-        glyph.ScaleFactor = factor
-        glyph.GlyphTransform = 'Transform2'
-        glyph.GlyphMode = mode
-
-        return glyph
-
-    def createColorTransferFunction(self, var, colors=None, nan_color=[1., 1., 1.], nan_opacity=None,
-                                    auto_rescale_range_mode="Never"):
-        """ Create a color transfer function/color map
-        """
-        # get color transfer function/color map
-        fct = pv.GetColorTransferFunction(var)
-        if auto_rescale_range_mode:
-            fct.AutomaticRescaleRangeMode = auto_rescale_range_mode
-        if colors:
-            fct.RGBPoints = colors
-        if nan_color:
-            fct.NanColor = nan_color
-        if nan_opacity is not None:
-            fct.NanOpacity = nan_opacity
-        fct.ScalarRangeInitialized = 1.0
-
-        return fct
-
-    def createOpacityTransferFunction(self, var, points=None):
-        """ Create an opacity transfer function/color map
-        """
-        # get color transfer function/color map
-        fct = pv.GetOpacityTransferFunction(var)
-        if points:
-            fct.Points = points
-        fct.ScalarRangeInitialized = 1
-
-        return fct
-
-    def createDisplay(self, reader, renderView, array_name, color_transfert_fct, line_width=None, scale_factor=0.3,
-                      glyph_type="Box", opacity_fct=None):
-        """ Create a 'Display'
-        """
-        # Show data from reader
-        display = pv.Show(reader, renderView)
-
-        display.Representation = 'Surface'
-        display.ColorArrayName = array_name
-        display.LookupTable = color_transfert_fct
-        if line_width:
-            display.LineWidth = line_width
-        display.OSPRayScaleArray = 'Load'
-        display.OSPRayScaleFunction = 'PiecewiseFunction'
-        display.SelectOrientationVectors = 'Load'
-        if scale_factor:
-            display.ScaleFactor = scale_factor
-        display.SelectScaleArray = 'Load'
-        if glyph_type:
-            display.GlyphType = glyph_type
-        display.GlyphTableIndexArray = 'Load'
-        display.GaussianRadius = 0.015
-        display.SetScaleArray = array_name
-        display.ScaleTransferFunction = 'PiecewiseFunction'
-        display.OpacityArray = array_name
-        display.OpacityTransferFunction = 'PiecewiseFunction'
-        display.DataAxesGrid = 'GridAxesRepresentation'
-        display.SelectionCellLabelFontFile = ''
-        display.SelectionPointLabelFontFile = ''
-        display.PolarAxes = 'PolarAxesRepresentation'
-        if opacity_fct:
-            display.ScalarOpacityFunction = opacity_fct
-            display.ScalarOpacityUnitDistance = 0.8601532551232605
-
-        # init the 'GridAxesRepresentation' selected for 'DataAxesGrid'
-        display.DataAxesGrid.XTitleFontFile = ''
-        display.DataAxesGrid.YTitleFontFile = ''
-        display.DataAxesGrid.ZTitleFontFile = ''
-        display.DataAxesGrid.XLabelFontFile = ''
-        display.DataAxesGrid.YLabelFontFile = ''
-        display.DataAxesGrid.ZLabelFontFile = ''
-
-        # init the 'PolarAxesRepresentation' selected for 'PolarAxes'
-        display.PolarAxes.PolarAxisTitleFontFile = ''
-        display.PolarAxes.PolarAxisLabelFontFile = ''
-        display.PolarAxes.LastRadialAxisTextFontFile = ''
-        display.PolarAxes.SecondaryRadialAxesTextFontFile = ''
-
-        return display
-
-    @abc.abstractmethod
-    def saveView(self, reader):
-        """ Save view
-        """
-        pass
-
-    def createViews(self):
-        """ Create views
-        """
-        # Disable automatic camera reset on 'Show'
-        pv._DisableFirstRenderCameraReset()
-
-        # Create render view
-        renderView = self.createRenderView([900, 900])
-
-        # Activate render view
-        pv.SetActiveView(renderView)
-
-        # Create ExodusII reader
-        reader = self.createExodusIIReader("Volume", "Load")
-
-        # Create sqrt(load) calculator to optimize visuals
-        sqrt_load = self.createCalculator(reader, "sqrt", "Load")
-
-        # Create sqrt(load) glyph
-        glyph = self.createGlyph(sqrt_load, factor=0.05)
-
-        # Instantiate volume colors and points
-        volume_colors = [223.48540319420192,
-                         0.231373,
-                         0.298039,
-                         0.752941,
-                         784.8585271892204,
-                         0.865003,
-                         0.865003,
-                         0.865003,
-                         1346.2316511842387,
-                         0.705882,
-                         0.0156863,
-                         0.14902]
-        volume_points = [223.48540319420192,
-                         0.0,
-                         0.5,
-                         0.0,
-                         1346.2316511842387,
-                         1.0,
-                         0.5,
-                         0.0]
-        # Create color transfert functions
-        volumeLUT = self.createColorTransferFunction(
-            "Volume",
-            volume_colors,
-            [1., 1., 1.],
-            0.0)
-        volumePWF = self.createOpacityTransferFunction(
-            "Volume",
-            volume_points)
-
-        readerDisplay = self.createDisplay(
-            reader,
-            renderView,
-            ['CELLS', 'Volume'],
-            volumeLUT,
-            4.0,
-            None,
-            None,
-            volumePWF)
-
-        # Instantiate load colors and points
-        load_colors = [0.0,
-                       0.231373,
-                       0.298039,
-                       0.752941,
-                       130.73569142337513,
-                       0.865003,
-                       0.865003,
-                       0.865003,
-                       261.47138284675026,
-                       0.705882,
-                       0.0156863,
-                       0.14902]
-        load_points = [0.0,
-                       0.0,
-                       0.5,
-                       0.0,
-                       261.47138284675026,
-                       1.0,
-                       0.5,
-                       0.0]
-
-        # Create color transfert functions
-        loadLUT = self.createColorTransferFunction(
-            "Load",
-            load_colors,
-            [1., 1., 1.],
-            None,
-            "Never")
-        loadPWF = self.createOpacityTransferFunction(
-            "Load",
-            load_points)
-
-        # Create displays
-        glyphDisplay = self.createDisplay(
-            glyph,
-            renderView,
-            ['POINTS', 'Load'],
-            loadLUT,
-            None,
-            0.005)
-
-        # Activate glyph source
-        pv.SetActiveSource(glyph)
-
-        return reader
diff --git a/src/lbaf/Applications/every20phases.yaml b/src/lbaf/Applications/every20phases.yaml
new file mode 100644
index 0000000..65120b3
--- /dev/null
+++ b/src/lbaf/Applications/every20phases.yaml
@@ -0,0 +1,92 @@
+# Docs
+# param_name_in_conf [type]    Description
+# work_model                   work model to be used
+#    name [str]                in LoadOnly, AffineCombination
+#    parameters: [dict]        optional parameters specific to each work model
+# algorithm                    balancing algorithm to be used
+#    name [str]                in InformAndTransfer, BruteForce
+#    parameters: [dict]        parameters specitic to each algorithm:
+#      InformAndtransfer:
+#        criterion [str]       in Tempered (default), StrictLocalizer
+#        n_iterations [int]    number of load-balancing iterations
+#        deterministic_transfer [bool] (default: False) for deterministic transfer
+#        n_rounds [int]        number of information rounds
+#        fanout [int]          information fanout index
+#        order_strategy [str]  ordering of objects for transfer
+#    	                       in arbitrary (default), element_id, increasing_times,
+#                              decreasing_times, increasing_connectivity,
+#                              fewest_migrations, small_objects
+#      BruteForce: 
+#        skip_transfer [bool]  (default: False) skip transfer phase
+#      PhaseStepper 
+# logging_level [str]          set to `info`, `debug`, `warning` or `error`
+# x_procs [int]                number of procs in x direction for rank visualization
+# y_procs [int]                number of procs in y direction for rank visualization
+# z_procs [int]                number of procs in z direction for rank visualization
+# data_stem [str]              base file name of VT load logs
+# phase_ids [list]             list of ids of phase to be read in VT load logs
+# map_file [str]               base file name for VT object/proc mapping
+# file_suffix [str]            file suffix of VT data files (default: "json")
+# output_dir [str]             output directory (default: '.')
+# terminal_background [str]    background color for terminal output
+# generate_meshes [bool]       generate mesh outputs (default: False)
+# generate_multimedia [bool]   generate multimedia visualization (default: False)
+# n_objects [int]              number of objects
+# n_mapped_ranks [int]         number of initially mapped processors
+# communication_degree [int]   object communication degree (no communication if 0)
+# time_sampler                 description of object times sampler: 
+#    name [str]                in uniform, lognormal
+#    parameters [list]         parameters e.g. 1.0,10.0 for lognormal
+# volume_sampler               description of object communication volumes sampler: 
+#    name [str]                in uniform, lognormal
+#    parameters [list]         parameters e.g. 1.0,10.0 for lognormal
+
+# Specify input
+from_data:
+  data_stem: "../../../data/nolb-8color-16nodes-every50phases/stats"
+  phase_ids:
+    - 2
+    - 52
+    - 102
+    - 152
+    - 202
+    - 252
+    - 302
+    - 352
+    - 402
+    - 452
+    - 502
+    - 552
+    - 602
+    - 652
+    - 702
+    - 752
+    - 802
+    - 852
+    - 902
+    - 952
+
+# Specify work model
+work_model:
+  name: AffineCombination
+  parameters:
+    alpha: 1.
+    beta: 1.0e-8
+    gamma: 0.
+
+# Specify algorithm
+algorithm:
+  name: PhaseStepper
+
+# Specify output
+#logging_level: debug
+terminal_background: light
+generate_multimedia: True
+output_dir: ../../../output
+output_file_stem: output_file
+n_ranks: 32
+generate_meshes:
+  x_ranks: 8
+  y_ranks: 4
+  z_ranks: 1
+  object_jitter: 0.5
diff --git a/src/lbaf/IO/lbsMeshBasedVisualizer.py b/src/lbaf/IO/lbsMeshBasedVisualizer.py
new file mode 100644
index 0000000..7e86cc3
--- /dev/null
+++ b/src/lbaf/IO/lbsMeshBasedVisualizer.py
@@ -0,0 +1,643 @@
+from logging import Logger
+import os
+import math
+import numbers
+import random
+import vtk
+
+from .lbsGridStreamer import GridStreamer
+from ..Model.lbsPhase import Phase
+
+
+class MeshBasedVisualizer:
+    """A class to visualize LBAF results via mesh files and VTK views."""
+
+    def __init__(self, logger: Logger, phases: list, grid_size: list, object_jitter=0.0, output_dir='.',
+                 output_file_stem="LBAF_out", distributions=None, statistics=None, resolution=1.):
+        """ Class constructor:
+            phases: list of Phase instances
+            grid_size: iterable containing grid sizes in each dimension
+            object_jitter: coefficient of random jitter with magnitude < 1
+            f: file name stem
+            r: grid_resolution value
+            output_dir: output directory
+        """
+        # Assign logger to instance variable
+        self.__logger = logger
+
+        if distributions is None:
+            distributions = {}
+        if statistics is None:
+            statistics = {}
+
+        # Make sure that Phase instances were passed
+        if not all([isinstance(p, Phase) for p in phases]):
+            self.__logger.error(
+                "Mesh writer expects a list of Phase instances as input")
+            raise SystemExit(1)
+        self.__phases = phases
+
+        # Ensure that all phases have the same number of ranks
+        n_r = phases[0].get_number_of_ranks()
+        if not all([p.get_number_of_ranks() == n_r for p in phases[1:]]):
+            self.__logger.error(
+                f"All phases must have {n_r} ranks as the first one")
+            raise SystemExit(1)
+        self.__n_ranks =  n_r
+        
+        # Ensure that specified grid resolution is correct
+        if not isinstance(resolution, numbers.Number) or resolution <= 0.:
+            self.__logger.error("Grid resolution must be a positive number")
+            raise SystemExit(1)
+        self.__grid_resolution = float(resolution)
+
+        # Determine available dimensions for object placement in ranks
+        self.__grid_size = grid_size
+        self.__rank_dims = [
+            d for d in range(3) if self.__grid_size[d] > 1]
+        self.__max_o_per_dim = 0
+
+        # Compute constant per object jitter
+        self.__jitter_dims = {
+            i: [(random.random() - 0.5) * object_jitter
+                if d in self.__rank_dims
+                else 0.0 for d in range(3)]
+            for i in self.__phases[0].get_object_ids()}
+
+        # Initialize maximum edge volume
+        self.__max_volume = 0.0
+
+        # Compute object time range
+        self.__time_range = [math.inf, 0.0]
+        for p in self.__phases:
+            for r in p.get_ranks():
+                for o in r.get_objects():
+                    # Update time range when necessary
+                    time = o.get_time()
+                    if time > self.__time_range[1]:
+                        self.__time_range[1] = time
+                    if time < self.__time_range[0]:
+                        self.__time_range[0] = time
+
+        # Assemble file and path names from constructor parameters
+        self.__rank_file_name = f"{output_file_stem}_rank_view.e"
+        self.__object_file_name = f"{output_file_stem}_object_view"
+        self.__output_dir = output_dir
+        if self.__output_dir is not None:
+            self.__rank_file_name = os.path.join(
+                self.__output_dir, self.__rank_file_name)
+            self.__object_file_name = os.path.join(
+                self.__output_dir, self.__object_file_name)
+            self.__visualization_file_name = os.path.join(
+                self.__output_dir, output_file_stem)
+
+        # Retrieve and verify rank attribute distributions
+        dis_l = distributions.get("load", [])
+        dis_w = distributions.get("work", [])
+        if not (n_dis := len(dis_l)) == len(dis_w):
+            self.__logger.error(
+                f"Both load and work distributions must have {n_dis} entries")
+            raise SystemExit(1)
+        self.__distributions = distributions
+        self.__work_range = (
+            min(min(dis_w, key=min)), max(max(dis_w, key=max)))
+
+        # Create attribute data arrays for rank loads and works
+        self.__loads, self.__works = [], []
+        for _ in range(n_dis):
+            # Create and append new load and work point arrays
+            l_arr, w_arr = vtk.vtkDoubleArray(), vtk.vtkDoubleArray()
+            l_arr.SetName("Load")
+            w_arr.SetName("Work")
+            l_arr.SetNumberOfTuples(self.__n_ranks)
+            w_arr.SetNumberOfTuples(self.__n_ranks)
+            self.__loads.append(l_arr)
+            self.__works.append(w_arr)
+
+        # Iterate over ranks and create rank mesh points
+        self.__rank_points = vtk.vtkPoints()
+        self.__rank_points.SetNumberOfPoints(self.__n_ranks)
+        for i in range(self.__n_ranks):
+            # Insert point based on Cartesian coordinates
+            self.__rank_points.SetPoint(i, [
+                self.__grid_resolution * c
+                for c in self.global_id_to_cartesian(
+                    i, self.__grid_size)])
+
+            # Set point attributes from distribution values
+            for l, (l_arr, w_arr) in enumerate(zip(self.__loads, self.__works)):
+                l_arr.SetTuple1(i, dis_l[l][i])
+                w_arr.SetTuple1(i, dis_w[l][i])
+
+        # Iterate over all possible rank links and create edges
+        self.__rank_lines = vtk.vtkCellArray()
+        index_to_edge = {}
+        edge_index = 0
+        for i in range(self.__n_ranks):
+            for j in range(i + 1, self.__n_ranks):
+                # Insert new link based on endpoint indices
+                line = vtk.vtkLine()
+                line.GetPointIds().SetId(0, i)
+                line.GetPointIds().SetId(1, j)
+                self.__rank_lines.InsertNextCell(line)
+
+                # Update flat index map
+                index_to_edge[edge_index] = frozenset([i, j])
+                edge_index += 1
+
+        # Number of edges is fixed due to vtkExodusIIWriter limitation
+        n_e = int(self.__n_ranks * (self.__n_ranks - 1) / 2)
+        self.__logger.info(f"Creating rank mesh with {self.__n_ranks} points and {n_e} edges")
+
+        # Create attribute data arrays for edge sent volumes
+        self.__volumes = []
+        for i, sent in enumerate(self.__distributions["sent"]):
+            # Reduce directed edges into undirected ones
+            u_edges = {}
+            for k, v in sent.items():
+                u_edges[frozenset(k)] = u_edges.setdefault(frozenset(k), 0.) + v
+
+            # Create and append new volume array for edges
+            v_arr = vtk.vtkDoubleArray()
+            v_arr.SetName("Largest Directed Volume")
+            v_arr.SetNumberOfTuples(n_e)
+            self.__volumes.append(v_arr)
+            
+            # Assign edge volume values
+            self.__logger.debug(f"\titeration {i} edges:")
+            for e, edge in index_to_edge.items():
+                v = u_edges.get(edge, float("nan"))
+                v_arr.SetTuple1(e, v)
+                if v > self.__max_volume:
+                    self.__max_volume = v
+                self.__logger.debug(f"\t{e} {edge}): {v}")
+
+        # Create and populate field arrays for statistics
+        self.__field_data = {}
+        for stat_name, stat_values in statistics.items():
+            # Skip non-list entries
+            if not isinstance(stat_values, list):
+                continue
+
+            # Create one singleton for each value of each statistic
+            for v in stat_values:
+                s_arr = vtk.vtkDoubleArray()
+                s_arr.SetNumberOfTuples(1)
+                s_arr.SetTuple1(0, v)
+                s_arr.SetName(stat_name)
+                self.__field_data.setdefault(stat_name, []).append(s_arr)
+
+    @staticmethod
+    def global_id_to_cartesian(flat_id, grid_sizes):
+        """ Map global index to its Cartesian grid coordinates."""
+        # Sanity check
+        n01 = grid_sizes[0] * grid_sizes[1]
+        if flat_id < 0 or flat_id >= n01 * grid_sizes[2]:
+            return None, None, None
+
+        # Compute successive Euclidean divisions
+        k, r = divmod(flat_id, n01)
+        j, i = divmod(r, grid_sizes[0])
+
+        # Return Cartesian coordinates
+        return i, j, k
+
+    def create_object_mesh(self, phase: Phase, object_mapping: set):
+        """ Map objects to polygonal mesh."""
+        # Retrieve number of mesh points and bail out early if empty set
+        n_o = phase.get_number_of_objects()
+        if not n_o:
+            self.__logger.warning("Empty list of objects, cannot write a mesh file")
+            return
+
+        # Compute number of communication edges
+        n_e = int(n_o * (n_o - 1) / 2)
+        self.__logger.info(
+            f"Creating object mesh with {n_o} points and {n_e} edges")
+
+        # Create point array for object times
+        t_arr = vtk.vtkDoubleArray()
+        t_arr.SetName("Time")
+        t_arr.SetNumberOfTuples(n_o)
+
+        # Create bit array for object migratability
+        b_arr = vtk.vtkBitArray()
+        b_arr.SetName("Migratable")
+        b_arr.SetNumberOfTuples(n_o)
+
+        # Create and size point set
+        points = vtk.vtkPoints()
+        points.SetNumberOfPoints(n_o)
+
+        # Iterate over ranks and objects to create mesh points
+        ranks = phase.get_ranks()
+        point_index, point_to_index, sent_volumes = 0, {}, []
+        for rank_id, objects in enumerate(object_mapping):
+            # Determine rank offsets
+            offsets = [
+                self.__grid_resolution * c
+                for c in self.global_id_to_cartesian(
+                    rank_id, self.__grid_size)]
+
+            # Compute local object block parameters
+            n_o_rank = len(objects)
+            n_o_per_dim = math.ceil(n_o_rank ** (
+                1. / len(self.__rank_dims)))
+            if n_o_per_dim > self.__max_o_per_dim:
+                self.__max_o_per_dim = n_o_per_dim
+            o_resolution = self.__grid_resolution / (n_o_per_dim + 1.)
+
+            # Iterate over objects and create point coordinates
+            self.__logger.debug(
+                f"Arranging a maximum of {n_o_per_dim} objects per dimension in {self.__rank_dims}")
+            rank_size = [n_o_per_dim
+                         if d in self.__rank_dims
+                         else 1 for d in range(3)]
+            centering = [0.5 * o_resolution * (n_o_per_dim - 1.)
+                         if d in self.__rank_dims
+                         else 0.0 for d in range(3)]
+
+            # Order objects of current rank
+            r = ranks[rank_id]
+            objects_list = sorted(objects, key=lambda x: x.get_id())
+            ordered_objects = {o: 0 for o in objects_list if r.is_sentinel(o)}
+            for o in objects_list:
+                if not r.is_sentinel(o):
+                    ordered_objects[o] = 1
+
+            # Add rank objects to points set
+            for i, (o, m) in enumerate(ordered_objects.items()):
+                # Insert point using offset and rank coordinates
+                points.SetPoint(point_index, [
+                    offsets[d] - centering[d] + (
+                        self.__jitter_dims[o.get_id()][d] + c) * o_resolution
+                    for d, c in enumerate(self.global_id_to_cartesian(
+                        i, rank_size))])
+                time = o.get_time()
+                t_arr.SetTuple1(point_index, time)
+                b_arr.SetTuple1(point_index, m)
+
+                # Update sent volumes
+                for k, v in o.get_sent().items():
+                    sent_volumes.append((point_index, k, v))
+
+                # Update maps and counters
+                point_to_index[o] = point_index
+                point_index += 1
+
+        # Summarize edges
+        edges = {
+            (tr[0], point_to_index[tr[1]]): tr[2]
+            for tr in sent_volumes}
+
+        # Iterate over all possible links and create edges
+        lines = vtk.vtkCellArray()
+        index_to_edge = {}
+        edge_index = 0
+        for i in range(n_o):
+            for j in range(i + 1, n_o):
+                # Insert new link based on endpoint indices
+                line = vtk.vtkLine()
+                line.GetPointIds().SetId(0, i)
+                line.GetPointIds().SetId(1, j)
+                lines.InsertNextCell(line)
+
+                # Update flat index map
+                index_to_edge[edge_index] = (i, j)
+                edge_index += 1
+
+        # Create and append volume array for edges
+        v_arr = vtk.vtkDoubleArray()
+        v_arr.SetName("Volume")
+        v_arr.SetNumberOfTuples(n_e)
+
+        # Assign edge volume values
+        self.__logger.debug(f"\tedges:")
+        for e in range(n_e):
+            v_arr.SetTuple1(e, edges.get(index_to_edge[e], float("nan")))
+            self.__logger.debug(f"\t{e} {index_to_edge[e]}): {v_arr.GetTuple1(e)}")
+
+        # Create and return VTK polygonal data mesh
+        pd_mesh = vtk.vtkPolyData()
+        pd_mesh.SetPoints(points)
+        pd_mesh.GetPointData().SetScalars(t_arr)
+        pd_mesh.GetPointData().AddArray(b_arr)
+        pd_mesh.SetLines(lines)
+        pd_mesh.GetCellData().SetScalars(v_arr)
+        return pd_mesh
+
+    @staticmethod
+    def create_color_transfer_function(attribute_range, scheme=None):
+        """ Create a color transfer function given attribute range."""
+
+        # Create color transfer function
+        ctf = vtk.vtkColorTransferFunction()
+        ctf.SetNanColorRGBA(1., 1., 1., 0.)
+
+        # Set color transfer function depending on chosen scheme
+        if scheme == "blue_to_red":
+            ctf.SetColorSpaceToDiverging()
+            mid_point = (attribute_range[0] + attribute_range[1]) * .5
+            ctf.AddRGBPoint(attribute_range[0], .231, .298, .753)
+            ctf.AddRGBPoint(mid_point, .865, .865, .865)
+            ctf.AddRGBPoint(attribute_range[1], .906, .016, .109)
+        elif scheme == "white_to_black":
+            ctf.AddRGBPoint(attribute_range[0], 1.0, 1.0, 1.0)
+            ctf.AddRGBPoint(attribute_range[1], 0.0, 0.0, 0.0)
+        else:
+            mid_point = (attribute_range[0] + attribute_range[1]) * .5
+            ctf.AddRGBPoint(attribute_range[0], .431, .761, .161)
+            ctf.AddRGBPoint(mid_point, .98, .992, .059)
+            ctf.AddRGBPoint(attribute_range[1], 1.0, .647, 0.0)
+
+        # Return color transfer function
+        return ctf
+
+    @staticmethod
+    def create_scalar_bar_actor(mapper, title, x, y):
+        """ Create scalar bar with default and custom parameters."""
+
+        # Instantiate scalar bar linked to given mapper
+        scalar_bar_actor = vtk.vtkScalarBarActor()
+        scalar_bar_actor.SetLookupTable(mapper.GetLookupTable())
+
+        # Set default parameters
+        scalar_bar_actor.SetOrientationToHorizontal()
+        scalar_bar_actor.SetNumberOfLabels(2)
+        scalar_bar_actor.SetHeight(0.08)
+        scalar_bar_actor.SetWidth(0.4)
+        scalar_bar_actor.SetLabelFormat("%.3E")
+        scalar_bar_actor.SetBarRatio(0.3)
+        scalar_bar_actor.DrawTickLabelsOn()
+        for text_prop in (
+            scalar_bar_actor.GetTitleTextProperty(),
+            scalar_bar_actor.GetLabelTextProperty()):
+            text_prop.SetColor(0.0, 0.0, 0.0)
+            text_prop.ItalicOff()
+            text_prop.BoldOff()
+            text_prop.SetFontFamilyToArial()
+
+        # Set custom parameters
+        scalar_bar_actor.SetTitle(title)
+        position = scalar_bar_actor.GetPositionCoordinate()
+        position.SetCoordinateSystemToNormalizedViewport()
+        position.SetValue(x, y, 0.0)
+
+        # Return created scalar bar actor
+        return scalar_bar_actor
+
+    def create_rendering_pipeline(self, iteration: int, pid: int, edge_width: int, glyph_factor: float, win_size: int,
+                                  object_mesh):
+        """ Create VTK-based pipeline all the way to render window."""
+        # Create rank mesh for current phase
+        rank_mesh = vtk.vtkPolyData()
+        rank_mesh.SetPoints(self.__rank_points)
+        rank_mesh.SetLines(self.__rank_lines)
+        rank_mesh.GetPointData().SetScalars(self.__works[iteration])
+
+        # Create renderer with parallel projection
+        renderer = vtk.vtkRenderer()
+        renderer.SetBackground(1.0, 1.0, 1.0)
+        renderer.GetActiveCamera().ParallelProjectionOn()
+
+        # Create square glyphs at ranks
+        rank_glyph = vtk.vtkGlyphSource2D()
+        rank_glyph.SetGlyphTypeToSquare()
+        rank_glyph.SetScale(.95)
+        rank_glyph.FilledOn()
+        rank_glyph.CrossOff()
+        rank_glypher = vtk.vtkGlyph2D()
+        rank_glypher.SetSourceConnection(rank_glyph.GetOutputPort())
+        rank_glypher.SetInputData(rank_mesh)
+        rank_glypher.SetScaleModeToDataScalingOff()
+
+        # Lower glyphs slightly for visibility
+        z_lower = vtk.vtkTransform()
+        z_lower.Translate(0.0, 0.0, -0.01)
+        trans = vtk.vtkTransformPolyDataFilter()
+        trans.SetTransform(z_lower)
+        trans.SetInputConnection(rank_glypher.GetOutputPort())
+
+        # Create mapper for rank glyphs
+        rank_mapper = vtk.vtkPolyDataMapper()
+        rank_mapper.SetInputConnection(trans.GetOutputPort())
+        rank_mapper.SetLookupTable(
+            self.create_color_transfer_function(self.__work_range))
+        rank_mapper.SetScalarRange(self.__work_range)
+
+        # Create rank work and its scalar bar actors
+        rank_actor = vtk.vtkActor()
+        rank_actor.SetMapper(rank_mapper)
+        work_actor = self.create_scalar_bar_actor(
+            rank_mapper, "Rank Work", 0.55, 0.9)
+        renderer.AddActor(rank_actor)
+        renderer.AddActor2D(work_actor)
+
+        # Create white to black look-up table
+        bw_lut = vtk.vtkLookupTable()
+        bw_lut.SetTableRange((0.0, self.__max_volume))
+        bw_lut.SetSaturationRange(0, 0)
+        bw_lut.SetHueRange(0, 0)
+        bw_lut.SetValueRange(1, 0)
+        bw_lut.SetNanColor(1.0, 1.0, 1.0, 0.0)
+        bw_lut.Build()
+
+        # Create mapper for inter-object edges
+        edge_mapper = vtk.vtkPolyDataMapper()
+        edge_mapper.SetInputData(object_mesh)
+        edge_mapper.SetScalarModeToUseCellData()
+        edge_mapper.SetScalarRange((0.0, self.__max_volume))
+        edge_mapper.SetLookupTable(bw_lut)
+
+        # Create communication volume and its scalar bar actors
+        edge_actor = vtk.vtkActor()
+        edge_actor.SetMapper(edge_mapper)
+        edge_actor.GetProperty().SetLineWidth(edge_width)
+        # edge_actor.GetProperty().SetOpacity(1.0)
+        volume_actor = self.create_scalar_bar_actor(
+            edge_mapper, "Inter-Object Volume", 0.05, 0.05)
+        renderer.AddActor(edge_actor)
+        renderer.AddActor2D(volume_actor)
+
+        # Compute square root of object times
+        sqrtT = vtk.vtkArrayCalculator()
+        sqrtT.SetInputData(object_mesh)
+        sqrtT.AddScalarArrayName("Time")
+        sqrtT_str = "sqrt(Time)"
+        sqrtT.SetFunction(sqrtT_str)
+        sqrtT.SetResultArrayName(sqrtT_str)
+        sqrtT.Update()
+        sqrtT_out = sqrtT.GetOutput()
+        sqrtT_out.GetPointData().SetActiveScalars("Migratable")
+
+        # Glyph sentinel and migratable objects separately
+        glyph_actors = []
+        for k, v in {0.0: "Square", 1.0: "Circle"}.items():
+            # Threshold by Migratable status
+            thresh = vtk.vtkThreshold()
+            thresh.SetInputData(sqrtT_out)
+            thresh.ThresholdBetween(k, k)
+            thresh.Update()
+            thresh_out = thresh.GetOutput()
+            if not thresh_out.GetNumberOfPoints():
+                continue
+            thresh_out.GetPointData().SetActiveScalars(
+                sqrtT_str)
+
+            # Glyph by square root of object times
+            glyph = vtk.vtkGlyphSource2D()
+            getattr(glyph, f"SetGlyphTypeTo{v}")()
+            glyph.SetResolution(32)
+            glyph.SetScale(1.0)
+            glyph.FilledOn()
+            glyph.CrossOff()
+            glypher = vtk.vtkGlyph3D()
+            glypher.SetSourceConnection(glyph.GetOutputPort())
+            glypher.SetInputData(thresh_out)
+            glypher.SetScaleModeToScaleByScalar()
+            glypher.SetScaleFactor(glyph_factor)
+            glypher.Update()
+            glypher.GetOutput().GetPointData().SetActiveScalars("Time")
+
+            # Raise glyphs slightly for visibility
+            z_raise = vtk.vtkTransform()
+            z_raise.Translate(0.0, 0.0, 0.01)
+            trans = vtk.vtkTransformPolyDataFilter()
+            trans.SetTransform(z_raise)
+            trans.SetInputData(glypher.GetOutput())
+
+            # Create mapper and actor for glyphs
+            glyph_mapper = vtk.vtkPolyDataMapper()
+            glyph_mapper.SetInputConnection(trans.GetOutputPort())
+            glyph_mapper.SetLookupTable(
+                self.create_color_transfer_function(
+                    self.__time_range, "blue_to_red"))
+            glyph_mapper.SetScalarRange(self.__time_range)
+            glyph_actor = vtk.vtkActor()
+            glyph_actor.SetMapper(glyph_mapper)
+            renderer.AddActor(glyph_actor)
+
+        # Create and add unique scalar bar for object time
+        time_actor = self.create_scalar_bar_actor(
+            glyph_mapper, "Object Time", 0.55, 0.05)
+        renderer.AddActor2D(time_actor)
+
+        # Create text actor to indicate iteration
+        text_actor = vtk.vtkTextActor()
+        text_actor.SetInput(f"Phase ID: {pid}\nIteration: {iteration}")
+        text_prop = text_actor.GetTextProperty()
+        text_prop.SetColor(0.0, 0.0, 0.0)
+        text_prop.ItalicOff()
+        text_prop.BoldOff()
+        text_prop.SetFontFamilyToArial()
+        text_prop.SetFontSize(72)
+        position = text_actor.GetPositionCoordinate()
+        position.SetCoordinateSystemToNormalizedViewport()
+        position.SetValue(0.1, 0.9, 0.0)
+        renderer.AddActor(text_actor)
+
+        # Create and return render window
+        renderer.ResetCamera()
+        render_window = vtk.vtkRenderWindow()
+        render_window.AddRenderer(renderer)
+        render_window.SetWindowName("LBAF")
+        render_window.SetSize(win_size, win_size)
+        return render_window
+
+    def generate(self, gen_meshes, gen_mulmed):
+        """ Generate mesh and multimedia outputs."""
+
+        # Write ExodusII rank mesh when requested
+        if gen_meshes:
+            # Create grid streamer
+            streamer = GridStreamer(
+                self.__rank_points,
+                self.__rank_lines,
+                self.__field_data,
+                [self.__loads, self.__works],
+                self.__volumes,
+                lgr=self.__logger)
+
+            # Write to ExodusII file when possible
+            if streamer.Error:
+                self.__logger.warning(
+                    f"Failed to instantiate a grid streamer for file {self.__rank_file_name}")
+            else:
+                self.__logger.info(
+                    f"Writing ExodusII file: {self.__rank_file_name}")
+                writer = vtk.vtkExodusIIWriter()
+                writer.SetFileName(self.__rank_file_name)
+                writer.SetInputConnection(streamer.Algorithm.GetOutputPort())
+                writer.WriteAllTimeStepsOn()
+                writer.Update()
+
+        # Determine whether phase must be updated
+        update_phase = True if len(
+            objects := self.__distributions.get("objects", set())
+            ) == len(self.__phases) else False
+
+        # Iterate over all object distributions
+        phase = self.__phases[0]
+        for iteration, object_mapping in enumerate(objects):
+            # Update phase when required
+            if update_phase:
+                phase = self.__phases[iteration]
+
+            # Create object mesh
+            object_mesh = self.create_object_mesh(phase, object_mapping)
+
+            # Write to VTP file when requested
+            if gen_meshes:
+                file_name = f"{self.__object_file_name}_{iteration:02d}.vtp"
+                self.__logger.info(f"Writing VTP file: {file_name}")
+                writer = vtk.vtkXMLPolyDataWriter()
+                writer.SetFileName(file_name)
+                writer.SetInputData(object_mesh)
+                writer.Update()
+
+            # Generate visualizations when requested
+            if gen_mulmed:
+                if len (self.__rank_dims) > 2:
+                    self.__logger.warning(
+                        "Visualization generation not yet implemented in 3-D")
+                    continue
+
+                # Compute visualization parameters
+                self.__logger.info(
+                    f"Generating 2-D visualization for iteration {iteration}:")
+                win_size = 800
+                self.__logger.info(
+                    f"\tnumber of pixels: {win_size}x{win_size}")
+                edge_width = 0.1 * win_size / max(self.__grid_size)
+                self.__logger.info(
+                    f"\tcommunication edges width: {edge_width:.2g}")
+                glyph_factor = self.__grid_resolution / (
+                    (self.__max_o_per_dim + 1)
+                    * math.sqrt(self.__time_range[1]))
+                self.__logger.info(
+                    f"\tobject glyphs scaling: {glyph_factor:.2g}")
+
+                # Run visualization pipeline
+                render_window = self.create_rendering_pipeline(
+                    iteration,
+                    phase.get_id(),
+                    edge_width,
+                    glyph_factor,
+                    win_size,
+                    object_mesh)
+                render_window.Render()
+
+                # Convert window to image
+                w2i = vtk.vtkWindowToImageFilter()
+                w2i.SetInput(render_window)
+                w2i.SetScale(3)
+                # w2i.SetInputBufferTypeToRGBA()
+
+                # Output PNG file
+                file_name = f"{self.__visualization_file_name}_{iteration:02d}.png"
+                self.__logger.info(f"Writing PNG file: {file_name}")
+                writer = vtk.vtkPNGWriter()
+                writer.SetInputConnection(w2i.GetOutputPort())
+                writer.SetFileName(file_name)
+                writer.SetCompressionLevel(2)
+                writer.Write()
diff --git a/src/lbaf/IO/lbsMeshWriter.py b/src/lbaf/IO/lbsMeshWriter.py
deleted file mode 100644
index f278cb8..0000000
--- a/src/lbaf/IO/lbsMeshWriter.py
+++ /dev/null
@@ -1,327 +0,0 @@
-from logging import Logger
-import os
-import sys
-import math
-import numbers
-import random
-import vtk
-
-from .lbsGridStreamer import GridStreamer
-from ..Model.lbsPhase import Phase
-
-
-class MeshWriter:
-    """A class to write LBAF results to mesh files via VTK layer."""
-
-    def __init__(
-        self,
-        n_r: int,
-        grid_size,
-        object_jitter: float,
-        logger: Logger,
-        f="lbs_out",
-        r=1.,
-        output_dir=None
-    ):
-        """ Class constructor:
-            n_r: number of ranks
-            grid_size: iterable containing grid sizes in each dimension
-            object_jitter: coefficient of random jitter with magnitude < 1
-            f: file name stem
-            r: grid_resolution value
-            output_dir: output directory
-        """
-        # Assign logger to instance variable
-        self.__logger = logger
-
-        # Ensure that specified grid resolution is correct
-        if not isinstance(r, numbers.Number) or r <= 0.:
-            self.__logger.error("Grid resolution must be a positive number")
-            raise SystemExit(1)
-        self.__grid_resolution = float(r)
-
-        # Keep track of mesh properties
-        self.__n_p =  n_r
-        self.__grid_size = grid_size
-        self.__object_jitter = object_jitter
-
-        # Assemble file and path names from constructor parameters
-        self.__rank_file_name = f"{f}_rank_view.e"
-        self.__object_file_name = f"{f}_object_view"
-        self.__output_dir = output_dir
-        if self.__output_dir is not None:
-            self.__rank_file_name = os.path.join(
-                self.__output_dir,
-                self.__rank_file_name)
-            self.__object_file_name = os.path.join(
-                self.__output_dir,
-                self.__object_file_name)
-
-    @staticmethod
-    def global_id_to_cartesian(flat_id, grid_sizes):
-        """ Map global index to its Cartesian grid coordinates."""
-        # Sanity check
-        n01 = grid_sizes[0] * grid_sizes[1]
-        if flat_id < 0 or flat_id >= n01 * grid_sizes[2]:
-            return None, None, None
-
-        # Compute successive Euclidean divisions
-        k, r = divmod(flat_id, n01)
-        j, i = divmod(r, grid_sizes[0])
-
-        # Return Cartesian coordinates
-        return i, j, k
-
-    def write_rank_view_file(self, ranks: list, distributions: dict, statistics: dict):
-        """ Map ranks to grid and write ExodusII file."""
-        # Number of edges is fixed due to vtkExodusIIWriter limitation
-        n_e = int(self.__n_p * (self.__n_p - 1) / 2)
-        self.__logger.info(f"Creating rank view mesh with {self.__n_p} points and {n_e} edges")
-
-        # Create and populate global field arrays for statistics
-        global_stats = {}
-        for stat_name, stat_values in statistics.items():
-            # Skip non-list entries
-            if not isinstance(stat_values, list):
-                continue
-
-            # Create one singleton for each value of each statistic
-            for v in stat_values:
-                s_arr = vtk.vtkDoubleArray()
-                s_arr.SetNumberOfTuples(1)
-                s_arr.SetTuple1(0, v)
-                s_arr.SetName(stat_name)
-                global_stats.setdefault(stat_name, []).append(s_arr)
-
-        # Create attribute data arrays for rank loads and works
-        loads, works = [], []
-        for _, _ in zip(distributions["load"], distributions["work"]):
-            # Create and append new load and work point arrays
-            l_arr, w_arr = vtk.vtkDoubleArray(), vtk.vtkDoubleArray()
-            l_arr.SetName("Load")
-            w_arr.SetName("Work")
-            l_arr.SetNumberOfTuples(self.__n_p)
-            w_arr.SetNumberOfTuples(self.__n_p)
-            loads.append(l_arr)
-            works.append(w_arr)
-
-        # Iterate over ranks and create mesh points
-        points = vtk.vtkPoints()
-        points.SetNumberOfPoints(self.__n_p)
-        for i, r in enumerate(ranks):
-            # Insert point based on Cartesian coordinates
-            points.SetPoint(i, [
-                self.__grid_resolution * c
-                for c in self.global_id_to_cartesian(
-                    r.get_id(), self.__grid_size)])
-            for l, (l_arr, w_arr) in enumerate(zip(loads, works)):
-                l_arr.SetTuple1(i, distributions["load"][l][i])
-                w_arr.SetTuple1(i, distributions["work"][l][i])
-
-        # Iterate over all possible links and create edges
-        lines = vtk.vtkCellArray()
-        index_to_edge = {}
-        edge_index = 0
-        for i in range(self.__n_p):
-            for j in range(i + 1, self.__n_p):
-                # Insert new link based on endpoint indices
-                line = vtk.vtkLine()
-                line.GetPointIds().SetId(0, i)
-                line.GetPointIds().SetId(1, j)
-                lines.InsertNextCell(line)
-
-                # Update flat index map
-                index_to_edge[edge_index] = frozenset([i, j])
-                edge_index += 1
-
-        # Create attribute data arrays for edge sent volumes
-        volumes = []
-        for i, sent in enumerate(distributions["sent"]):
-            # Reduce directed edges into undirected ones
-            u_edges = {}
-            for k, v in sent.items():
-                u_edges[frozenset(k)] = u_edges.setdefault(frozenset(k), 0.) + v
-
-            # Create and append new volume array for edges
-            v_arr = vtk.vtkDoubleArray()
-            v_arr.SetName("Largest Directed Volume")
-            v_arr.SetNumberOfTuples(n_e)
-            volumes.append(v_arr)
-            
-            # Assign edge volume values
-            self.__logger.debug(f"\titeration {i} edges:")
-            for e in range(n_e):
-                v_arr.SetTuple1(e, u_edges.get(index_to_edge[e], float("nan")))
-                self.__logger.debug(f"\t {e} {index_to_edge[e]}): {v_arr.GetTuple1(e)}")
-
-        # Create grid streamer
-        streamer = GridStreamer(points, lines, global_stats, [loads, works], volumes, lgr=self.__logger)
-
-        # Write to ExodusII file when possible
-        if streamer.Error:
-            self.__logger.error(f"Failed to instantiate a grid streamer for file {self.__rank_file_name}")
-            raise SystemExit(1)
-        else:
-            self.__logger.info(f"Writing ExodusII file: {self.__rank_file_name}")
-            writer = vtk.vtkExodusIIWriter()
-            writer.SetFileName(self.__rank_file_name)
-            writer.SetInputConnection(streamer.Algorithm.GetOutputPort())
-            writer.WriteAllTimeStepsOn()
-            writer.Update()
-
-    def write_object_view_file(self, phases: list, distributions: dict):
-        """ Map objects to grid and write ExodusII file."""
-        # Determine available dimensions for object placement in ranks
-        rank_dims = [d for d in range(3) if self.__grid_size[d] > 1]
-
-        # Compute constant per object jitter
-        jitter_dims = {
-            i: [(random.random() - 0.5) * self.__object_jitter
-                if d in rank_dims else 0.0 for d in range(3)]
-            for i in phases[0].get_object_ids()}
-
-        # Determine whether phase must be updated
-        update_phase = True if len(distributions["objects"]
-            ) == len(phases) else False
-        
-        # Iterate over all object distributions
-        phase = phases[0]
-        for iteration, object_mapping in enumerate(distributions["objects"]):
-            # Update phase when required
-            if update_phase:
-                phase = phases[iteration]
-
-            # Retrieve number of mesh points and bail out early if empty set
-            n_o = phase.get_number_of_objects()
-            if not n_o:
-                self.__logger.warning("Empty list of objects, cannot write a mesh file")
-                return
-
-            # Compute number of communication edges
-            n_e = int(n_o * (n_o - 1) / 2)
-            self.__logger.info(
-                f"Creating object view mesh with {n_o} points, " +
-                f"{n_e} edges, and jitter factor: {self.__object_jitter}")
-            
-            # Create point array for object times
-            t_arr = vtk.vtkDoubleArray()
-            t_arr.SetName("Time")
-            t_arr.SetNumberOfTuples(n_o)
-
-            # Create bit array for object migratability
-            b_arr = vtk.vtkBitArray()
-            b_arr.SetName("Migratable")
-            b_arr.SetNumberOfTuples(n_o)
-
-            # Create and size point set
-            points = vtk.vtkPoints()
-            points.SetNumberOfPoints(n_o)
-            
-            # Iterate over ranks and objects to create mesh points
-            ranks = phase.get_ranks()
-            point_index, point_to_index, sent_volumes = 0, {}, []
-            for rank_id, objects in enumerate(object_mapping):
-                # Determine rank offsets
-                offsets = [
-                    self.__grid_resolution * c
-                    for c in self.global_id_to_cartesian(rank_id, self.__grid_size)]
-
-                # Iterate over objects and create point coordinates
-                n_o_rank = len(objects)
-                n_o_per_dim = math.ceil(n_o_rank ** (1. / len(rank_dims)))
-                self.__logger.debug(f"Arranging a maximum of {n_o_per_dim} objects per dimension in {rank_dims}")
-                o_resolution = self.__grid_resolution / (n_o_per_dim + 1.)
-                rank_size = [n_o_per_dim if d in rank_dims else 1 for d in range(3)]
-                centering = [0.5 * o_resolution * (n_o_per_dim - 1.)
-                             if d in rank_dims else 0.0 for d in range(3)]
-
-                # Order objects of current rank
-                r = ranks[rank_id]
-                objects_list = sorted(objects, key=lambda x: x.get_id())
-                ordered_objects = {o: 0 for o in objects_list if r.is_sentinel(o)}
-                for o in objects_list:
-                    if not r.is_sentinel(o):
-                        ordered_objects[o] = 1
-
-                # Add rank objects to points set
-                for i, (o, m) in enumerate(ordered_objects.items()):
-                    # Insert point using offset and rank coordinates
-                    points.SetPoint(point_index, [
-                        offsets[d] - centering[d] + (
-                            jitter_dims[o.get_id()][d] + c) * o_resolution
-                        for d, c in enumerate(self.global_id_to_cartesian(
-                            i, rank_size))])
-                    t_arr.SetTuple1(point_index, o.get_time())
-                    b_arr.SetTuple1(point_index, m)
-
-                    # Update sent volumes
-                    for k, v in o.get_sent().items():
-                        sent_volumes.append((point_index, k, v))
-
-                    # Update maps and counters
-                    point_to_index[o] = point_index
-                    point_index += 1
-
-            # Summarize edges
-            edges = {
-                (tr[0], point_to_index[tr[1]]): tr[2]
-                for tr in sent_volumes}
-
-            # Iterate over all possible links and create edges
-            lines = vtk.vtkCellArray()
-            index_to_edge = {}
-            edge_index = 0
-            for i in range(n_o):
-                for j in range(i + 1, n_o):
-                    # Insert new link based on endpoint indices
-                    line = vtk.vtkLine()
-                    line.GetPointIds().SetId(0, i)
-                    line.GetPointIds().SetId(1, j)
-                    lines.InsertNextCell(line)
-
-                    # Update flat index map
-                    index_to_edge[edge_index] = (i, j)
-                    edge_index += 1
-
-            # Create and append volume array for edges
-            v_arr = vtk.vtkDoubleArray()
-            v_arr.SetName("Volume")
-            v_arr.SetNumberOfTuples(n_e)
-
-            # Assign edge volume values
-            self.__logger.debug(f"\titeration {iteration} edges:")
-            for e in range(n_e):
-                v_arr.SetTuple1(e, edges.get(index_to_edge[e], float("nan")))
-                self.__logger.debug(f"\t {e} {index_to_edge[e]}): {v_arr.GetTuple1(e)}")
-
-            # Create VTK polygonal data mesh
-            pd_mesh = vtk.vtkPolyData()
-            pd_mesh.SetPoints(points)
-            pd_mesh.GetPointData().SetScalars(t_arr)
-            pd_mesh.GetPointData().AddArray(b_arr)
-            pd_mesh.SetLines(lines)
-            pd_mesh.GetCellData().SetScalars(v_arr)
-
-            # Write to VTP file
-            file_name = f"{self.__object_file_name}_{iteration:02d}.vtp"
-            self.__logger.info(f"Writing VTP file: {file_name}")
-            writer = vtk.vtkXMLPolyDataWriter()
-            writer.SetFileName(file_name)
-            writer.SetInputData(pd_mesh)
-            writer.Update()
-
-    def write(self, phases: list, distributions: dict, statistics: dict):
-        """ Write rank and object ExodusII files."""
-
-        # Make sure that Phase instances were passed
-        if not all([isinstance(p, Phase) for p in phases]):
-            self.__logger.error(
-                "Mesh writer expects a list of Phase instances as input")
-            raise SystemExit(1)
-
-        # Write rank view file with global per-rank statistics
-        self.write_rank_view_file(
-            phases[0].get_ranks(), distributions, statistics)
-
-        # Write object view file
-        self.write_object_view_file(phases, distributions)
diff --git a/src/lbaf/Model/lbsPhase.py b/src/lbaf/Model/lbsPhase.py
index 6e21357..e4b258f 100644
--- a/src/lbaf/Model/lbsPhase.py
+++ b/src/lbaf/Model/lbsPhase.py
@@ -1,6 +1,5 @@
 from logging import Logger
 import random as rnd
-import sys
 
 from .lbsObject import Object
 from .lbsRank import Rank
@@ -14,15 +13,15 @@ class Phase:
     """ A class representing the state of collection of ranks with objects at a given round
     """
 
-    def __init__(self, logger: Logger, t: int = 0, file_suffix="json"):
+    def __init__(self, logger: Logger, pid: int = 0, file_suffix="json"):
         # Initialize empty list of ranks
         self.__ranks = []
 
         # Initialize null number of objects
         self.__n_objects = 0
 
-        # Default time-step/phase of this phase
-        self.__phase_id = t
+        # Index of this phase
+        self.__phase_id = pid
 
         # Assign logger to instance variable
         self.__logger = logger
@@ -34,6 +33,10 @@ class Phase:
         # Data files suffix(reading from data)
         self.__file_suffix = file_suffix
 
+    def get_id(self):
+        """ Retrieve index of this phase."""
+        return self.__phase_id
+
     def get_number_of_ranks(self):
         """ Retrieve number of ranks belonging to phase."""
         return len(self.__ranks)
@@ -46,10 +49,6 @@ class Phase:
         """ Retrieve IDs of ranks belonging to phase."""
         return [p.get_id() for p in self.__ranks]
 
-    def get_phase_id(self):
-        """ Retrieve the time-step/phase for this phase."""
-        return self.__phase_id
-
     def get_number_of_objects(self):
         """ Return number of objects."""
         return self.__n_objects


Test patch:
diff --git a/tests/test_lbs_phase.py b/tests/test_lbs_phase.py
index e9472fc..27af752 100644
--- a/tests/test_lbs_phase.py
+++ b/tests/test_lbs_phase.py
@@ -41,7 +41,7 @@ class TestConfig(unittest.TestCase):
         ranks = sorted([rank.get_id() for rank in self.phase.get_ranks()])
         self.assertEqual(ranks, [0, 1, 2, 3])
         self.assertEqual(sorted(self.phase.get_rank_ids()), [0, 1, 2, 3])
-        self.assertEqual(self.phase.get_phase_id(), 0)
+        self.assertEqual(self.phase.get_id(), 0)
 
     def test_lbs_phase_edges(self):
         file_prefix = os.path.join(self.data_dir, 'synthetic_lb_stats_compressed', 'data')