WriteVTK node numbering for 27 node Lagrange Hexahedron

jman87 · January 28, 2023, 3:36pm

I’m trying to export a mesh containing 27-node Lagrange Hexahedron elements using WriteVTK, but I can’t seem to get the node ordering correct.

Another ParaView user suggested I get the ordering from ParaView by going to Sources>Data Objects>Unstructured Cell Types and then generating a Lagrange Hexahedron with cell order of 2, and then turn on the point ID’s using Edit>Find Data. This produces the following:

However, when I export my mesh with WriteVTK using this numbering scheme it produces the following for a single HEX27:

One issue I see is that the numbering in the .vtu file generated by WriteVTK seems to have a different numbering pattern from what ParaView is saying is the order in the other plot above. I’ve tried several permutations to see if I can get it to work but no luck yet.

If someone could confirm if the second-order Lagrange Hexahedron does in fact work for WriteVTK that would be my first request. Then, if it does, if someone could clarify what the proper order of the nodes should be that would be most helpful.

Thanks.

jipolanco · January 28, 2023, 4:18pm

Can you provide a minimal example of your code, to make it easier to identify the issue?

Not sure if you’re already taking this into account, but note that connectivity indices in WriteVTK are one-based (indices start at 1) following the convention in Julia, and as opposed to indexing in VTK/ParaView which is zero-based. So the lower-left corner of your first image would correspond to point id 1 in WriteVTK.

jman87 · January 28, 2023, 5:12pm

As I was posting it I figured that would be your first question. LOL

I’ll see if I can set something simple up. It’s got too many dependencies tangled up in it right now.

I do know about the base 1, which is one of the things I like about Julia. Python still hasn’t learned how to count!

jipolanco · January 28, 2023, 5:58pm

I think the answer is in a relatively recent change in the VTK libraries, since version 9.1. As described in the release notes, the order of the nodes of VTK_LAGRANGE_HEXAHEDRON was changed for consistency with VTK_QUADRATIC_HEXAHEDRON. I’m guessing your initial image uses the new ordering. But the files created by WriteVTK are by default backwards-compatible with previous VTK versions, and so they use the old ordering.

What you can try is to open your .vtu file generated by WriteVTK in a text editor, and change the VTKFile version from 1.0 to 2.2 (in the second line of the file). It should look something like this:

<VTKFile type="UnstructuredGrid" version="2.2" byte_order="LittleEndian" header_type="UInt64">

Then ParaView should recognise the file as having the “new” node ordering, which is the one you used. If this works, I’ll make a patch to WriteVTK to use the latest version (2.2) by default.

jman87 · January 28, 2023, 6:19pm

I just threw the following script together to show what I’m seeing. I haven’t tried your suggestion yet. Here’s my code and the resulting mesh shown in ParaView.

using WriteVTK

points = permutedims(
    [
      1     1     1   
     -1     1     1   
     -1     1    -1   
      1     1    -1   
      1    -1     1   
     -1    -1     1   
     -1    -1    -1   
      1    -1    -1   
      0     1     1   
     -1     1     0   
      0     1    -1   
      1     1     0   
      0    -1     1   
     -1    -1     0   
      0    -1    -1   
      1    -1     0   
      1     0     1   
     -1     0     1   
     -1     0    -1   
      1     0    -1   
      0     0     0   
      0     1     0   
      0    -1     0   
      0     0     1   
     -1     0     0   
      0     0    -1   
      1     0     0   
    ]
)

# My Element Node Numbering
elem = [ 1, 2, 3, 4,        # Vertex   Y =  1
         5, 6, 7, 8,        # Vertex   Y = -1
         9,10,11,12,        # MidEdge  Y =  1
        13,14,15,16,        # MidEdge  Y = -1
        17,18,19,20,        # MidEdge  Y =  0
        21,                 # Body     X=Y=Z=0
        22,23,24,25,26,27]  # Face Centroids

# The following shows the mapping of ParaView's node numbering
# generated via Sources to my node numbering:
#   VTK NID          My NID
#  0  1  3  2  =>  7  8  4  3    # Bottom   Z = -1
#  4  5  7  6  =>  6  5  1  2    # Top      Z =  1
#  8 11 12  9  => 15 20 11 19    # MidEdge  Z = -1
# 22 25 26 23  => 13 17  9 18    # MidEdge  Z =  1
# 13 15 21 19  => 14 16 12 10    # MidEdge  Z =  0
# 16           => 25             # Face     X = -1
# 18           => 27             # Face     X =  1
# 14           => 23             # Face     Y = -1
# 20           => 22             # Face     Y =  1
# 10           => 26             # Face     Z = -1
# 24           => 24             # Face     Z =  1
# 17           => 21             # Body    X=Y=Z = 0

# From what I've output from ParaView, this is what 
# I THINK the permutation for VTK should be:
perm27 = [ 7,  8,  4,  3,
           6,  5,  1,  2,
          15, 20, 11, 19,  
          13, 17,  9, 18,
          14, 16, 12, 10,
          25, 27, 23, 22, 26, 24, 
          21]

# Define cells for WriteVTK (imitating my larger program's logic)
cells = Vector{MeshCell}(undef,1)
for eid = 1:1
    cells[eid] = MeshCell(VTKCellTypes.VTK_LAGRANGE_HEXAHEDRON, elem[perm27,eid])
end

vtk_model = vtk_grid("HEX27_ex", points, cells)

vtk_file = vtk_save(vtk_model)

Resulting Mesh:

jman87 · January 28, 2023, 6:44pm

@jipolanco

I tried changing the version from 1.0 to 2.2 as you suggested but that makes ParaView throw the following error:

FWIW, I’m running ParaView 5.11.0 on a M2 MBP.

Also, for completeness/transparency, I did find a bug in the way I was doing my big codes’ permutation, hence the difference in the images of the mesh. But as the code I posted above shows, there’s still something goofy going on.

Did they change anything else format-wise that could be causing the error from ParaView when switching to the 2.2 format?

jman87 · January 28, 2023, 7:22pm

I think I’ve gotten my issue fixed. Thanks for the help @jipolanco. I was not aware of the 1.0 to 2.2 change. Once I dug into that it made a lot more sense. I’m going to write a summary of the node ordering that seems to be what is expected for 1.0 and for 2.2 to have as a reference and to help anyone else that may have this problem.

jman87 · January 28, 2023, 7:56pm

For generating a 27-node Lagrangian hexahedral element using Julia’s WriteVTK, there are two distinct node orderings available, corresponding to VTKFile versions 1.0 and 2.2, respectively. I have summarized the respective node orders below using the node numbering shown in the following figure, which was generated in ParaView using Sources>Data Objects>Unstructured Cell Types as described above.

Node IDs of a 2nd-order 27-node Lagrange Hexahedron

For VTKFile 1.0 the node numbering is as follows:

 # VTKFile Version 1.0
 0  1  3  2  # Vertex   Z = -1
 4  5  7  6  # Vertex   Z =  1
 8 11 12  9  # MidEdge  Z = -1
22 25 26 23  # MidEdge  Z =  1
13 15 19 21  # MidEdge  Z =  0
16 18        # Face     X = -1,  X = 1
14 20        # Face     Y = -1,  Y = 1
10 24        # Face     Z = -1,  Z = 1
17           # Body     X=Y=Z = 0

For VTKFile 2.2 the node numbering is as follows:

# VTKFile Version 2.2
 0  1  3  2  # Vertex   Z = -1
 4  5  7  6  # Vertex   Z =  1
 8 11 12  9  # MidEdge  Z = -1
22 25 26 23  # MidEdge  Z =  1
13 15 21 19  # MidEdge  Z =  0
16 18        # Face     X = -1,  X = 1
14 20        # Face     Y = -1,  Y = 1
10 24        # Face     Z = -1,  Z = 1
17           # Body     X=Y=Z = 0

Note that the only change in the order is for the last two MidEdge Z=0 nodes, which get swapped. This makes it so that the MidEdge node ordering for version 2.2 is consistent between Z=-1, Z=0, and Z=1, whereas for version 1.0 the Z=0 MidEdge node order was different from Z=+/-1.

As further demonstration of this, I have provided sample VTK files for each of these below.

The following is a mesh of a single HEX27 element for VTKFile 1.0:

<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
  <UnstructuredGrid>
    <Piece NumberOfPoints="27" NumberOfCells="1">
      <PointData>
      </PointData>
      <CellData>
      </CellData>
      <Points>
        <DataArray type="Int64" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="1.7320508075688772">
          1 1 1 -1 1 1
          -1 1 -1 1 1 -1
          1 -1 1 -1 -1 1
          -1 -1 -1 1 -1 -1
          0 1 1 -1 1 0
          0 1 -1 1 1 0
          0 -1 1 -1 -1 0
          0 -1 -1 1 -1 0
          1 0 1 -1 0 1
          -1 0 -1 1 0 -1
          0 0 0 0 1 0
          0 -1 0 0 0 1
          -1 0 0 0 0 -1
          1 0 0
          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0
            </Value>
            <Value index="1">
              1.7320508076
            </Value>
          </InformationKey>
          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0
            </Value>
            <Value index="1">
              1.7320508076
            </Value>
          </InformationKey>
        </DataArray>
      </Points>
      <Cells>
        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="26">
          6 7 3 2 
          5 4 0 1 
          14 19 10 18
          12 16  8 17 
          13 15  9 11 
          24 26 
          22 21
          25 23 
          20
        </DataArray>
        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="27" RangeMax="27">
          27
        </DataArray>
        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="72" RangeMax="72">
          72
        </DataArray>
      </Cells>
    </Piece>
  </UnstructuredGrid>
</VTKFile>

The following is a mesh of a single HEX27 element for VTKFile 2.2:

<VTKFile type="UnstructuredGrid" version="2.2" byte_order="LittleEndian" header_type="UInt64">
  <UnstructuredGrid>
    <Piece NumberOfPoints="27" NumberOfCells="1">
      <PointData>
      </PointData>
      <CellData>
      </CellData>
      <Points>
        <DataArray type="Int64" Name="Points" NumberOfComponents="3" format="ascii" RangeMin="0" RangeMax="1.7320508075688772">
          1 1 1 -1 1 1
          -1 1 -1 1 1 -1
          1 -1 1 -1 -1 1
          -1 -1 -1 1 -1 -1
          0 1 1 -1 1 0
          0 1 -1 1 1 0
          0 -1 1 -1 -1 0
          0 -1 -1 1 -1 0
          1 0 1 -1 0 1
          -1 0 -1 1 0 -1
          0 0 0 0 1 0
          0 -1 0 0 0 1
          -1 0 0 0 0 -1
          1 0 0
          <InformationKey name="L2_NORM_FINITE_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0
            </Value>
            <Value index="1">
              1.7320508076
            </Value>
          </InformationKey>
          <InformationKey name="L2_NORM_RANGE" location="vtkDataArray" length="2">
            <Value index="0">
              0
            </Value>
            <Value index="1">
              1.7320508076
            </Value>
          </InformationKey>
        </DataArray>
      </Points>
      <Cells>
        <DataArray type="Int64" Name="connectivity" format="ascii" RangeMin="0" RangeMax="26">
          6 7 3 2 
          5 4 0 1 
          14 19 10 18
          12 16  8 17 
          13 15 11  9
          24 26 
          22 21
          25 23 
          20
        </DataArray>
        <DataArray type="Int64" Name="offsets" format="ascii" RangeMin="27" RangeMax="27">
          27
        </DataArray>
        <DataArray type="UInt8" Name="types" format="ascii" RangeMin="72" RangeMax="72">
          72
        </DataArray>
      </Cells>
    </Piece>
  </UnstructuredGrid>
</VTKFile>

jipolanco · January 30, 2023, 12:55pm

I just released a new WriteVTK version (v1.17.0) which allows setting the version of the VTK file format.

So in your original script you can just do:

vtk_model = vtk_grid("HEX27_ex", points, cells; vtkversion = :latest)

(or equivalently, vtkversion = v"2.2") to make sure that ParaView uses the newest node ordering. See the updated docs for more details.

jman87 · January 31, 2023, 6:19pm

@jipolanco Is it possible to set the version globally? So only do it once? Or do we need to include it in each statement?

jipolanco · February 1, 2023, 7:21am

I don’t know what you mean by “in each statement”, but the version should be set whenever you call vtk_grid. So basically once for each VTK file that you want to create. There is no global option to set the VTK file version.

juliohm · February 1, 2023, 11:14am

I would like to take the opportunity to ask if you could share a typical pseudocode where these higher order nodes are used. We have plans to add them to our geometries in Meshes.jl but didn’t decide the API yet.

In particular, we would like to introduce all these Quadratic, BiQuadratic, Cubic, and HigherOrder constructs as a wrapper type that takes a Meshes.jl geometry and provides an iterator of nodes.

Do you have suggestions on what would be a good API for consuming these nodes that are not part of the geometry per se, but are used for FEM and polynomial approximation?

jman87 · February 3, 2023, 3:55pm

There is no global option to set the VTK file version.

That’s all I was asking. And that’s fine. For my scripts I’ll just make a global version parameter that I can change and use in subsequent vtk_grid calls. Thanks for clarifying.

jman87 · February 3, 2023, 4:41pm

Hi @juliohm,

I’m probably not the best person to answer that question, but let me offer some rambling thoughts that are hopefully at least moderately interesting, and might even be helpful! It might also be good for you to open another topic to pull others into the discussion and explore this in more depth.

First, let me start with a little context. My background is primarily computational solid mechanics, with most of that being in explicit dynamics. My interest in Julia leans more toward ways to use it for “production work” and less toward research, though there’s obviously a lot of overlap. I definitely do not consider myself a Julia expert. My code tends to look more like Fortran95 written in Julia, though I’ve started incorporating a bit more of the multiple dispatch and custom struct features here lately. In recent years, we’ve made some good strides at developing a suite of second-order elements that work extremely well for explicit solid dynamics. Here’s a link to the latest development, recently published in CMAME: “Second-order pyramid element formulations suitable for lumped-mass explicit methods in nonlinear solid mechanics”

The full suite of elements, consisting of a HEX27, TET15, WEDGE21, and PYR19, are implemented in several government FE solvers, and are also supported in CUBIT for meshing, and mostly supported in ParaView and VTK (I say mostly because the pyramid is often missing the centroid node but that’s easily added). I have not, however, found much in the open source community supporting these element types, either in terms of meshing or FE solvers for explicit.

So one of my hobbies in my spare time has been putting together a block mesher similar to the approach taken by the government code InGrid, which was later incorporated into LS-DYNA’s LS-PrePost and also the commercial meshing software, TrueGrid. This block meshing approach is extremely powerful and can be very effective for production work. And the Julia interface seems to be the perfect place for such an open-source block mesher to exist.

In terms of file formats, it would be nice to support all the usual suspects, but my focus has been VTK since it’s easy to work with and the WriteVTK package has helped me speed this along tremendously. Longterm though, I think the ExodusII format is definitely needed, and the other usual suspects could be included (LS-DYNA, Abaqus, Ansys, etc.).

I’d be happy to discuss more if you’re interested. Hope this was helpful in some way.

juliohm · February 3, 2023, 5:37pm

Thank you for your reply @jman87. I am mostly interested in understanding how you consume these nodal information in algorithms. I am considering two APIs:

Given a Meshes.Geometry, we can create wrapper types Quadratic, Bilinear, … that take the geometry and provide additional nodes.
Alternatively, and my currently preferred API, we could simply introduce a new function nodes(geometry, :Quadratic) that defaults to the geometric nodes (i.e. vertices) of the geometry. Currently vertices(hexa) returns the 8 vertices needed to define the geometry, but we could add virtual nodes into this new nodes function.

juliohm · February 3, 2023, 5:40pm

Here is a MWE so that you have an idea of what I mean:

julia> using Meshes

julia> hexa = CartesianGrid(10,10,10)[1]
Hexahedron{3,Float64}
  └─Point(0.0, 0.0, 0.0)
  └─Point(1.0, 0.0, 0.0)
  └─Point(1.0, 1.0, 0.0)
  └─Point(0.0, 1.0, 0.0)
  └─Point(0.0, 0.0, 1.0)
  └─Point(1.0, 0.0, 1.0)
  └─Point(1.0, 1.0, 1.0)
  └─Point(0.0, 1.0, 1.0)

julia> vertices(hexa)
8-element Vector{Point3}:
 Point(0.0, 0.0, 0.0)
 Point(1.0, 0.0, 0.0)
 Point(1.0, 1.0, 0.0)
 Point(0.0, 1.0, 0.0)
 Point(0.0, 0.0, 1.0)
 Point(1.0, 0.0, 1.0)
 Point(1.0, 1.0, 1.0)
 Point(0.0, 1.0, 1.0)

The proposal would consist of adding a new function nodes that defaults to vertices. Whenever users provide a second argument like :Quadratic it does something else.

jman87 · February 3, 2023, 5:52pm

Sure, that makes sense. And I do probably need to look at Meshes.jl a bit more, as it may be I can reduce some of the overhead for my block mesher by leveraging Meshes.jl.

The way I’ve been envisioning it (possibly wrongly) is that the quadratic label would be associated with the element type and the node data would simply be a list of coordinates. I need to ponder this a bit more for my application.

But yes, I think you’ve got a good idea. The one important thing is that the higher order nodes need to be tied to the geometry. Some meshes I’ve seen can enhance first-order geometry with higher-order nodes, but the new nodes are simply linearly interpolated between the first-order nodes - thus, the curvature is still only modeled at a first-order resolution. As long as your new nodes are accounting for any geometrical features I think you should be fine.

juliohm · February 3, 2023, 6:16pm

Yes, I think it is just a historical artifact the fact that geometries are defined with variations like QuadraticHexa, etc. As far as I understand it, these are normal Hexahedron with additional virtual nodes, so no need for extra types. Also, consider contributing some of these ideas to Meshes.jl if you are planning to add geometric functionality instead of developing a separate package. We are adding features all the time over there.