I am trying to create a mesh comprising of three subdomains. On viewing the mesh with Gmsh it shows the correct rendering…but the PyVista engine is only tagging the surfaces and not the volumes. Not sure if it is a PyVista problem or the model generation.
import pyvista
from dolfinx.io import gmshio, XDMFFile
from dolfinx import mesh, fem, plot, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl,gmsh
import numpy as np
scale =1
L = 600.0*scale
W = 200.0*scale
H=50*scale
delta = W / L
airL=600*scale
airW=300*scale
airT=200*scale
airG=50*scale
model_rank = 0
mesh_comm = MPI.COMM_WORLD
def gmsh_cantilever (model:gmsh.model, name:str) ->gmsh.model:
airbox=model.occ.addBox(0,0,0,airL,airW,airT,1)
cantilever=model.occ.addBox(0,airW/2-W/2,airG+H,L,W,H,2)
base=model.occ.addBox(0,airW/2-W/2,0,L,W,H,3)
dom0,dom1=model.occ.fragment([(3,1)],[(3,i) for i in range(2,4)],removeTool=False,removeObject=True)
model.occ.removeAllDuplicates()
model.occ.synchronize()
model.addPhysicalGroup(3,[dom0[-1][1]],tag=1001)
model.setPhysicalName(3,1001,"Air")
model.addPhysicalGroup(3,[dom0[0][1]],tag=1002)
model.setPhysicalName(3,1002,"Cantilever")
model.addPhysicalGroup(3,[dom0[1][1]],tag=int(1003))
model.setPhysicalName(3,1003,"Base")
lcar1=0.1*scale
lcar2=0.05*scale
model.mesh.generate(dim=3)
gmsh.fltk.run()
return model
gmsh.initialize()
model=gmsh.model()
model =gmsh_cantilever(model,"Structure")
domain, ct, facet_tags = gmshio.model_to_mesh(model, mesh_comm, model_rank, gdim=3)
gmsh.finalize()
material_tags = np.unique(ct.values)
airregion = ct.find(1001)
beamregion = ct.find(1002)
baseregion = ct.find(1003)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_meshtags(ct, domain.geometry)
plotter = pyvista.Plotter()
tdim = domain.topology.dim
topology, cell_types, geometry = plot.vtk_mesh(domain)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
num_local_cells = domain.topology.index_map(tdim).size_local
grid.cell_data["Marker"] = ct.values[ct.indices]
grid.cell_data["Opacity"] = np.where(ct.values[ct.indices]==1001,0.5,1)
grid.set_active_scalars("Marker")
actor = plotter.add_mesh(grid, show_edges=True, opacity="Opacity", cmap="coolwarm")
plotter.view_xy()
plotter.show()