Unexpected behavior of SpatialCoordinate, when using create_submesh

mesh
1

Dear Community,

when creating a submesh using dolfinx.mesh.create_submesh and visualizing the spatial coordinates I find that the mesh is somehow destroyed/ the triangles appear to be wrongly associated.
The phenomenon that I see here also appears, when using

x = ufl.SpatialCoordinate(submesh)

I created a MWE for this:

import numpy as np
from mpi4py import MPI
import dolfinx
from dolfinx import fem
from dolfinx.mesh import locate_entities, meshtags

# 1. Create base mesh
domain = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 20, 20)
tdim = domain.topology.dim
# 2. Tag subdomain
def left_half(x):
    return x[0]<= 0.5 + 1e-14

cells = locate_entities(domain, tdim, left_half)
ct = meshtags(domain, tdim, cells, np.full(len(cells), 20, dtype=np.int32))

# 3. Create SubMesh
submesh, submesh_to_mesh, vertex_map = dolfinx.mesh.create_submesh(
    domain, tdim, ct.find(20))[0:3]

# 4. Function space on SubMesh
Vb = fem.functionspace(submesh, ("Lagrange", 1))

# 5. Plot the grid
import pyvista
topology, cell_types,_ = dolfinx.plot.vtk_mesh(Vb)
grid = pyvista.UnstructuredGrid(topology, cell_types, submesh.geometry.x)
xvals = submesh.geometry.x
grid.point_data["x0"] = xvals[:, 0]
plotter = pyvista.Plotter()
plotter.add_mesh(grid, scalars="x0", show_edges=True)
plotter.show()

See also the following plot:

coordinates
coordinates1024×768 130 KB

Interestingly, this does not appear, when changing the definition left_half to

def left_half(x):
    return x[1]<= 0.5 + 1e-14

I appreciate any help or insights on what might be the problem or what I am doing wrong here.
Best,
Marvin

2

This is just a visualisation issue. In your MWE you’re confusing a function space’s VTK representation with the mesh’s VTK representation.

The following creates a VTK mesh of the function space Vb so you’d expect to use the corresponding geometry points defined by that function space. So instead of

topology, cell_types,_ = dolfinx.plot.vtk_mesh(Vb)
grid = pyvista.UnstructuredGrid(topology, cell_types, submesh.geometry.x)
xvals = submesh.geometry.x
grid.point_data["x0"] = xvals[:, 0]

you’d opt for something like

topology, cell_types, Vb_x = dolfinx.plot.vtk_mesh(Vb)
grid = pyvista.UnstructuredGrid(topology, cell_types, Vb_x)
grid.point_data["x0"] = Vb_x[:, 0]

This is designed for visualising functions defined on arbitrary order function spaces such that the VTK mesh may be finer than the original FE mesh.

The reason the function space’s coordinates do not line up with the mesh geometry is due to a reordering based on the degrees of freedom in that function space.

If you wanted to purely plot the mesh itself you could do:

topology, cell_types, x = dolfinx.plot.vtk_mesh(submesh)
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.point_data["x0"] = submesh.geometry.x[:, 0]
assert np.all(np.isclose(x, submesh.geometry.x))  # Sanity check both geometries are the same.

By making this changes I see for the function space VTK plot and the mesh VTK plot, respectively:

image
image840×562 109 KB

image
image840×562 109 KB

I have some examples of plotting here.

3 Likes
3

@nate’s answer explain the faulty visualization that you observe. It is a little annoying, though, that the pyvista code that one has to write to debug ones code is itself bug sensitive :sweat_smile:

In that regard, if I may shamelessly self-promote a different visualization package (pip install pyvista4dolfinx), the following also works:

You can replace…

import pyvista
topology, cell_types,_ = dolfinx.plot.vtk_mesh(Vb)
grid = pyvista.UnstructuredGrid(topology, cell_types, submesh.geometry.x)
xvals = submesh.geometry.x
grid.point_data["x0"] = xvals[:, 0]
plotter = pyvista.Plotter()
plotter.add_mesh(grid, scalars="x0", show_edges=True)
plotter.show()

by…

from pyvista4dolfinx import plot
plot(Vb.mesh,show=True)
1 Like