1
$\begingroup$

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:

enter image description here

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:

enter image description here

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. enter image description here

$\endgroup$
  • $\begingroup$ I just tried to render a VTU file in ParaView 5.4.0 64 bits in Linux Mint 18.2 and it worked, you can see the image here. $\endgroup$ – nicoguaro Aug 16 '17 at 17:30
  • $\begingroup$ It seems that your image link is not working. Also, I can't really see how you built the VTK file from that XML file. Did you save each cell as a list of PolyVertex objects? Did you save edges or connectivity? The source code might be more helpful... I've updated the question to include it. $\endgroup$ – user3482876 Aug 16 '17 at 17:44
  • $\begingroup$ Weird thing about the image. Here it is again. $\endgroup$ – nicoguaro Aug 16 '17 at 17:53
  • $\begingroup$ Regarding the file creation, I did not write it myself, I just exported an unstructured mesh from ParaView for tests purposes. $\endgroup$ – nicoguaro Aug 16 '17 at 17:54
  • 2
    $\begingroup$ 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 integration point data is normally dealt with is to extrapolate it to the element nodes. Then you can average the values from the different elements connecting to each node. Alternatively, you can define multiple nodes at the same location (one for each element) to represent the discontinuous results (this is the approach used by Deal II, for example). $\endgroup$ – Bill Greene Aug 16 '17 at 18:20
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.