Plot mesh parallel, show partitioning

[Stein]

I’d like to pyvista-plot my mesh, when running on multiple processors, with a flag for showing mesh partitioning or not (the ‘or not’ is where I struggle).

I am following How to collect the global mesh without writing a file - #3 by edgar and Pyvista, how to gather data to plot on rank=0 - #3 by LiborKudela, but would like to avoid creating a FunctionSpace if not necessary. So I got the following MWE:

import pyvista,dolfinx

show_partitioning = True

domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,[(0,0),(1,1)], (10,10))
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(domain)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_cells_local_geom = domain.geometry.index_map().size_local

# Gather data
num_dofs_per_cell = topology[0]
global_topology = domain.comm.gather(
    topology[: (num_dofs_per_cell + 1) * num_cells_local], root=0
)
if show_partitioning:
    global_geometry = domain.comm.gather(geometry[:, :], root=0)
else:
    # FAULTY LINE IS HERE:
    global_geometry = domain.comm.gather(geometry[:num_cells_local_geom, :], root=0)
global_ct = domain.comm.gather(cell_types[:num_cells_local], root=0)

if domain.comm.rank == 0:
    root_geoms = [np.vstack(global_geometry)]  if not show_partitioning else global_geometry
    root_tops = [np.concatenate(global_topology)]  if not show_partitioning else global_topology
    root_cts = [np.concatenate(global_ct)]  if not show_partitioning else global_ct
    
    plotter = pyvista.Plotter()
    colors = ['k','r','b','g']
    for color, root_top, root_ct, root_geom in zip(colors,root_tops, root_cts, root_geoms):
        grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
        plotter.add_mesh(grid, style="wireframe", color=color, line_width=2)
    plotter.view_xy()
    plotter.show()

This works nicely for --np 1 (obviously) and works nicely with --np 2 and show_partitioning = True, but shows something a wrong with --np 2 and show_partitioning = False.

I know the error must be in one particular line (highlighted in the MWE), but can’t figure out what it should be. A hint is much appreciated!

[dokken]

I am not quite sure what you are trying to achieve here.
You have ignored the crucial:

# Map to global dof indices
global_dofs = V.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
# Overwrite topology
topology[topology_dofs] = global_dofs

How to collect the global mesh without writing a file

Here is an example code of how to do it: # Copyright (C) 2023 Jørgen S. Dokken # # SPDX-License-Identifier: MIT # # Code to plot mesh with pyvista in parallel by gather data on a root process # # Last changed: 2023-10-23 import dolfinx from mpi4py import MPI import pyvista import numpy as np pyvista.start_xvfb(0.2) mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2)) u = dolfinx.fem.Function(V) u.interpolate(lambda x: x[0] + 2*…

This is crucial as one maps the local topology indices to a single global numbering.
For more about this see for instance section 6.2 of:

DOLFINx: The next generation FEniCS problem solving environment

which tries to explain global ordering.

[Stein]

That code snipped doesn’t quite work for me since I want to operate directly on the mesh, not on a functionspace:

I’m trying to adopt

topology, cell_types, geometry = dolfinx.plot.vtk_mesh(domain)

instead of

topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)

[dokken]

The geometry also has an index map, I.e mesh.geometry.index_map.local_to_global that you have to use to map the topology array retrieved by plot.vtk_mesh to its global indices before gather them.

[Stein]

I believe I tried that before, as I figured the geometry object would behave almost identical to a Function in a Functionspace. But I’m getting the error message:

global_dofs = domain.geometry.index_map.local_to_global(topology[topology_dofs].copy())
AttributeError: 'nanobind.nb_method' object has no attribute 'local_to_global'

I’m running v0.8.0. Should I upgrade to 0.10.0? Comparing dolfinx.cpp.mesh — DOLFINx 0.8.0 documentation to dolfinx.mesh — DOLFINx 0.10.0.0 documentation it looks like the interface changed (but not from v0.8.0 to v0.9.0 though?).

In your view, what are the primary resources for developing ones understanding of Fenics’ inner workings? I notice I’m often limited to manipulating snippets that I find online. I also notice that FenicsX asks its users to have a deeper understanding of the software than Fenics. The paper you cite is a good start, I had seen it a while back but I’ll give it a second read.

[dokken]

@Stein

global_dofs = domain.geometry.index_map.local_to_global(topology[topology_dofs].copy())

Try
global_dofs = domain.geometry.index_map().local_to_global(topology[topology_dofs].copy())

@Stein

In your view, what are the primary resources for developing ones understanding of Fenics’ inner workings? I notice I’m often limited to manipulating snippets that I find online. I also notice that FenicsX asks its users to have a deeper understanding of the software than Fenics. The paper you cite is a good start, I had seen it a while back but I’ll give it a second read.

The paper is a good start.
I would also consider
http://jsdokken.com/FEniCS-workshop/README.html
and
https://scientificcomputing.github.io/mpi-tutorial/notebooks/dolfinx_MPI_tutorial.html
that goes through some notions.
Also the short note I made on mesh creation might be of help
http://jsdokken.com/dolfinx_docs/meshes.html

[Stein]

Ah damn, sorry, that was a novice mistake. That seems to have done it, thanks.

Current working MWE:

import pyvista,dolfinx

show_partitioning = False # or True

domain = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD,[(0,0),(1,1)], (10,10))
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(domain)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_cells_local_geom = domain.geometry.index_map().size_local
num_dofs_per_cell = topology[0]

# Gather data
if show_partitioning:
    global_geometry = domain.comm.gather(geometry[:, :], root=0)
else:
    global_geometry = domain.comm.gather(geometry[:num_cells_local_geom, :], root=0)

    # Map topology to global dof indices
    topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell+1)) != 0
    global_dofs = domain.geometry.index_map().local_to_global(topology[topology_dofs].copy())
    topology[topology_dofs] = global_dofs

global_topology = domain.comm.gather(topology[: (num_dofs_per_cell + 1) * num_cells_local], root=0)
global_ct = domain.comm.gather(cell_types[:num_cells_local], root=0)

if domain.comm.rank == 0:
    root_geoms = [np.vstack(global_geometry)]  if not show_partitioning else global_geometry
    root_tops = [np.concatenate(global_topology)]  if not show_partitioning else global_topology
    root_cts = [np.concatenate(global_ct)]  if not show_partitioning else global_ct
    
    plotter = pyvista.Plotter()
    colors = ['k','r','b','m']
    for color, root_top, root_ct, root_geom in zip(colors,root_tops, root_cts, root_geoms):
        grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
        plotter.add_mesh(grid, style="wireframe", color=color, line_width=2)
    plotter.view_xy()
    plotter.show()

And thanks for the addititional resources. I’ll take a look.

It seems like a lot of great resources are available (scifem, oasix, comet, ambit, donfin_dg, etc), but they are not too easily found if one doesn’t know they exist. Maybe it’s an option to create a separate “Category” on the Discourse where people can promote their work?

[dokken]

I’ll make a category for that, as it lowers the threshold for promoting. We have had internal discussions on how to promote external libraries.