VTK Lagrange elements in 3D elements is not rendered correctly

[jorgensd]

MWE using dolfinx:

from mpi4py import MPI
import vedo
import numpy as np
import dolfinx
import dolfinx.io

mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 1, 1, 1, dolfinx.cpp.mesh.CellType.hexahedron)
geo = mesh.geometry.x
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
topo = dolfinx.cpp.mesh.entities_to_geometry(mesh, mesh.topology.dim, np.arange(num_cells, dtype=np.int32), False)
perm_vtk = dolfinx.cpp.io.perm_vtk(mesh.topology.cell_type, topo.shape[1])
dolfin_to_vtk = np.zeros(topo.shape[1], dtype=np.int32)
for i in range(topo.shape[1]):
    dolfin_to_vtk[perm_vtk[i]] = i

dolfinx.io.VTKFile("mesh.pvd").write(mesh)
topo = topo[:, dolfin_to_vtk]
mesh_vedo = vedo.Mesh([geo, topo])
p = vedo.Plotter(shape=(1, 1), N=1, pos=(0, 0))
p.add(mesh_vedo)
p.show()
p.screenshot("mesh.png")
print(topo,"\n", geo)

yields:

[[1 4 6 2 0 5 7 3]]
 [[1. 0. 0.]
 [0. 0. 0.]
 [0. 1. 0.]
 [1. 1. 0.]
 [0. 0. 1.]
 [1. 0. 1.]
 [0. 1. 1.]
 [1. 1. 1.]]

and

while the correponding VTU file contains;

<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid"  version="0.1" >
<UnstructuredGrid>
<Piece  NumberOfPoints="8" NumberOfCells="1">
<Points>
<DataArray  type="Float64"  NumberOfComponents="3"  format="ascii">1 0 0  0 0 0  0 1 0  1 1 0  0 0 1  1 0 1  0 1 1  1 1 1  </DataArray>
</Points>
<Cells>
<DataArray  type="Int32"  Name="connectivity"  format="ascii">1 4 6 2 0 5 7 3  </DataArray>
<DataArray  type="Int32"  Name="offsets"  format="ascii">8 </DataArray>
<DataArray  type="Int8"  Name="types"  format="ascii">72 </DataArray>
</Cells>
</Piece>
</UnstructuredGrid>
</VTKFile>

Which is the same as the vedo mesh has, but yields

[jorgensd]

I actually get a similar issue with tetrahedron, so is there something Im missing?

[marcomusy]

I tried to install dolfinx but the instructions in https://github.com/FEniCS/dolfinx do not work.
In my ubuntu 20.04.1 I had to
sudo apt install libboost-timer-dev libboost-filesystem-dev

but then I get:

-- Test PETSC_TEST_RUNS with shared library linking - Success
-- HDF5: Using hdf5 compiler wrapper to determine C configuration
CMake Error at CMakeLists.txt:155 (message):
  Found serial HDF5 build, MPI HDF5 build required, try setting HDF5_DIR


-- Configuring incomplete, errors occurred!

trying:
sudo apt-get install libhdf5-mpich-dev
did not fix it.

Nonetheless I could reproduce the issue by copying the output you provided.
I'll get back to you as soon as i find out what's wrong..

[jorgensd]

I mostly use docker so all dependencies are right, docker pull dolfinx/dolfinx based on https://github.com/FEniCS/dolfinx/blob/master/docker/Dockerfile
The relevant part for you is:

apt-get install libmpich-dev libhdf5-mpich-dev

All dependencies are in the Dockerfile. If you rely on h5py, you should add:

export export HDF5_MPI="ON" && \
    export HDF5_DIR="/usr/lib/x86_64-linux-gnu/hdf5/mpich/" && \
    export CC=mpicc && \ 
    pip3 install --no-cache-dir --no-binary=h5py h5py

[marcomusy]

Thanks for the dockerfile instructions I'll give it a try later!

There are basically two separate issues with the original question.

1. The vedo.Mesh class is meant to only manage polygonal surface meshes so the code you have is merely creating a polygon of 8 edges. The class to be used should instead be UGrid:

import vedo

ug = vedo.UGrid('b.vtu').c('green').alpha(0.5)
vpts = vedo.Points(ug.points(), r=10).c('red')

vedo.show(ug,
          vpts,
          vpts.labels('id'),
          axes=1,
)

b.zip

2. Note that the vtu file you posted (even though paraview seems to be able to automagically fix it?) seems vtk-wise incorrect in the cell orientation (order of point insertion):
   

At the minute UGrid cannot be initialized by a sequence of points and connectivity (not implemented - but easy to add).
I may also add a HexMesh class analogous to the existing TetMesh.

[jorgensd]

I agree that the ordering seems strange, I think there is a bug in dolfinx after changing from FIAT to basix. I will get back to you with something with the correct ordering. It would be great if it was possible to extend UGrid to take in the topology, geometry, and celltype and return a mesh using the VTK arbitrary ordered lagrange elements: https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/
Note that there is a bug for second order hexahedral elements, see: https://gitlab.kitware.com/vtk/vtk/-/issues/17746

[marcomusy]

It would be great if it was possible to extend UGrid to take in the topology, geometry, and celltype and return a mesh using the VTK arbitrary ordered lagrange elements

by reading at the links I think it should be doable!

[jorgensd]

@marcomusy Btw. The bug with basix only affects second order geometries.
If we consider the data from vtu above, we can visualize the real ordering (when taking the ordering of the topology into account:

import vedo
import numpy as np
geo = np.array([[1., 0., 0.],
                [0., 0., 0.],
                [0., 1., 0.],
                [1., 1., 0.],
                [0., 0., 1.],
                [1., 0., 1.],
                [0., 1., 1.],
                [1., 1., 1.]])
topo = np.array([1, 4, 6, 2, 0, 5, 7, 3])
geo2 = np.zeros((len(topo), 3))
for i in range(len(topo)):
    geo2[i] = geo[topo[i], :]

pts = vedo.Points(geo2)
vedo.show(pts, pts.labels("id"))
vedo.screenshot("mesh.png")

yielding

which is the counter clockwise order (I admit that the choice of axis is weird, we will have a fix in dolfinx for this tomorrow).

[marcomusy]

..maybe i'm getting confused here but if I rotate by 90deg (which cannot change handedness) it actually looks clockwise - or anyway opposite to the green cube above:

import vedo
import numpy as np
geo = np.array([[1., 0., 0.],
                [0., 0., 0.],
                [0., 1., 0.],
                [1., 1., 0.],
                [0., 0., 1.],
                [1., 0., 1.],
                [0., 1., 1.],
                [1., 1., 1.]])
topo = np.array([1, 4, 6, 2, 0, 5, 7, 3])
geo2 = np.zeros((len(topo), 3))
for i in range(len(topo)):
    geo2[i] = geo[topo[i], :]

pts = vedo.Points(geo2).rotateY(-90)
vedo.show(pts, pts.labels("id"), axes=1)

[jorgensd]

If its clockwise or anti clockwise doesn't really matter, as one is just the reflection of the other, and should be able to render nicely (as it does in Paraview.:)

[jorgensd]

So the plan would be the following:

• Add support in UGrid for topology, geometry and celltype (integer identifier from VTK https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html?)
• Support arbitrary order lagrangian elements.

[marcomusy]

If its clockwise or anti clockwise doesn't really matter, as one is just the reflection of the other, and should be able to render nicely (as it does in Paraview.:)

True. What matters for the VTK readers is likely to be the order of insertion of the points which happens before the connectivity is loaded, and it's used to define the orientation. My guess is that paraview does it the other way around... but i'm not sure (i asked for it in the vtk discourse).

So the plan would be the following

If I understand, at the end of the day you will need something that generates polygons (and associated data) for visualization of the different elements, as all the analysis would be already made upstream.

[jorgensd]

Im not very familiar with VTK and vedo,but I figured that it would be something along the lines described in: https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/ to UGrid.
Please correct me if Im wrong.

[marcomusy]

Yes.
although I'm not super familiar with FEM :) my understanding is that you need to visualize your solutions as they are defined on the vertices and/or cells of the space discretization you choose. If that's the case extending the UGrid class as you say seems the right way to go.
Eventually what you want - i guess - is to be able to plot your mesh + solution with either ray-casting or by polygons (either surface-only or filled/shriked) ... maybe being able to slice through it etc.. E.g.

there are a few examples in https://vedo.embl.es/ (volumetric)

[jorgensd]

There are a few ways to think about this. In the simplest case, the mesh data might only be on vertices.
But lets say that we for instance use a mesh consisting of second order lagrangian elements, we would be able to define point data corresponding to the midpoint node of the facets as well (in the case of tetrahedron). For hexahedron, an additional mesh point will be added inside the cell. In any case, Currently, I would say we can assume a one to one correspondence between the nodes in the mesh and data-points (degrees of freedom) we would like to visualize. With a bit of hacking by duplication of grid nodes, it is quite straightforward to add discontinuous fields to such a visualization strategy (shown with pyvista in FEniCS/dolfinx#1161

The idea is to deprecate the matplotlib module in dolfin-x and rely on vedo as a backend for visualization.

[marcomusy]

Sounds good!
So what is the best way to proceed? Maybe you can workout on the dolfin-x side some example that generates different types of meshes and solutions you wish to visualize? Already in vtk format or just as numpy arrays?
Probably it's better if it's already in a vtk format to a) be more general, less dependent on vedo, b) have more control over it (and i guess you already have some working code)

[jorgensd]

I havent had a lot of time to look at this yet. However, if you want to start working on this, you could have a look at:
FEniCS/dolfinx#1161
where I extract data to use it with pyvista. As you can see, there are a few use-cases we would like to cover:

• Creating an arbitrary order Lagrangian mesh using geometry and topology (and vtk celltype) as input.
• Being able to attach point data and vector data to the geometry.

[marcomusy]

Hi @jorgensd , I'm trying the following:

sudo docker pull dolfinx/dolfinx
sudo docker run -ti dolfinx/dolfinx
pip3 install vedo
apt update
apt install nano

git clone https://github.com/FEniCS/dolfinx.git

cd dolfinx
git pull origin pull/1161/head

cd cpp
mkdir build
cd build
cmake ..
make install -j3
cd ../../python
pip3 install .

and it installs fine. But when import dolfinx i get:

> import dolfinx
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.8/dist-packages/dolfinx/__init__.py", line 31, in <module>
    from .cpp import __version__
ImportError: /usr/local/lib/python3.8/dist-packages/dolfinx/cpp.cpython-38-x86_64-linux-gnu.so: undefined symbol: _ZN7dolfinx2io5cells17get_vtk_cell_typeERKNS_4mesh4MeshEi

error msg seems has to do with method get_vtk_cell_type
I'm I doing something wrong?

[jorgensd]

Its tricky to reinstall dolfinx inside the prebuilt image. I would suggest using the dolfinx/dev-env image and use the following commands:

pip3  install vedo
# Install python components
pip3 install git+https://github.com/FEniCS/basix.git --upgrade && \
	pip3 install git+https://github.com/FEniCS/ufl.git --upgrade && \
	pip3 install git+https://github.com/FEniCS/ffcx.git --upgrade && \
	rm -rf /usr/local/include/dolfin /usr/local/include/dolfin.h

# Build C++ layer
git clone https://github.com/FEniCS/dolfinx.git && \
	 cd dolfinx/ && \
        git checkout  dokken/pyvista-demo && \
	 mkdir -p build && \
	 cd build && \
	 cmake -G Ninja -DCMAKE_BUILD_TYPE=Relase ../cpp/ && \
	 ninja -j3 install

cd dolfinx/python && \
pip3 -v install .

[marcomusy]

It still complains about MPI:

> cmake -G Ninja -DCMAKE_BUILD_TYPE=Relase ../
-- Could NOT find MPI_CXX (missing: MPI_CXX_WORKS) (Required is at least version "3")
CMake Error at /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:146 (message):
  Could NOT find MPI (missing: MPI_CXX_FOUND) (found suitable version "3.1",
  minimum required is "3")
Call Stack (most recent call first):
  /usr/share/cmake-3.16/Modules/FindPackageHandleStandardArgs.cmake:393 (_FPHSA_FAILURE_MESSAGE)
  /usr/share/cmake-3.16/Modules/FindMPI.cmake:1688 (find_package_handle_standard_args)
  CMakeLists.txt:89 (find_package)


-- Configuring incomplete, errors occurred!

but I've checked that I have both apt install libmpich-dev libhdf5-mpich-dev
It also complained about scotch, but this was fixed with apt install libscotch-dev

[jorgensd]

Let me have a look, and I'll give you a complete set of tested instructions

[jorgensd]

I forgot to add the PETSC_ARCH, i.e.

export  PETSC_ARCH="linux-gnu-real-32"
pip3  install vedo
# Install python components
pip3 install git+https://github.com/FEniCS/basix.git --upgrade && \
	pip3 install git+https://github.com/FEniCS/ufl.git --upgrade && \
	pip3 install git+https://github.com/FEniCS/ffcx.git --upgrade && \
	rm -rf /usr/local/include/dolfin /usr/local/include/dolfin.h

# Build C++ layer
git clone https://github.com/FEniCS/dolfinx.git && \
	 cd dolfinx/ && \
        git checkout  dokken/pyvista-demo && \
	 mkdir -p build && \
	 cd build && \
	 cmake -G Ninja -DCMAKE_BUILD_TYPE=Relase ../cpp/ && \
	 ninja -j3 install
source /usr/local/share/dolfinx/dolfinx.conf

cd ../python && \
pip3 -v install .

[marcomusy]

Thanks Jorgen! I'll get back to you as soon as i have something working smoothly.

[jorgensd]

@marcomusy Im looking forward to it! Tools such as vedo would make my life so much easier (after spending one weekend trying to generalize our matplotlib support to quad/hexes, and failing massively).

[marcomusy]

..it starts to take shape, although i'll need a couple of days more to finalize it..

all the element are created with an interface that looks like:

ug = UGrid([pts, cells, cellstypes])
ug.alpha(0.5).lineWidth(2).lighting('off')
show(ug)

[marcomusy]

.. while i'm at this - for my research - I'm also implementing a mesher, to mesh arbitrary concave domains, which seems to work pretty well and fast (~50k triangles/sec) with simple syntax:

from vedo import Spline, show
from vedo.pyplot import histogram

shape = Spline([[0, 0],
                [1, 0],
                [1.1, 4],
                [1.0, 1.5],
                [0.2, 5],
                [-1.0, 3.0],
                [0.4, 2.7],
                [-1.0, 2.4],
               ],
               closed=True).c('red4')

msh = shape.tomesh(resLine=400)  # resample boundary

msh.smoothLaplacian().addQuality().addScalarBar3D()
histo = histogram(msh.getCellArray('Quality'),
                  aspect=3/4, c='RdYlBu', xtitle='triangle mesh quality')
show(shape.ps(3), msh, histo, N=3, sharecam=False)

maybe dolfinx people will find it useful too

[jorgensd]

That is quite interesting! Are you planning to generalize it to other cell types (quads) or 3D?

[marcomusy]

oh i haven't thought about it! Not sure about quads, but for 3d it might be possible..

[jorgensd]

oh i haven't thought about it! Not sure about quads, but for 3d it might be possible..

Let me know how you get along with things:) i've spent the day updating my dolfin-x tutorial after we removed matplotlib. When vedo is ready, Im ready to add it to some of the examples, such as: https://jorgensd.github.io/dolfinx-tutorial/chapter3/em.html or https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code1.html
Could also add a meshing demo with vedo if you have a good idea for an example application.

[marcomusy]

Indeed it can be extended to quads (if this was what you meant):

with the bonus of allowing for variable resolution along x and/or y axes. In this case it's slightly slower (~30k quads/sec).

I also extended it to 3D convex shape to generate linear tets. I'll check if it's doable with concave shapes.

About the rest i'll give it a try today and tomorrow - i'll keep you informed :)

[marcomusy]

@jorgensd ..first checks with tet meshes..

import dolfinx
import dolfinx.io
import dolfinx.plot
import numpy as np
from mpi4py import MPI


try:
    import vedo
except ModuleNotFoundError:
    print("vedo is required for this demo")
    exit(0)


mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 4, 3, 5, cell_type=dolfinx.cpp.mesh.CellType.tetrahedron)

# Extract mesh data from dolfin-X
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
cell_entities = np.arange(num_cells, dtype=np.int32)
vtk_cells, cell_types = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim, cell_entities)

scals = np.linspace(0, 1, num=len(cell_types)).astype(float)
m = vedo.TetMesh([mesh.geometry.x, vtk_cells]).addCellArray(scals,'myscal').c('jet')

# m = vedo.UGrid([mesh.geometry.x, vtk_cells, cell_types]).tomesh(shrink=0.9)
# m.cutWithPlane(normal=[1,1,-0.5]).bc('tomato')

vedo.show(m, axes=1, offscreen=1, elevation=-50).screenshot()

I wonder if i should keep a specialized class for tetrahedral meshes or "unify" it as UGrid...

[jorgensd]

I would prefer if we have a general UGrid case (because we can use triangular, quadrilateral, hexahedron and tetrahedron elements with dolfinx), and having a common interface would help usability.