Volume rendering 3D VTK *.vtu UnstructuredGrid file in Paraview

Question

I would like to volumetrically render 3D scalar data in Paraview, and I'm not sure if my inability to do so is incorrect usage of VTK or Paraview.

I have built a *.vtu VTK unstructured grid file containing 2 cells consisting of 10 PolyVertex objects with scalar data associated with each point. These are meant to represent points where the value of the numerical solution is known in space (quadrature points for a finite element, for example). I can load the file into Paraview and view the output as points:

No problem. However, when I choose the "Volume" representation, rather than points, the points simply disappear, without any volume rendering. I am looking for paraview to linearly interpolate the solution between the points for each cell, the way it would if I provide an Image (Uniform Rectilinear Grid) from a data file:

However, I seem to be unable to find documentation for how to do this. I imagine the finite element community must commonly render unstructured volume data, so this is surprising.

And the source code used to write out the VTK file in Python:

class VtkPolyVertCloud(object):
    """ save each finite element as a set of polyvertices, but lose cell information """

    def __init__(self):

        # geometry
        self.points= vtk.vtkPoints()
        self.grid = vtk.vtkUnstructuredGrid()

        # data
        self.values = vtk.vtkDoubleArray()
        self.values.SetName('point_values_array')

        self.grid.SetPoints(self.points)
        self.grid.GetPointData().SetScalars(self.values)

    def add_polyVertex_cell(self, points, data):
        """
        adds points according to user-supplied numpy arrays, for convenience and to eliminate loops
        in calling code

        @param points: numpy array of 3d point coords -- points.shape = (npoints, 3)
        @param data: scalar-valued data belonging to each point -- data.shape = (npoints,)
        """
        npts = points.shape[0]
        assert(points.shape[1] == 3)             # make sure 3d points passed in
        assert(data.shape[0] == npts) # make sure same number of data, points

        pv = vtk.vtkPolyVertex()
        pv.GetPointIds().SetNumberOfIds(npts)
        for idx, point in enumerate(points):
            pointID = self.points.InsertNextPoint(point)
            pv.GetPointIds().SetId(idx, pointID)
            self.values.InsertNextValue(data[idx])

        self.grid.InsertNextCell(pv.GetCellType(), pv.GetPointIds())

and the calling code:

def test_vtkPolyVertexCloud_writeToFile():
    """ adds a set of polyvertices meant to represent a finite element """
    pc = vtku.VtkPolyVertCloud()
    points, data = get_random_points_and_data(10)
    pc.add_polyVertex_cell(points, data)
    pc.add_polyVertex_cell(points + 1, data)

    # write
    fn = 'test_PolyVertexCloud.vtu'
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName(fn)
    writer.SetInputData(pc.grid)
    writer.Write()

Update: I took heed of the accepted answer below and did the following: 1. Performed a spatial Delaunay triangulation on each of my finite elements (the numerical solution is known at the nodes of each finite element). The triangulation is fast since there aren't that many points even for a high-order finite element. 2. Constructed a VTK file where each cell is a tetrahedron from the spatial Delaunay triangulation on each element.

Paraview is able to volumetrically plot this.

Answer 1

It is not surprising that VTK can't generate a volume rendering for the poly_vertex cell type since there is no topology associated with that type of cell.

The way finite element integration point results is often dealt with is described below:

It is common in finite element methods to calculate results (e.g. stresses) at the element integration points. But, typically, for plotting purposes, the result quantities are needed at the nodes. One way to achieve this is, for each element, to fit a parametric function to the integration point data. Then this function can be evaluated at each of the element nodes. The exact form of this function depends on the particular element topology and the number of integration points but it is often similar to the shape functions used for the element.

In general, since this procedure is applied element-by-element, the result function will be discontinuous when more than one element is connected to a node. Most visualization libraries (e.g. VTK) don't directly support discontinuous nodal results so one of the following approaches is typically used.

1. You can simply average the element-nodal values for all the elements attached to a particular node.

2. You can create "extra" nodes for each element attached to a particular node. These nodes are at the same location as the actual node. But having duplicate nodes allows each element to have its own, unique set nodes to which results can be attached. This is the approach used by Deal II, for example.