Python Filter Scripts

Python Filter Scripts

Expressions

Domain Id Expression

#
# [VisIt Python Expression]
#
# Labels zones with the id of domain they live in.
#
# CLI Usage (replace <mesh_name> with the name of your mesh) :
# >>DefinePythonExpression("domain_id",file="visit_domain_id_expr.vpe",args=["<mesh_name>"])
#
#

class DomainIdExpr(SimplePythonExpression):
    def __init__(self):
        SimplePythonExpression.__init__(self)
        self.name = "DomainIdExpr"
        self.description = "Labeling zones with domin id ..."
        self.output_is_point_var = False
        self.output_dimension = 1
    def derive_variable(self,ds_in,domain_id):
        ncells = ds_in.GetNumberOfCells()
        res = vtk.vtkFloatArray()
        res.SetNumberOfComponents(1)
        res.SetNumberOfTuples(ncells)
        for i in xrange(ncells):
            res.SetTuple1(i,float(domain_id))
        return res

py_filter = DomainIdExpr

Normalize Expression

#
# [VisIt Python Expression]
#
# Normalize a scalar variable to [0-1] using (val-min)/(max-min)
#

class PyNormalize(PythonExpression):
    def __init__(self):
        """
        Constructor.
        """
        PythonExpression.__init__(self)
        self.name = "PyNormalize"
        self.description = "Normalize Expression"
    def execute(self,data_sets, domain_ids):
        data_arrays =  [ ds_in.GetCellData().GetArray(self.input_var_names[0]) for ds_in in data_sets]
        vmin, vmax = self.__cal_min_max(data_arrays)
        vmin = mpicom.min(vmin)
        vmax = mpicom.max(vmax)
        res_arrays = self.__normalize(data_arrays,vmin,vmax)
        return self.__construct_result_dsets(data_sets,domain_ids,res_arrays)
    def __cal_min_max(self,data_arrays):
        vmin, vmax = 1e100,-1e100
        for darr in data_arrays:
            nvals = darr.GetNumberOfTuples()
            vals = [ darr.GetTuple1(i) for i in xrange(nvals)]
            vmin = min(vmin,min(vals))
            vmax = max(vmax,max(vals))
        return vmin, vmax
    def __normalize(self,data_arrays,vmin,vmax):
        res = []
        vdiff = vmax - vmin
        if vdiff == 0: vdiff = 1
        for darr in data_arrays:
            nvals = darr.GetNumberOfTuples()
            vals = [ darr.GetTuple1(i) for i in xrange(nvals)]
            carr = vtk.vtkFloatArray()
            carr.SetNumberOfComponents(1)
            carr.SetNumberOfTuples(nvals)
            for i in xrange(nvals):
                carr.SetTuple1(i,(vals[i]-vmin)/vdiff)
            res.append(carr)
        return res
    def __construct_result_dsets(self,data_sets,domain_ids,res_arrays):
        res_sets = []
        for i in xrange(len(data_sets)):
            res_arrays[i].SetName(self.output_var_name)
            rset = data_sets[i].NewInstance()
            rset.ShallowCopy(data_sets[i])
            rset.GetCellData().AddArray(res_arrays[i])
            res_sets.append(rset)
        return res_sets,domain_ids

py_filter = PyNormalize

Zone Center Expression

This expression computes a vector that is the average of all of the nodes in a zone, giving a zone center. Create a new expression called "zonecenter" and set its output type to vector expression. Click on the Python tab and pass in the name of the mesh on which to operate in the arguments. Then paste this script into the Python Expression Script editor.

class ZoneCenterExpr(SimplePythonExpression):
    def __init__(self):
        SimplePythonExpression.__init__(self)
        self.name = "ZoneCenterExpr"
        self.description = "Zone center"
        self.output_is_point_var = False
        self.output_dimension = 3
    def derive_variable(self,ds_in,domain_id):
        ncells = ds_in.GetNumberOfCells()
        res = vtk.vtkFloatArray()
        res.SetNumberOfComponents(3)
        res.SetNumberOfTuples(ncells)
        for i in xrange(ncells):
            cell = ds_in.GetCell(i)
            npts = cell.GetNumberOfPoints()
            inv_npts = 1. / npts
            x,y,z = 0,0,0
            for j in xrange(npts):
                pt = ds_in.GetPoint(cell.GetPointId(j))
                x += pt[0]
                y += pt[1]
                z += pt[2]
            x *= inv_npts
            y *= inv_npts
            z *= inv_npts
            res.SetTuple(i, (x,y,z))
        return res
 
py_filter = ZoneCenterExpr

If you want just a single component of the zone center, you can create another expression to obtain it. Create another expression X and set its value to:

zonecenter[0]

ZoneCenterExpression.png

Calling from Python

If you are defining a Python filter in a VisIt CLI Python script, you can save your Python expression script to an external file and then get VisIt to load it via the DefinePythonExpression function.

meshName="mesh"
DefinePythonExpression("zonecenter",file="zonecenter.pye",args=[meshName],type="vector")

Connect the dots

If you have a set of particles represented by points, you can connect them with lines using the following Python expression.

Known issues:

  • The resulting expression will plot as white with the Pseudocolor plot.
    • This can be corrected by copying over a variable with same name as the expression. The input_var_names[0] value is the name of the mesh. Can either add another variable input, in which case we can use input_var_names[1] as the name of the variable, or we can hard-code a name like "connect".
class ConnectTheDots(PythonExpression):
    def __init__(self):
        """
        Constructor.
        """
        PythonExpression.__init__(self)
        self.name = "ConnectTheDots"
        self.description = "Connect the dots Expression"

    def execute(self,data_sets, domain_ids):
        grids = []
        for i in xrange(len(data_sets)):
            ugrid = vtk.vtkUnstructuredGrid()
            ugrid.SetPoints(data_sets[i].GetPoints())
            npts = ugrid.GetPoints().GetNumberOfPoints()
            ugrid.Allocate(npts-1)
           # Connect the points into lines.
            for j in xrange(npts-1):
                ids = vtk.vtkIdList()
                ids.InsertNextId(j)
                ids.InsertNextId(j+1)
                ugrid.InsertNextCell(3, ids)
            # Try and copy over variables, field data
            d = data_sets[i].GetPointData().GetArray(self.input_var_names[0]) 
            ugrid.GetPointData().AddArray(d)
            ugrid.GetPointData().SetScalars(d)
            ugrid.GetPointData().SetActiveScalars(self.input_var_names[0])
            ugrid.GetFieldData().DeepCopy(data_sets[i].GetFieldData())
            grids = grids + [ugrid]
        return grids, domain_ids
 
py_filter = ConnectTheDots

Unique Points

Sometimes you get unstructured meshes that have duplicated points and certain algorithms in VisIt don't work well with that. You instead want to combine the like points into single, unique points and then change the connectivity of the mesh to use the unique point list. This is a Python expression that does just that:

  • NOTE: the filter does fix up the connectivity but it is not successfully getting the point variables set up as scalars on the output mesh. This can probably be fixed -- but do it later...
class UniquePoints(PythonExpression):
    def __init__(self):
        """
        Constructor.
        """
        PythonExpression.__init__(self)
        self.name = "UniquePoints"
        self.description = "Combine unique points expression"
 
    def execute(self,data_sets, domain_ids):
        grids = []
        for i in xrange(len(data_sets)):
            # Make unique points.
            npts = data_sets[i].GetNumberOfPoints()
            pts = {}
            old2new = []
            new2old = []
            newpts = []
            pIdx = 0
            for j in xrange(npts):
                p = data_sets[i].GetPoints().GetPoint(j)
                key = "%g,%g,%g" % p
                if key not in pts:
                    pts[key] = pIdx
                    new2old = new2old + [pIdx]
                    newpts = newpts + list(p)
                    pIdx = pIdx + 1
                old2new = old2new + [pts[key]]
            newvtkpts = vtk.vtkPoints()
            newnpts = len(newpts)/3
            for j in xrange(newnpts):
               newvtkpts.InsertNextPoint(newpts[j*3], newpts[j*3+1], newpts[j*3+2])

            # Make a new unstructured grid using the new points.
            ugrid = vtk.vtkUnstructuredGrid()
            ugrid.SetPoints(newvtkpts)
            npts = newvtkpts.GetNumberOfPoints()
            ugrid.Allocate(npts*5)
            for j in xrange(data_sets[i].GetNumberOfCells()):
                cell = data_sets[i].GetCell(j)
                cellids = cell.GetPointIds()
                ids = vtk.vtkIdList()
                ncellpts = cellids.GetNumberOfIds()
                for k in xrange(ncellpts):
                    ids.InsertNextId(old2new[cellids.GetId(k)])
                ugrid.InsertNextCell(cell.GetCellType(), ids)

            # Try and copy over variables, field data
            # (FIXME: this code gets the point values we care about and puts them
            #         in new arrays. I'm not getting the scalar BS right though.
            #         Debug more later...
            # )
            varidx = 0
            nvars = data_sets[i].GetPointData().GetNumberOfArrays()
            for varidx in xrange(nvars):
                orig = data_sets[i].GetPointData().GetArray(varidx)
                newvar = orig.NewInstance()
                newvar.SetName(orig.GetName())
                newvar.SetNumberOfComponents(orig.GetNumberOfComponents())
                newvar.SetNumberOfTuples(newnpts)
                for k in xrange(newnpts):
                    newvar.SetTuple(k, orig.GetTuple(new2old[k]))
                ugrid.GetPointData().AddArray(newvar)
                ugrid.GetPointData().SetScalars(newvar)
            ugrid.GetPointData().SetActiveScalars(self.input_var_names[0])
            ugrid.GetFieldData().DeepCopy(data_sets[i].GetFieldData())
            ugrid.GetCellData().DeepCopy(data_sets[i].GetCellData())

            grids = grids + [ugrid]
        return grids, domain_ids
 
py_filter = UniquePoints

Line Intersection Test

class LineIntersectionExpr(SimplePythonExpression):
    def __init__(self):
        SimplePythonExpression.__init__(self)
        self.name = "LineIntersectionExpr"
        self.description = "Identify cells that intersect a line"
        self.output_is_point_var  = False
        self.output_dimension = 1
    def modify_contract(self,contract):
        pass
    def derive_variable(self,ds_in,domain_id):
        # ds_in is a vtk dataset, we want
        # to create and return a new vtkDataArray
        # that contains the result
        cell_vals = ds_in.GetCellData().GetArray(self.input_var_names[0])
        ncells = ds_in.GetNumberOfCells()
        res = vtk.vtkFloatArray()
        res.SetNumberOfComponents(1)
        res.SetNumberOfTuples(ncells)
        # Here is our line start and end points
        p1 = [0.5, .5, -1.0]
        p2 = [0.5, .5, 1.0]
        tolerance = 0.001;
        for i in xrange(ncells):
            cell = ds_in.GetCell(i)
            t = vtk.mutable(0)
            x = [0.0, 0.0, 0.0]
            pcoords = [0.0, 0.0, 0.0]
            sub_id = vtk.mutable(0)
            isect_res = cell.IntersectWithLine(p1, p2, tolerance, t, x, pcoords, sub_id);
            res.SetTuple1(i,isect_res)
        return res

py_filter = LineIntersectionExpr

Queries

Export Cell Values to OKC

# file:visit_cell_export.vpq
# Export Cells to OKC Point Dataset
#
class CellExportQuery(SimplePythonQuery):
    def __init__(self):
        SimplePythonQuery.__init__(self)
        self.name = "CellExportQuery"
        self.description = "OKC Cell Export Query"
        self.obase = "cell.export.output"
    def pre_execute(self):
        self.active_domains = []
    def execute_chunk(self,ds_in,domain_id):
        ncells = ds_in.GetNumberOfCells()
        if ncells > 0:
            nvars = len(self.input_var_names)
            self.active_domains.append(domain_id)
            var_arrays = [ ds_in.GetCellData().GetArray(vname) for vname in self.input_var_names]
            f = open("%s.%04d.okc" % (self.obase,domain_id),"w")
            f.write("%d %d 1\n" % (3+nvars,ncells))
            f.write("cell_x\n")
            f.write("cell_y\n")
            f.write("cell_z\n")
            for vname in self.input_var_names:
                f.write("%s\n" % vname)
            for i in xrange(3+nvars):
                f.write("-1e30 1e30 10\n") # we could provide the actual extents here ...
            for i in xrange(ncells):
                cell = ds_in.GetCell(i)
                bounds = cell.GetBounds()
                cx = (bounds[1] + bounds[0])/2.0
                cy = (bounds[3] + bounds[2])/2.0
                cz = (bounds[5] + bounds[4])/2.0
                f.write("%s\t%s\t%s" % (str(cx),str(cy),str(cz)))
                for var_array in var_arrays:
                    f.write("\t%s" % str(var_array.GetTuple1(i)))
                f.write("\n")
            f.close()
    def post_execute(self):
        # create .visit file to group domains
        if mpicom.parallel:
            # collective on domains list
            res = mpicom.gather(self.active_domains)
            self.active_domains = []
            for r in res:
                self.active_domains.extend(r)
        if mpicom.rank() == 0:
            f = open(self.obase + ".visit","w")
            f.write("!NBLOCKS %d\n" % len(self.active_domains))
            for domain_id in self.active_domains:
                f.write("%s.%04d.okc\n" % (self.obase,domain_id))
            f.close()
            self.set_result_text("OKC Cell Export Complete. (%s.visit)" % self.obase)

py_filter = CellExportQuery

Want to take it for a spin? Try this:

  • Save the above script as 'visit_cell_export.vpq'
  • Open 'multi_rect2d.silo'
  • Create a Pseudocolor Plot of 'd'
  • From the cli type:
PythonQuery(file="visit_cell_export.vpq",vars=["default","p"])

Logical Extents Query

For DB Plugins that provide "base_index" field data:

#
# Simple Query to Calculate the Logical Extents
#

class LogicalExtentsQuery(SimplePythonQuery):
    def __init__(self):
        SimplePythonQuery.__init__(self)
        self.name = "LogicalExtentsQuery"
        self.description = "Logical Extents Query"
        self.msg = ""
    def pre_execute(self):
        self.ijk = [-1,-1,-1]
    def execute_chunk(self,ds_in,domain_id):
        bidx_array = ds_in.GetFieldData().GetArray("base_index")
        bidx = [bidx_array.GetComponent(0,i) for i in range(3)]
        dims = ds_in.GetDimensions()
        self.ijk[0] = max(self.ijk[0],dims[0] + bidx[0])
        self.ijk[1] = max(self.ijk[1],dims[1] + bidx[1])
        self.ijk[2] = max(self.ijk[2],dims[2] + bidx[2])
    def post_execute(self):
        self.ijk = mpicom.max(self.ijk)
        if mpicom.rank() == 0:
            self.set_result_text(self.msg + "\nLogical Extents =  %s " % str(self.ijk))
            self.set_result_value(self.ijk)

py_filter = LogicalExtentsQuery