Add Discontinuous Galerkin Finite Element Mesh Support to VTK/ParaView
As part of an ongoing collaboration between Sandia National Lab and Kitware, it is desired to upgrade VTK and ParaView to fully support Discontinuous Galerkin parametrizations read in by the Ioss reader.
Project Requirements
Phase I
-
Upgrade Ioss reader to support prototype DG Exodus mesh -
The info records syntax has changed. Adapt parsing methods to follow the latest pattern.
-
-
Support for H^Grad elements -
Linear bases (hex, quad, tri elements) -
Quadratic bases (hex, quad, tri elements) -
Higher order bases (hex, quad, tri elements)
-
-
Support for H^Div elements (interpolation onto Lagrange elements of proper order) -
Linear bases (hex, quad, tri elements) -
Quadratic bases (hex, quad, tri elements) -
Higher order bases (hex, quad, tri elements)
-
-
Support for H^Curl elements (interpolation onto Lagrange elements of proper order) -
Linear bases (hex, quad, tri elements) -
Quadratic bases (hex, quad, tri elements) -
Higher order bases (hex, quad, tri elements)
-
Phase II
-
Upgrade rendering support for exotic element types (above simple interpolation to Lagrange elements) -
Upgrade rendering hidden face removal for disconnected DG elements -
Upgrade VTK data model to support mismatched field order and mesh order
Background
VTK's Datamodel Challenge
In a classical Galerkin finite element mesh, control variables defining solution fields on the mesh are often associated with physical features of the mesh. For instance, in a linear element, control variables will be associated with the vertices of the element. Thus, when multiple elements share a physical feature, they also share control variables associated with those physical features. This ensures continuity of solution fields across elements. It is this collection of assumptions that forms the foundation of VTK's data model. No matter what elements make up a mesh, each element is made up of vertices, and the connectivity of vertices for the mesh is identical to any connectivity of control variables defining solution fields on the mesh.
This model works well for many finite element applications. However, it falls apart in a number of areas. First, for a simple geometry, it is often advantageous to use a linear mesh to define the domain, but it may be advantageous or necessary to use higher order basis functions to define solution fields, leading to a mismatch in the connectivity of the mesh and solution fields. Second, some problems require different basis functions to be used for different fields. A popular example is the stable velocity-pressure pair where pressure is described by a field one polynomial degree lower than the velocity. Here, two fields will require separate connectivity arrays to organize control variables across elements which may both be different from that describing the physical mesh. Third, some basis functions have little to no relationship with geometry at all. Modal, or spectral elements have control variables that have no physical distribution at all, only describing weights on basis functions over an element. This problem is complicated further by basis functions which take the form of vector fields over an element which could be represented by different polynomial bases in each cardinal direction. The basis function rabbit hole is remarkably deep!
In VTK's current form, the only way to visualize fields constructed by exotic basis functions is to covert to something which satisfies all of The VTK requirements above. For instance, one can always interpolate onto Lagrange polynomial elements of an appropriate order and mosey along with visualization. This can be time consuming for a user to do this themselves, and may be exceedingly wasteful with memory and computational power depending on the size of the user's problem and their particular requirements.
In order to begin addressing this, the first step will be to disconnect the mesh connectivity array from fields in VTK. This could be done by allowing for an arbitrary collection of connectivity arrays, and then assigning one of them to the geometry and the rest to fields. The second key will be to expand how shaders utilize these connectivities and draw arbitrary order fields on mismatched order meshes. For more exotic fields such as those defined by vector-valued basis functions, new shaders must be introduced to properly render.
Discontinuous Galerkin Meshes
Moving one step away from classical Galerkin methods, it can be desirable to break that cross-element continuity, which is where Discontinuous Galerkin comes in. Here, each element has its own collection of control variables, sharing none with neighboring elements. This immediately runs into the challenge in VTK where control variables do not align with physical components of the mesh (which in this case are typically shared across elements). (Sandia: if you would like to expand this section for other readers, that would be great!)
Current Exodus Implementation
Currently, in order to encode discontinuous Galerkin (DG) fields on a mesh, the Exodus file format employs a novel hack. In the Exodus file, a DG field is written as a multi-component cell-centered field where each component represents an element control variable. For instance, a field described by linear H^Grad DG elements over a quad mesh will encode each scalar field by a 4-component cell-centered vector field. A field described by a quadratic H^Grad DG element over a hex mesh will encode each scalar field by a 27-component cell-centered vector field. For vector fields, each component of the vector field is described as a scalar field for the purpose of writing. Metadata for each DG field is then encoded in the info-records section of the Exodus file.
Current Ioss Reader Implementation
At the time of this writing, the Ioss reader has been upgraded to handle the initial prototype Exodus DG file format. This is done via a few modifications. First, the info-records section of the Exodus file is parsed into a data structure which stores information about what fields on which blocks are DG fields. Second, on DG blocks of elements, the mesh is expanded so that no element shares vertices with any other element. This is done so that the eventual DG fields share the same connectivity with the mesh elements. Second, if the DG fields are quadratic, the mesh is upgraded to a quadratic mesh from the linear mesh stored in the Exodus file. This is done so that on each element, the number of vertices match the number of control variables associated with that element. We do this by interpolating extra vertices on edges, faces and the volume center of each element. Next, the routine GetFields
was modified to filter out DG fields, and a new GetDGFields
routine was added to properly parse the DG fields.