Submesh from parallel mesh

[Chaitanya_Apte]

import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

import matplotlib.pyplot as plt
import dolfinx

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.model.add("square")

gmsh.model.occ.addRectangle(0,0,0,1,1,1)
gmsh.model.occ.addRectangle(1,0,0,1,1,2)
gmsh.model.occ.addRectangle(0,1,0,1,1,3)
gmsh.model.occ.addRectangle(1,1,0,1,1,4)

gmsh.model.occ.fragment([(2,1)],[(2,2),(2,3),(2,4)],removeObject=True,removeTool=True)

gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)

for i in lines:
    gmsh.model.mesh.setTransfiniteCurve(i[1],5)

for i in surf:
    gmsh.model.mesh.setTransfiniteSurface(i[1])
    gmsh.model.mesh.setRecombine(i[0],i[1])


gmsh.model.addPhysicalGroup(2,[1],101)
gmsh.model.setPhysicalName(2,101,"d1")

gmsh.model.addPhysicalGroup(2,[2],102)
gmsh.model.setPhysicalName(2,102,"d2")

gmsh.model.addPhysicalGroup(2,[3],103)
gmsh.model.setPhysicalName(2,103,"d3")

gmsh.model.addPhysicalGroup(2,[4],104)
gmsh.model.setPhysicalName(2,104,"d4")


gmsh.model.mesh.generate(2)

# if '-nopopup' not in sys.argv:
#     gmsh.fltk.run()
gmsh.write("mesh.msh")
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


####################################################################################################################


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)


# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = np.concatenate((markers.find(101),markers.find(102)))
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = np.concatenate((markers.find(103),markers.find(104)))
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


####################################################################################################################




V = fem.functionspace(domain, ("Lagrange", 1))

# Create local VTK mesh data structures
topology, cell_types, geometry = plot.vtk_mesh(V)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_dofs_local = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
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

# Gather data
global_topology = domain.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=0)
global_geometry = domain.comm.gather(geometry[:V.dofmap.index_map.size_local,:], root=0)
global_ct = domain.comm.gather(cell_types[:num_cells_local])


global_def = domain.comm.gather(uh.x.array[:num_dofs_local])


if rank == 0:
    # Stack data
    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)
        
    root_def = np.concatenate(global_def)

    
    domain,markers,facets = io.gmshio.read_from_msh("mesh.msh",MPI.COMM_SELF)
    combined_indices = []
    for tag in [101,104]:
       combined_indices.append(markers.find(tag))
    combined_indices = np.sort(np.hstack(combined_indices))

    submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)
    V = fem.functionspace(domain, ("Lagrange", 1))

    uh = fem.Function(V)
    uh.x.array[:] = root_def
    Vs = fem.functionspace(submesh, ("Lagrange", 1))
    uDs = fem.Function(Vs)
    uDs.interpolate(uh,None,entity_map)
    
    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(Vs)
    
    u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
    u_grid.point_data["u"] = uDs.x.array
    u_grid.set_active_scalars("u")
    u_plotter = pyvista.Plotter()
    u_plotter.add_mesh(u_grid, show_edges=True)
    u_plotter.view_xy()
    if not pyvista.OFF_SCREEN:
        u_plotter.show()

I want to visualize a subdomain. So i create a submesh and plot only the required subdomains from markers. I dont know how to reconstruct domain, markers, facets from parallelized mesh, so I am re-reading the whole mesh on root process.
But the above code gives different outputs in serial(as expected) and parallel.

[dokken]

See: Add `transfer_meshtags_to_submesh` · Issue #3096 · FEniCS/dolfinx · GitHub which can transfer meshtags from parent to submesh.

[Chaitanya_Apte]

@dokken

See: Add transfer_meshtags_to_submesh · Issue #3096 · FEniCS/dolfinx · GitHub which can transfer meshtags from parent to submesh.

Okay.

Also, can this step be eliminated? I dont want to save the mesh file and re-read it.

@Chaitanya_Apte

    domain,markers,facets = io.gmshio.read_from_msh("mesh.msh",MPI.COMM_SELF)

Can i assemble domain, markers and facets from all processes and use them as dolfinx.mesh objects just like in a serial code.

[dokken]

Not really sure what you are trying to do, please clarify.

In the previous post I showed how to transfer meshtags to the sub-mesh. Then you can output that for visualization/debugging.

If you want to use pyvista to plot in parallel, you should look at:

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*…

[Chaitanya_Apte]

To create submesh does the domain which goes as input to the function has to be the entire mesh?? :

@Chaitanya_Apte

submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)

Because when i try to run the below code on 4 processes (mpirun -n 4 python3 test.py), it keeps on running forever without any output.

import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

import matplotlib.pyplot as plt
import dolfinx

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.model.add("square")

gmsh.model.occ.addRectangle(0,0,0,1,1,1)
gmsh.model.occ.addRectangle(1,0,0,1,1,2)
gmsh.model.occ.addRectangle(0,1,0,1,1,3)
gmsh.model.occ.addRectangle(1,1,0,1,1,4)

gmsh.model.occ.fragment([(2,1)],[(2,2),(2,3),(2,4)],removeObject=True,removeTool=True)

gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)

for i in lines:
    gmsh.model.mesh.setTransfiniteCurve(i[1],5)

for i in surf:
    gmsh.model.mesh.setTransfiniteSurface(i[1])
    gmsh.model.mesh.setRecombine(i[0],i[1])


gmsh.model.addPhysicalGroup(2,[1],101)
gmsh.model.setPhysicalName(2,101,"d1")

gmsh.model.addPhysicalGroup(2,[2],102)
gmsh.model.setPhysicalName(2,102,"d2")

gmsh.model.addPhysicalGroup(2,[3],103)
gmsh.model.setPhysicalName(2,103,"d3")

gmsh.model.addPhysicalGroup(2,[4],104)
gmsh.model.setPhysicalName(2,104,"d4")


gmsh.model.mesh.generate(2)

# if '-nopopup' not in sys.argv:
#     gmsh.fltk.run()
#gmsh.write("mesh.msh")
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


####################################################################################################################


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)


# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = np.concatenate((markers.find(101),markers.find(102)))
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = np.concatenate((markers.find(103),markers.find(104)))
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


####################################################################################################################


combined_indices = []
for tag in [101,104]:
    combined_indices.append(markers.find(tag))
combined_indices = np.sort(np.hstack(combined_indices))

if len(combined_indices) != 0:
 
    submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)
    
    Vs = fem.functionspace(submesh, ("Lagrange", 1))
    uDs = fem.Function(Vs)
    uDs.interpolate(uh,None,entity_map)

    
    # Create local VTK mesh data structures
    topology, cell_types, geometry = plot.vtk_mesh(Vs)
    num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
    num_dofs_local = Vs.dofmap.index_map.size_local * Vs.dofmap.index_map_bs
    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 = Vs.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
    # Overwrite topology
    topology[topology_dofs] = global_dofs

    # Gather data
    global_topology = domain.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=0)
    global_geometry = domain.comm.gather(geometry[:Vs.dofmap.index_map.size_local,:], root=0)
    global_ct = domain.comm.gather(cell_types[:num_cells_local])

    global_def = domain.comm.gather(uDs.x.array[:num_dofs_local])
    
else:
    global_topology = None
    global_ct = None
    global_geometry = None
    global_def = None
    


if rank == 0:
    # Stack data
    
    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)
        
    root_def = np.concatenate(global_def)
            
    u_grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
    u_grid.point_data["u"] = root_def
 
    u_plotter = pyvista.Plotter()
    u_plotter.add_mesh(u_grid, show_edges=True)

    if not pyvista.OFF_SCREEN:
        u_plotter.show()

Works fine when running in serial (python3 test.py)

[dokken]

@Chaitanya_Apte

if len(combined_indices) != 0:
 
    submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)

This should not be within an if-test. You should call this even if combined_indices is an empty list.

[Chaitanya_Apte]

removing the if-test still gives the same issue.

Surprisingly, it works when parallelizing on 2 or 3 processes
but issue persists for 4 or more processes

here is the updated code:

import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

import matplotlib.pyplot as plt
import dolfinx

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.model.add("square")

gmsh.model.occ.addRectangle(0,0,0,1,1,1)
gmsh.model.occ.addRectangle(1,0,0,1,1,2)
gmsh.model.occ.addRectangle(0,1,0,1,1,3)
gmsh.model.occ.addRectangle(1,1,0,1,1,4)

gmsh.model.occ.fragment([(2,1)],[(2,2),(2,3),(2,4)],removeObject=True,removeTool=True)

gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)

for i in lines:
    gmsh.model.mesh.setTransfiniteCurve(i[1],20)

for i in surf:
    gmsh.model.mesh.setTransfiniteSurface(i[1])
    gmsh.model.mesh.setRecombine(i[0],i[1])


gmsh.model.addPhysicalGroup(2,[1],101)
gmsh.model.setPhysicalName(2,101,"d1")

gmsh.model.addPhysicalGroup(2,[2],102)
gmsh.model.setPhysicalName(2,102,"d2")

gmsh.model.addPhysicalGroup(2,[3],103)
gmsh.model.setPhysicalName(2,103,"d3")

gmsh.model.addPhysicalGroup(2,[4],104)
gmsh.model.setPhysicalName(2,104,"d4")


gmsh.model.mesh.generate(2)

# if '-nopopup' not in sys.argv:
#     gmsh.fltk.run()
#gmsh.write("mesh.msh")
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


####################################################################################################################


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)


# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = np.concatenate((markers.find(101),markers.find(102)))
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = np.concatenate((markers.find(103),markers.find(104)))
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


####################################################################################################################


combined_indices = []
for tag in [101,104]:
    combined_indices.append(markers.find(tag))
combined_indices = np.sort(np.hstack(combined_indices))



submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)

Vs = fem.functionspace(submesh, ("Lagrange", 1))
uDs = fem.Function(Vs)
uDs.interpolate(uh,None,entity_map)


# Create local VTK mesh data structures
topology, cell_types, geometry = plot.vtk_mesh(Vs)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_dofs_local = Vs.dofmap.index_map.size_local * Vs.dofmap.index_map_bs
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 = Vs.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
# Overwrite topology
topology[topology_dofs] = global_dofs

# Gather data
global_topology = domain.comm.gather(topology[:(num_dofs_per_cell+1)*num_cells_local], root=0)
global_geometry = domain.comm.gather(geometry[:Vs.dofmap.index_map.size_local,:], root=0)
global_ct = domain.comm.gather(cell_types[:num_cells_local])

global_def = domain.comm.gather(uDs.x.array[:num_dofs_local])
    


if rank == 0:
    # Stack data
    
    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)
        
    root_def = np.concatenate(global_def)
            
    u_grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
    u_grid.point_data["u"] = root_def
 
    u_plotter = pyvista.Plotter()
    u_plotter.add_mesh(u_grid, show_edges=True)

    if not pyvista.OFF_SCREEN:
        u_plotter.show()

[dokken]

The deadlock happens at the interpolation stage, and has been fixed in:
Make standalone function for interpolation over non-matching meshes by jorgensd · Pull Request #3177 · FEniCS/dolfinx · GitHub and can be found on the
ghcr.io/fenics/dolfinx/lab:nightly image.

The following code (shown below with some minor changes compared to yours runs fine on 4 processes:

import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

import matplotlib.pyplot as plt
import dolfinx

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 0)
gmsh.model.add("square")

gmsh.model.occ.addRectangle(0, 0, 0, 1, 1, 1)
gmsh.model.occ.addRectangle(1, 0, 0, 1, 1, 2)
gmsh.model.occ.addRectangle(0, 1, 0, 1, 1, 3)
gmsh.model.occ.addRectangle(1, 1, 0, 1, 1, 4)

gmsh.model.occ.fragment(
    [(2, 1)], [(2, 2), (2, 3), (2, 4)], removeObject=True, removeTool=True
)

gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)

for i in lines:
    gmsh.model.mesh.setTransfiniteCurve(i[1], 20)

for i in surf:
    gmsh.model.mesh.setTransfiniteSurface(i[1])
    gmsh.model.mesh.setRecombine(i[0], i[1])


gmsh.model.addPhysicalGroup(2, [1], 101)
gmsh.model.setPhysicalName(2, 101, "d1")

gmsh.model.addPhysicalGroup(2, [2], 102)
gmsh.model.setPhysicalName(2, 102, "d2")

gmsh.model.addPhysicalGroup(2, [3], 103)
gmsh.model.setPhysicalName(2, 103, "d3")

gmsh.model.addPhysicalGroup(2, [4], 104)
gmsh.model.setPhysicalName(2, 104, "d4")


gmsh.model.mesh.generate(2)

# if '-nopopup' not in sys.argv:
#     gmsh.fltk.run()
# gmsh.write("mesh.msh")
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


####################################################################################################################


V = fem.functionspace(domain, ("Lagrange", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)


# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)

Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = np.concatenate((markers.find(101), markers.find(102)))
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = np.concatenate((markers.find(103), markers.find(104)))
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

problem = LinearProblem(
    a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
uh = problem.solve()


####################################################################################################################


combined_indices = []
for tag in [101, 104]:
    combined_indices.append(markers.find(tag))
combined_indices = np.sort(np.hstack(combined_indices))


submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(
    domain, domain.topology.dim, combined_indices
)

Vs = fem.functionspace(submesh, ("Lagrange", 1))
uDs = fem.Function(Vs)
uDs.interpolate(
    uh, cells0=entity_map, cells1=np.arange(len(entity_map), dtype=np.int32)
)

# Create local VTK mesh data structures
topology, cell_types, geometry = plot.vtk_mesh(Vs)
num_cells_local = domain.topology.index_map(domain.topology.dim).size_local
num_dofs_local = Vs.dofmap.index_map.size_local * Vs.dofmap.index_map_bs
if topology.shape[0] > 0:
    num_dofs_per_cell = topology[0]
else:
    num_dofs_per_cell = 0
# Get only dof indices
topology_dofs = (np.arange(len(topology)) % (num_dofs_per_cell + 1)) != 0
# Map to global dof indices
global_dofs = Vs.dofmap.index_map.local_to_global(topology[topology_dofs].copy())
# Overwrite topology
topology[topology_dofs] = global_dofs
# Gather data
global_topology = domain.comm.gather(
    topology[: (num_dofs_per_cell + 1) * num_cells_local], root=0
)
global_geometry = domain.comm.gather(
    geometry[: Vs.dofmap.index_map.size_local, :], root=0
)
global_ct = domain.comm.gather(cell_types[:num_cells_local])

global_def = domain.comm.gather(uDs.x.array[:num_dofs_local])

if rank == 0:
    # Stack data

    root_geom = np.vstack(global_geometry)
    root_top = np.concatenate(global_topology)
    root_ct = np.concatenate(global_ct)

    root_def = np.concatenate(global_def)

    u_grid = pyvista.UnstructuredGrid(root_top, root_ct, root_geom)
    u_grid.point_data["u"] = root_def
    pyvista.start_xvfb()
    u_plotter = pyvista.Plotter(off_screen=True)
    u_plotter.add_mesh(u_grid, show_edges=True)

    u_plotter.show()
    u_plotter.screenshot(f"u_{domain.comm.rank}_{domain.comm.size}.png")

[Chaitanya_Apte]

Need some answers for better understanding:

1. what does this mean?? because the above code doesnt work on my dolfinx(0.8.0)

@dokken

nightly

2. The below link is a suggested feature right? so is it possible to use it in dolfinx(0.8)? or will it be added in future version?

@dokken

See: Add transfer_meshtags_to_submesh · Issue #3096 · FEniCS/dolfinx · GitHub which can transfer meshtags from parent to submesh.

[dokken]

@Chaitanya_Apte

what does this mean?? because the above code doesnt work on my dolfinx(0.8.0)

This means that you have to use the docker image i stated above, or install the main branch of dolfinx.

@Chaitanya_Apte

2. The below link is a suggested feature right? so is it possible to use it in dolfinx(0.8)? or will it be added in future version?

You can use the exact code from the issue to transfer data from parent to submesh. There is no functionality in that code relying on features from the main branch.