How to pass quadratic quadrilateral mesh from `dolfinx.io.gmshio` to PyVista?

General
1

Hello!

Brief description (goal)

I am trying to pass a mesh with qudrangles to PyVista. The mesh has lagrangean quadratic elements (“Lagrange”, degree=2), and is made with Gmsh. Then, it is imported with dolfinx.io.gmshio.

What is going on

I only see triangular elements on PyVista. The same thing happens with dolfinx.mesh.create_unit_square.

Minimal working example (MWE)

Both of these show the same with all the data formats.

::: captioned-content
::: caption
Here, I am only showing the code with dolfinx.mesh.create_unit_square.
:::

import basix
from dolfinx.fem import functionspace
from mpi4py import MPI
from dolfinx import fem
from dolfinx import mesh
from dolfinx import io
import numpy as np
import dolfinx.fem.petsc as petsc
from petsc4py.PETSc import ScalarType
from dolfinx import plot
import pyvista
import ufl


fe_degree = 2

domain = mesh.create_unit_square(
      MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
gdim = domain.geometry.dim

element = basix.ufl.element(
    family="Lagrange", cell=basix.CellType.quadrilateral,
    # family="Lagrange", cell=domain.basix_cell(),
    shape=(gdim,), degree=fe_degree)
V = functionspace(domain, element)

topology, cell_types, geometry = plot.vtk_mesh(V)
vtu = pyvista.UnstructuredGrid(
  topology, cell_types, geometry)
pvpl = pyvista.Plotter(window_size=[200, 200])
pvpl.add_mesh(vtu, show_edges=True)
pvpl.view_xy()
pvpl.show()

:::

import basix
from dolfinx.fem import functionspace
from mpi4py import MPI
from dolfinx import fem
import gmsh
from dolfinx import io
import numpy as np
import dolfinx.fem.petsc as petsc
from petsc4py.PETSc import ScalarType
from dolfinx import plot
import pyvista
import ufl


fe_degree = 2
simu_name = "quadrangles"
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
    gmsh.initialize()
    gmsh.option.set_number("Geometry.CopyMeshingMethod", 2)
    gmsh.option.setNumber('General.Verbosity', 0)
    gmsh.model.mesh.setOrder(fe_degree)
    # Do your best at creating quadrangles
    gmsh.option.setNumber("Mesh.Algorithm", 8)
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    gmsh.option.setNumber("Mesh.RecombineAll", 1)

    gmsh.model.add(simu_name)
    gmsh.model.setCurrent(simu_name)
    occ = gmsh.model.occ
    lc = 1e-3
    quad_mesh = True
    quart_orig = (0, 0, 0)
    square_ratio = 0.4
    n_pts = (3, 3)
    div_factor = (0.5, 0.5)
    mesh_type = ["Bump", "Bump"]
    tag = gmsh.model.occ.add_rectangle(0, 0, 0, 1, 1, roundedRadius=0)
    occ.synchronize()
    gmsh.model.add_physical_group(2, [tag], 1)
    occ.synchronize()
    gmsh.model.mesh.generate(2)

mesh_gmsh, cell_tags, facet_tags = io.gmshio.model_to_mesh(
    gmsh.model, mesh_comm, model_rank, gdim=2)
if mesh_comm.rank == model_rank:
    gmsh.fltk.run()
    gmsh.finalize()

domain = mesh_gmsh
gdim = domain.geometry.dim

element = basix.ufl.element(
    family="Lagrange", cell=basix.CellType.quadrilateral,
    # family="Lagrange", cell=domain.basix_cell(),
    shape=(gdim,), degree=fe_degree)
V = functionspace(domain, element)

topology, cell_types, geometry = plot.vtk_mesh(V)
vtu = pyvista.UnstructuredGrid(
  topology, cell_types, geometry)
pvpl = pyvista.Plotter(window_size=[200, 200])
pvpl.add_mesh(vtu, show_edges=True)
pvpl.view_xy()
pvpl.save_graphic("q.svg")
pvpl.show()

with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

Gmsh

gmsh

PyVista

PyVista

ParaView

paraview

XDMF (can’t upload otherwise)

::: captioned-content
::: caption
output.xdmf
:::

<?xml version="1.0"?>
<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>
<Xdmf Version="3.0" xmlns:xi="https://www.w3.org/2001/XInclude">
  <Domain>
    <Grid Name="mesh" GridType="Uniform">
      <Topology TopologyType="Quadrilateral" NumberOfElements="64" NodesPerElement="4">
        <DataItem Dimensions="64 4" NumberType="Int" Format="HDF">output.h5:/Mesh/mesh/topology</DataItem>
      </Topology>
      <Geometry GeometryType="XY">
        <DataItem Dimensions="81 2" Format="HDF">output.h5:/Mesh/mesh/geometry</DataItem>
      </Geometry>
    </Grid>
  </Domain>
</Xdmf>

:::

My system

(I also tested on 0.7.0.dev0~r27554~.cfeffe0-1)

- FEniCSx software
  dolfinx: 0.8.0.dev0_r27568.a34b0c9-1
  basix: 0.8.0.dev0_r952.eebbc52-1
  ufl: 2023.3.0.dev0_r3568.f132188-1
  ffcx: 0.8.0.dev0_r7083.b4cb217-1

- Dependencies
  python: 3.11.5-2
  python-numpy: 1.26.0-1
  petsc: 3.19.5.1171.g37df9106526-1
  hdf5-openmpi: 1.14.2-1
  boost: 1.83.0-2
  adios2: 2.8.3-5
  scotch: 7.0.4-1
  pybind11: 2.11.1-1
  python-build: 1.0.1-1
  python-cffi: 1.15.1-4
  python-cppimport: 22.08.02.r6.g0849d17-1
  python-installer: 0.7.0-3
  python-mpi4py: 3.1.4-3
  python-pytest: 7.4.2-1
  python-scikit-build: 0.17.6-1
  python-setuptools: 1:68.0.0-1
  python-wheel: 0.40.0-3
  xtensor: 0.24.0-2
  xtensor-blas: 0.20.0-1

- Operating system
  Arch Linux: 6.5.4-arch2-1
2

Pyvista has an issue with edge visualization of higher order, see: VTK_LAGRANGE_HEXAHEDRON not supported with showing edges · Issue #947 · pyvista/pyvista · GitHub
Thus the mesh visualized is a quadrilateral mesh, pyvista just adds multiple extra edges to it.

2 Likes