Discontinuous Galerkin elements and other novel cell-types/function-spaces

Hi all! We (@jaswantp @c.wetterer-nelson and me) have been working on a prototype to support discontinuous Galerkin (DG) finite elements (and hopefully many other novel function spaces such as mimetic techniques, X-FEM, CW-complexes, etc.). This post describes the problem and discusses a work-in-progress merge request with our prototype solution.


Figure 1. A screenshot of our prototype mapper rendering two DG hexahedra that share 4 corner points but still have field discontinuities along the shared boundary. This is achieved with a custom mapper that delegates most of the work to “render responders” for each cell type.

Problem Statement

As many of you are aware, vtkDataObject subclasses make increasingly-strong assumptions about the nature of data each can process. In particular, vtkDataSet is the first class in the inheritance hierarchy that introduces the notion of data with spatial extents. However, this is not the only assumption present. Below we outline some additional (and problematic for DG elements) assumptions.

Fixed set of cell shapes

vtkDataSet has a strict enumeration of cell types; external libraries cannot add support for new cell types to VTK’s set of types. This has led to a proliferation of cell types in VTK with some overlap between them.

We would like a data-object that accepts cells defined external to VTK.

1–1 relationship between arrays and fields

vtkDataSet adds vtkPointData and vtkCellData array containers and assumes that arrays in these containers have a 1-to-1 mapping to spatial functions (i.e., mathematical maps from locations in space to scalar, vector, or tensor values) which we will call attributes or fields. (Note that sometimes the arrays in the vtkFieldData instance inherited from vtkDataObject are treated as functions constant over all of space).

DG elements don’t fit the pattern of these two (three including vtkFieldData) specific array containers, which we’ll discuss next. However, even if DG elements could be adapted to the existing function types (which the vtkFiniteElementFieldDistributor does to some extent by duplicating points on shared boundaries), there are use cases where this correspondence between arrays and fields fails. If DG cells are expected to each have a different polynomial degree (for the same field), then multiple arrays are required (to specify the degree per element + per field and possibly additional connectivity entries). Also, there are other use cases (e.g., NURBS elements that have knot vectors in addition to values defined on control polygons), where we may wish to use multiple arrays to specify a single function in space.

Only two or three types of fields are supported

The assumption that arrays in vtkFieldData, vtkCellData, and vtkPointData array containers define

  • spatially constant fields over all of space covered by cells;
  • spatially-constant fields within each cell as determined by a single array value for each cell (cell-data); or
  • spatially varying fields within each cell as determined by array values specified at “corner” points of cells (point-data)

respectively, does not match many finite element schemes. DG finite elements define fields that vary spatially within a cell but are not specified by values shared at boundaries (unlike vtkPointData arrays). Many codes use coefficients that are associated with edges or faces as opposed to cell corners. Some use generic degrees of freedom that are not associated with spatial locations or cell boundaries at all.

Field values must come from the same space as array values

Although not relevant to Lagrangian DG elements, some DG elements use function spaces such as H-curl. In these function spaces, the coefficients stored with each finite element are single (scalar) values, but the field evaluates to a vector or tensor value at each location in space. (The basis functions are vector/tensor-valued, so weighting them by scalars unique to each finite element produces a vector/tensor field using only scalar coefficients.)

We would like to specify the range of the spatial function separately from the array(s) used to defined the function.

Fixed set of visualization operations

Because vtkDataSet also provides an API to fetch cells into a uniform container (vtkCell and its subclasses), cells provide support for a basic set of visualization “operations” (isocontouring, cutting, clipping, evaluation, parameter location, boundary extraction). Other capabilities must be provided by filter algorithms that explicitly handle each cell type (or more frequently, complain when they encounter an unknown cell type).

We would like to separate the cell type from visualization operations that act on cells of that type.

Proposed solution

We are providing a prototype solution for comment that makes none of these assumptions. It provides a new subclass of vtkDataObject called vtkCellGrid. This subclass holds:

  • An arbitrary number of array containers rather than just two. Each container requires its arrays to have the same number of tuples but does not assume each container corresponds to a type of field.

  • vtkCellAttribute instances that have a 1-to-1 correspondence with mathematical fields (or attributes or functions) and may reference any number of arrays from any number of containers above.

    One special vtkCellAttribute instance must be provided to map cell reference-coordinates into world coordinates; this is called the physical or shape attribute since it maps the cell into physical space. This field must provide vector values at each parametric coordinate.

  • Subclasses of vtkCellMetadata. Each subclass corresponds to a different type of cell (hexahedral DG cell, tetrahedral NURBS cell, etc.).


Figure 2. An example of how discontinuous Galerkin elements might be supported. The light green rectangles on the left are array containers that hold dark green arrays of matching size. The light blue rectangles on the right are vtkCellAttribute instances corresponding to fields that reference arrays from the left. At top: the case where each field has the same interpolation order for all cells. At bottom: the case where each cell may have a different polynomial order.

Since one instance of any vtkCellMetadata-subclass is accepted, the vtkCellGrid can be extended by external libraries that add cells of new types. However, we also want external libraries to be able to extend VTK’s cell types by adding support for new visualization operations. To support this, we introduce vtkCellGridQuery, which is an object passed to a vtkCellGrid in order to perform an operation on (or obtain information about) cells of every type in the vtkCellGrid. External libraries then register instances of a vtkCellGridResponder that can respond to a particular query type for a particular cell type. For example:

vtkCellGrid* grid;
vtkNew<vtkCellGridBoundsQuery> boundsQuery;
if (grid->Query(boundsQuery))
{
  vtkBoundingBox bbox(boundsQuery->GetBounds());
}
else
{
  std::cerr << "Some cell type could not respond to this query.\n";
}

This way, a library that adds a new cell type can add responders for queries added by a different external library.

Prototype results

We have implemented the data-object classes described above plus the following in order to validate the approach:

  • vtkCellGridReader – a simple JSON-based reader that constructs instances of vtkCellGrid for testing purposes.
  • vtkCellGridBoundsQuery – A query to compute extents of all cells of all types in a vtkCellGrid. A responder for discontinuous Galerkin cells is provided (not specific to any particular cell shape). The query is used by vtkCellGrid internally to compute bounds for the grid.
  • vtkCellGridSidesQuery – A query to compute external sides (i.e., external faces of volumetric cells; external edges of surface cells) of cells in a vtkCellGrid. A responder for discontinuous Galerkin cells is provided (not specific to any particular cell shape). A filter named vtkCellGridExtractSurface uses this query to add an array of sides to an input cell-grid. The array is used by our prototype renderer below.
  • vtkOpenGLCellGridRenderRequest – A query to render external faces of a vtkCellGrid to an OpenGL renderer. This query is used by a cusom vtkMapper subclass named vtkOpenGLCellGridMapper. A responder for discontinuous Galerkin hexahedra (named vtkDGHexahedronOpenGLRenderer) is provided.

You can try the prototype out by checking out VTK merge request 9499.

  1. Build the branch of VTK from MR 9499.

  2. Copy the included test-prototype-dg-mapper.py file and test data (dgHex2d.dg) from the top of the source tree to your build tree.

  3. Run ./bin/vtkpython ./test-prototype-dg-mapper.py. Many debugging printouts are included. The following keys demonstrate different features of the prototype python demo:

    Key Action
    Ctrl+d Cycle through rendering L-1 and L-2 norms of parametric coordinates instead of coloring the hexahedra by field value.
    Ctrl+b Cycle through rendering of basis functions instead of coloring the hexahedra by field value.
    Ctrl+f Cycle through coloring the hexahedra by different example fields. One is continuous at the shared boundary, but 3 others are not.
    r Reset the camera as well as the color-by mode (to color by the currently-selected scalar field).



Figure 3. Cycling through parametric coordinates using the Ctrl+D key. The top row shows the r, s, and t coordinates, respectively. The bottom row shows the L-2 norm of the r+s coordinates and s+t coordinates, respectively (the t+r norm is solid red for the viewpoint in the screenshots).


Figure 4. Cycling through 4 of the 8 basis functions for the hexahedron (the other 4 are on the back side of the hexahedra from this perspective).

Questions

While the prototype is very minimal, we hope there is enough there to spur some discussion of what works and what doesn’t.

  1. Are there cell types you think this approach would not handle well? What assumptions does the prototype make that are unhelpful?
  2. The prototype does not have a lot of traditional functionality. Rather than comment on what is missing, please comment on how you think missing features should be approached.
  3. Although we did not have time to pursue it with this funding, we think this approach would accommodate NURBS cells more cleanly that the existing vtkBezierXXX cells (e.g., vtkBezierHexahedron) since it does not preclude 4-dimensional points and knot-vector arrays tied to the “shape” field. Do you have a preference on how 4-dimensional (x,y,z,w) point-coordinates should be handled? Do you see issues with NURBS this approach cannot handle?

We hope to continue development along these lines in the future, so your feedback is really appreciated.

3 Likes

You’ve been busy :slight_smile: this looks promising.

Before I actually dive into the implementation details, a couple of pragmatic questions come to mind:

  • I see that the code is added into the existing VTK modules (e.g., Common/DataModel). I would consider having compile-time flags to enable/disable inclusion of these DG classes, and/or segregate them into their own module(s). (We spent a fair amount of time modularizing VTK etc. - VTK continues to grow in size and it is very important that users can enable/disable code that they want/don’t want.)
  • It would be nice to see a few more algorithms implemented, for example isocontouring, probing, and/or streamlines (did I miss these in your MR?). It would help the community understand the practical usage of these classes. So a suggestion for addressing missing features: identify 4-5 important algorithms and then implement them (at a very rudimentary level), extending/adjusting APIs as necessary.
1 Like

Hi Will,

  • While I might argue to keep the base vtkCellGrid class in Common/DataModel, I agree it makes sense to move the other pieces into a couple other modules (maybe Filters/CellGrid and Rendering/CellGrid?).
  • I agree, having more algorithms would be nice. Our funding for this was specifically for handling the visual aspects of discontinuous cells with conforming boundaries, but
    • Isocontour is interesting because of discontinuities in DG cells – would you leave cracks at cell boundaries?
    • Streamline and/or Probe would take advantage of conformal shapes (shared corner points allowing neighbor traversal) that our current approach to DG cells cannot provide;
    • Filters to convert to/from unstructured grids (with accompanying conversion of fields to point or cell fields). This would make comparisons with existing algorithms easier to benchmark.

Besides filters, there are some other directions to push things that would be useful for validating the approach

  • Adding support for a few more cell shapes (DG quads, tetrahedra, and triangles – probably in that order)
  • Adding support for some other exotic cell types (NURBS patches, polyhedral cells) and field types (H-curl/H-div function spaces, algebraic splines) to ensure we stay general enough to handle everything.
  • Getting the mapper to deal with selection (rendering cell IDs, which may best be implicit).

The other thing I like about this approach is that the mapper pushes the raw arrays for the grid to the GPU (or, could allow GPU-based simulations to keep them there) rather than converting the volume mesh to a surface mesh for rendering. This should pave the way for fancier widgets for selecting volumetric cells, since they might modify the “sides” array used by the mapper to include internal cell boundaries on a selection curve/surface. Some prototypes that use the mapper this way might also help motivate the approach.

This is very interesting work and it seems very tricky to design a satisfying solution with all of VTK’s legacy. Hopefully this will work well :slight_smile:

I have a few questions:

  • How does vtkCellGrid behave in an SMP instance? Do you have a thread-safe way to get a cell? If not, could we fill a CellType passed as a parameter?
  • Creating cells is pretty expensive: lots of buffer copying. As a result, performance can be hit quite a lot from actually using vtkCell instances. An example of that is vtkClipDataSet, which uses vtkCell::Clip to generate the output cells, and is terribly slow compared to vtkTableBasedClipDataSet which uses switch statements on the input cell types in coordination with a clip table to decide how to clip the cells. I feel like the latter approach would not be possible with this design. Do you think there could be a way to optimize a clip filter for vtkCellGrid in a similar way somehow, or is it a necessary limitation that we are paying by linking against any custom cell type?
  • Do you think that it would be possible to somehow, in the long run, merge the implementation of every vtkDataSet into some more general API such as vtkCellGrid to avoid having to duplicate the filter implementation for both types? I don’t think it would be achievable in the short term, I’m asking long term, kind of like vtkPartitionedDataSet trying to replace vtkMultiBlockDataSet.
  • Ghost question: how do we handle ghosts in Galerkin meshes? Do we have ghost edges? Ghost anything? If so, is it time to unify the ghost framework? Currently HIDDENCELL and HIDDENPOINT don’t hold the same value for legacy reasons coming from VisIt. Any take on that?
1 Like

Hi Yohann,

Currently, there is no equivalent to vtkCell. The vtkCellMetadata class is a single instance representing all the cells of a given type in the entire vtkCellGrid. I would resist the urge to use class inheritance to force a common API for individual cells. By removing this level of granularity, we can remove a lot of the overhead you mention (copying connectivity/field data/etc. into a data structure just to iterate cells).

As far as SMP goes, the vtkCellGrid doesn’t do anything to help or hurt you; responders to vtkCellGridQuery classes implement an algorithm however they choose (and they should choose to use vtkSMPTools wherever it makes sense).

We could in theory make vtkCellGrid::Query() perform task parallelism by running responders in multiple threads, but that could complicate access to the vtkCellGridQuery instance that all the responders typically insert result/output data into. It would not be impossible to be efficient (for instance, the prototype vtkOpenGLCellGridRenderRequest class provides a GetState<StateType>(cellType) method that each cell type could access without fear of races/deadlocks). However, queries that require responders to modify shared state would need to provide atomics or locks for thread safety.

That is part of the design. Each cell type registers responders that can handle a given query class. When you run a query:

vtkCellGrid* grid;
vtkNew<vtkCellGridClipQuery> clipper;
clipper->SetImplicitFunction(fn);
grid->Query(clipper);

The grid iterates over all its cell types (instances of subclasses of vtkCellMetadata) and for each one, attempts to find a responder (i.e., vtkCellGridResponder<vtkCellGridQuerySubclass>) to the query class you pass it (from the list of registered responders, indexed on both cell type and query type). Each responder is then invoked on the query in turn. No copying is done. The responders can make calls on the query object to record their output. The base query class has virtual Initialize() and Finalize() methods called before-any/after-all responders have been invoked.

Yes, it is possible. Insert coin here. :slight_smile: In theory, even vtkTable or vtkGraph could become vtkCellGrid instances with different subclasses of vtkCellMetadata representing rows/columns or nodes/arcs, respectively.

Right now, there’s enough rope to make plenty of nooses.

That is currrently up to the individual cell classes.

  1. I think traditional DG cells will have slightly different needs than NURBS or polyhedral cells. Exotic cells (table rows? graph arcs?) may have significantly different needs. This has played out a bit with structured vs. unstructured grids in VTK with extents vs pieces (not directly related to ghost cells, but the way that decomposition is performed will affect ghosts and I do not want to force concepts about data decomposition into the base API yet).
  2. If possible, I would like to allow the arrays present in vtkCellGrid to match simulation codes as much as possible rather than force “translation” of ghosts – the query responders should adapt solver-specific array-data to the task at hand.
  3. It is important not to force a uniform API that is not completely generic.

So, I don’t want to deal with ghost cells/points/(edges/faces/graph-arcs/table-rows/table-columns) until we have implemented ghosts in 2-3 specific contexts first.

Thanks for your answer. I think I understand more how this class works now.

In theory, even vtkTable or vtkGraph could become vtkCellGrid instances with different subclasses of vtkCellMetadata representing rows/columns or nodes/arcs, respectively.

I rapidly read through your MR. Correct me if I’m wrong, but it seems that vtkCellGrid is a collection of cell types that you group by cell types. Each cell grid query type is linked to a responder that performs whatever you want done on each visited cells.

I’m trying to understand how this class could take over all of vtkDataSet inheritance tree.

If you want to have, say a vtkImageData, done with this representation, you’ll need to give it a different voxel cell type than on a vtkUnstructuredGrid represented as a vtkCellGrid, because we want to take advantage of the implicit structure of an image and do not need to carry a conn array for the connectivity.

Do you see us specializing instances of vtkCellGrid in such a way that you have a vtkImageData-like or a vtkUnstructuredGrid-like? I feel like it would conflict with vtkCellGrid::AddCellMetadata API for vtkImageData because we’re restricting ourselves to uniform grids in this instance. If we do not create specialized instances of such dataset types, how do you deal with filter selection in ParaView on filters that do not work on a certain dataset type? As an example, we can think of a convolution filter which is well defined on images, but maybe more complicated to generalize to wilder topologies.

Or do we assume more that you can put any cell types in a vtkCellGrid, uniform grids of voxels, triangles, hexahedrons, and that each cell being linked to a responder, we just try to respond to a query, and if no responder exists, we skip, and thus we never gray out filters in ParaView?

Yes. But each cell type may have a different responder for the same query.

To a degree, yes. It is possible that filters will not care about connectivity or other arrays that distinguish cell types (e.g., Calculator-like algorithms). There is nothing in that case to stop you from registering the same responder to multiple cell types. Or it may be that a responder can support several cell types.

For instance, the prototype provides a vtkDGSidesResponder for identifying the external faces (sides) of hexahedra in a mesh. But the same provider could easily handle any shape (it just hashes the connectivity entries for each boundary of a cell and looks for odd numbers of matching hashes to report out). My plan was to add support for another few cell shapes (quad, tet, tri, maybe pyramid) and use the same responder for them all. Coding the faces of each of those types is simple. And if some novel cell type appears that makes things difficult, a separate responder can be added.

Yes, something along those lines. (Edit no: I think I understand your question now. I would add new vtkCellMetadata subclasses but would probably not subclass vtkCellGrid.) I never claimed it was easy to develop, but I do claim it is flexible enough to accommodate nearly all the solver schemas for memory layout you might imagine. Figure 2 in the top post shows how DG cells might be developed first for fields that use the same order interpolant for all cells in a mesh (top) and then adapted to deal with varying the interpolant order on a per-cell basis. You can imagine either (a) abandoning the first cell type for the second or (b) keeping both around to deal with different solvers as efficiently as possible. Which you choose to do depends on (a) how stable the API is (maintenance costs) and (b) how much funding/developer time is available. For other cases (like your image vs unstructured grid), it is easier to see how both would be needed.

You could have a vtkCellGrid instance with both an unstructured grid and a structured grid inside it. (I am not saying that’s advisable, just that you could.) I would expect a NURBS cell type to deal with both rectilinear and curvilinear structured data much like image data today (and be more flexible as well). Many filters would not need to know about the knot vector since they just add new fields arrays on points or cells.

It is not in this prototype, but it would be easy to add a method to vtkCellGrid such as bool TestQuery(vtkCellGridQuery*);. That would simply check whether there is a responder registered for the given query that can handle every vtkCellMetadata object in a particular vtkCellGrid instance.

Not every vtkCellGrid needs to have every type of cell-metadata; also – even if a grid does have every type of cell – it may be that, due to the arrays present, there are zero cells of that type in the grid.

I don’t think we need to do that. We could even have a TestQuery() method provide feedback about what cell types are unsupported for a query; ParaView could then provide tooltips like “The triangle-strip filter does not support quadrilaterals,” for a disabled filter.

@Yohann_Bearzi I realized after responding above that you were asking about subclassing vtkCellGrid rather than subclassing vtkCellMetadata. I don’t see a reason to subclass vtkCellGrid. As my original answer above outlined, I think there are ways to discriminate whether filters apply to a particular grid’s cells without using data-object inheritance.

If you want to get meta/cute, you could even implement TestQuery() as a query (vtkCellGridApplicabilityQuery) that accepts a child query it tests against each cell type. This would let you accumulate feedback for users about compatibility between a filter and the available responders.

@jfausty @Francois_Mazen this may be of interest to you :slight_smile:

1 Like

Hi @dcthomp, @jaswantp and @c.wetterer-nelson,

First, as a developer of DG methods for finite elements in my previous work, I would like to say I am super excited to see this kind of support coming to VTK! It can be a real challenge to visualize solutions to DG type simulations in existing software for exactly all the reasons you guys specified in this thread.

Secondly, I think that the approach of letting users potentially provide their own cells is a good one. There is a whole zoology of element types out there which can differ in the functional spaces they discretize but also how they describe them (positioning of interpolation nodes, modal function definitions, etc.).

Lastly, on a more technical note, in the context of supporting broad p-adaptation, have you though of what would happen for meshes that support different polynomial orders for their geometry as well as their fields ? Indeed, it is possible to generate meshes that are essentially “linear” or straight throughout their domain and higher order only at their boundaries. In order to cut down on the memory cost of the connectivity and the points, you could have different connectivity arrays: one holds elements of order 1 while the others hold elements of higher order geometry, etc. This would operate much in the same way that you’ve designed your support for different order elements for your fields. Of course, the geometric polynomial order and the field polynomial order can still vary independently as long as p_{geom} \le p_{field}.

In any case, thanks for the forum for discussing this topic which looks really well thought out so far!

1 Like

Yes, we just didn’t have enough funding/time to implement this as part of the prototype.

Also, there are many ways that simulations deal with varying order. There is nothing in the prototype that prevents using a different cell type only at boundaries as you outlined, or a different approach such as Figure 2 (but applied to the shape attribute rather than fields used for coloring).

For another example, Szabó and Babuska describe CG (not DG) elements which share Legendre shape functions for higher orders on boundaries. In that case, you could have a single connectivity entry for a shared higher-order face, but it would be two numbers (the face ID and an orientation number tagging a symmetric group that matched the rotation+mirror required to put the face into the volume-element’s reference coordinates). So, if you had a quintic face requiring 4x4 = 16 coeffcients for the face’s degrees of freedom, you would only need 2 numbers in the connectivity array for the volume’s face-reference (the start of the 16 numbers plus the symmetric group ID). This is very compact, but also hard to develop+maintain. It also has coefficients that are associated with boundaries of elements but not with specific locations in space (i.e., Legendre != Lagrange). The approach we are proposing would allow this layout as well, but you have to be willing to spend the time to implement it :slight_smile: .

1 Like