How to collect the global mesh without writing a file

edgar · October 12, 2023, 5:24pm

Hello. If one runs the tutorial about subdomains in parallel (c9214c68438ab27d2fa8936307027c5209272162):

mpirun -np 2 python3 subdomains.py

The plots of the meshes are split (among processes)

Is there a way to do some sort of collect with MPI to show the whole mesh? I’m looking for a solution which works with Gmsh’s API for Python. Thanks!

For reference, the relevant (slightly modified) code of subdomains.py:

import gmsh
import numpy as np
import pyvista
from dolfinx.fem import (
    dirichletbc, Function, FunctionSpace, locate_dofs_geometrical,)
from dolfinx.mesh import create_unit_square, locate_entities
from dolfinx.plot import vtk_mesh
from dolfinx import io
from ufl import (TestFunction, TrialFunction)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType

bottom_marker = 1
top_marker = 2
left_marker = 1

gmsh.initialize()
proc = MPI.COMM_WORLD.rank
model = gmsh.model.occ
if proc == 0:
    # Create one rectangle for each subdomain
    model.add_rectangle(0, 0, 0, 1, 0.5, tag=1)
    model.add_rectangle(0, 0.5, 0, 1, 0.5, tag=2)

    # Fuse the two rectangles and keep their interface
    model.fragment([(2, 1)], [(2, 2)])
    model.synchronize()

    # Mark the top (2) and bottom (1) rectangle (2D)
    top, bottom = None, None
    for surf_dimtag in gmsh.model.get_entities(dim=2):
        com = model.get_center_of_mass(*surf_dimtag)
        if np.allclose(com, [0.5, 0.25, 0]):
            bottom = surf_dimtag[1]
        else:
            top = surf_dimtag[1]
    gmsh.model.add_physical_group(2, [bottom], bottom_marker)
    gmsh.model.add_physical_group(2, [top], top_marker)

    # Tag the left boundary (1D)
    left = []
    for line in gmsh.model.get_entities(dim=1):
        com = model.get_center_of_mass(line[0], line[1])
        if np.isclose(com[0], 0):
            left.append(line[1])
    gmsh.model.add_physical_group(1, left, left_marker)

    gmsh.model.mesh.generate(2)

mesh_gmsh, ct, ft = io.gmshio.model_to_mesh(
    gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
cell_tags, facet_tags = ct, ft
gmsh.finalize()

domain = mesh_gmsh


Q = FunctionSpace(domain, ("DG", 0))


def Omega_0(x):
    return x[1] <= 0.5

def Omega_1(x):
    return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = locate_entities(domain, domain.topology.dim, Omega_1)

kappa.x.array[cells_0] = np.full_like(cells_0, 1, dtype=ScalarType)
kappa.x.array[cells_1] = np.full_like(cells_1, 0.1, dtype=ScalarType)

V = FunctionSpace(domain, ("CG", 1))
u, v = TrialFunction(V), TestFunction(V)
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
bcs = [dirichletbc(ScalarType(1), dofs, V)]

num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
marker = np.zeros(num_cells_local, dtype=np.int32)
cells_0 = cells_0[cells_0<num_cells_local]
cells_1 = cells_1[cells_1<num_cells_local]
marker[cells_0] = 1
marker[cells_1] = 2
topology, cell_types, x = vtk_mesh(domain, domain.topology.dim, np.arange(num_cells_local, dtype=np.int32))

p = pyvista.Plotter(window_size=[800, 800])
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
grid.cell_data["Marker"] = marker
grid.set_active_scalars("Marker")
p.add_mesh(grid, show_edges=True)
if not pyvista.OFF_SCREEN:
    p.show()

My system

FEniCSx software
dolfinx: 0.7.0.dev0_r27554.cfeffe0-1
basix: 0.7.0.dev0_r945.1117a8d-1
ufl: 2023.2.0.dev0_r3562.77ae57c-1
ffcx: 0.7.0.dev0_r7077.1d27238-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

dokken · October 12, 2023, 5:33pm

What in particular do you want to do with Gmsh with the mesh post using it in dolfinx?

edgar · October 12, 2023, 6:01pm

It’s mostly to automate visualisation with PyVista. At the moment, I don’t have anything meaningful, but would prefer to generate plots of the whole thing programmatically instead of saving a file and opening it on ParaView. The comment about Gmsh is because the alternative is to save a mesh file, which I want to avoid.

dokken · October 12, 2023, 8:25pm

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*x[1]**2)


# Create local VTK mesh data structures
topology, cell_types, geometry = dolfinx.plot.vtk_mesh(V)
num_cells_local = mesh.topology.index_map(mesh.topology.dim).size_local
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
# VTK topology stored as (num_dofs_cell_0, dof_0,.... dof_(num_dofs_cell_0, num_dofs_cell_1, ....))
# We assume we only have one cell type and one dofmap, thus every `num_dofs_per_cell` is the same
num_dofs_per_cell = topology[0]
# Get only dof indices
topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell+1)) != 0

# 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
# Choose root
root = 0

# Gather data
global_topology = mesh.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=root)
global_geometry = mesh.comm.gather(geometry[:V.dofmap.index_map.size_local,:], root=root)
global_ct = mesh.comm.gather(cell_types[:num_cells_local])
global_vals = mesh.comm.gather(u.x.array[:num_dofs_local])

if mesh.comm.rank == root:
    # Stack data
    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)
    root_vals = np.concatenate(global_vals)

    # Plot as in serial
    pyvista.OFF_SCREEN = True
    grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
    grid.point_data["u"] = root_vals
    grid.set_active_scalars("u")
    u_plotter = pyvista.Plotter()
    u_plotter.add_mesh(grid, show_edges=True)
    u_plotter.view_xy()
    u_plotter.screenshot("figure.png")

I’ll add it to the list of add to my tutorial when I have time for it.

edgar · October 12, 2023, 10:35pm

Thanks! This is great.

jormoral · December 12, 2023, 1:46pm

Dear dokken,

This solution works for me until I use dolfinx.fem.VectorFunctionSpace instead of dolfinx.fem.FunctionSpace.

How would you modify the example code you provided if V = dolfinx.fem.VectorFunctionSpace(mesh, (“Lagrange”,2))?

edgar · December 12, 2023, 6:33pm

I hope that this helps: How to use the shape parameter of dolfinx.fem.FunctionSpace?, but the discourse interface is telling me that

If you have an unrelated issue, please start a new topic instead.

Topic		Replies	Views
Write and read back mesh in parallel with same meshtags dolfinx mesh , mpi	12	139	March 25, 2024
Visualisation of distributed mesh mesh	7	859	August 24, 2022
Reproducibility in parallel programming dolfinx	1	357	August 8, 2022
Gathering function on root process (MPI) issue with mixed element General	2	87	February 20, 2024
Issues with saving mesh connectivity from parellel Legacy DOLFIN mesh	7	200	December 15, 2023

How to collect the global mesh without writing a file

My system

Related Topics