How to build a quarter of a disk (easy mesh question)

mesh
1

Hi,
(I’m using a Docker compose with fenicsx and pyvista on my windows laptop)

I’m trying to simulate a PDE on a domain that is the first quarter of a disk (from 0 to pi/2 radians).
I built the mesh and everything went fine but pyvista couldn’t plot it, even though it can plot easily every mesh in the tutorial.
I built the mesh in the following way:

gmsh.model.add("diffusion_equation")

gmsh.model.occ.addPoint(0,0,0)#point1
gmsh.model.occ.addPoint(1,0,0)#point2
gmsh.model.occ.addPoint(0,1,0)#point3

gmsh.model.occ.addLine(1,2,1)#line1
gmsh.model.occ.addCircleArc(2,1,3,2)#line2
gmsh.model.occ.addLine(3,1,3)#line3

gmsh.model.occ.addCurveLoop([1,2,3],1)#loop1

#gmsh.model.occ.addSurfaceFilling(1,1)#plane1
gmsh.model.occ.addPlaneSurface([1],1)#plane1

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(1,[1],101)
gmsh.model.setPhysicalName(1,101, "t-p")
gmsh.model.addPhysicalGroup(1,[2],102)
gmsh.model.setPhysicalName(1,102, "t-a")
gmsh.model.addPhysicalGroup(1,[3],103)
gmsh.model.setPhysicalName(1,103, "t-c")
gmsh.model.addPhysicalGroup(2,[1],201)
gmsh.model.setPhysicalName(2,201, "b")

gdim = 2 #
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.05)

gmsh.model.mesh.generate(gdim)

filename = 'quarter_of_circle.msh'
gmsh.write(filename)

from dolfinx import io, mesh, cpp, fem
from mpi4py import MPI
import numpy as np
if MPI.COMM_WORLD.rank == 0:
    # Get mesh geometry
    geometry_data = io.extract_gmsh_geometry(gmsh.model)
    # Get mesh topology for each element
    topology_data = io.extract_gmsh_topology_and_markers(gmsh.model)


if MPI.COMM_WORLD.rank == 0:
    # Extract the cell type and number of nodes per cell and broadcast
    # it to the other processors 
    gmsh_cell_type = list(topology_data.keys())[0]    
    properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type)
    name, dim, order, num_nodes, local_coords, _ = properties
    cells = topology_data[gmsh_cell_type]["topology"]
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)
else:        
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
    cells, geometry_data = np.empty([0, num_nodes]), np.empty([0, gdim])
    
# Permute topology data from MSH-ordering to dolfinx-ordering
ufl_domain = io.ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = io.cell_perm_gmsh(cpp.mesh.to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]

# Create distributed mesh
domain = mesh.create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :gdim], ufl_domain)

from dolfinx import plot
import pyvista
tdim = domain.topology.dim
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

pyvista.set_jupyter_backend("pythreejs")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    figure = plotter.screenshot("fundamentals_mesh.png")

By looking at the gmsh file in the gmsh GUI it seems fine but pyvista can’t plot it.

Am I doing something wrong?

Thanks in advance!

2

There are several things that I find curious in your code.

should be:

    gmsh_cell_type = list(topology_data.keys())[-1]

as you want to extract the cell type for the cells, not the lowest tagged dimension (facets/intervals).

1 Like
3

Thanks for the help! I used a cell from your tutorial but I shouldn’t since in fact the problem is different, I need facets to impose boundary conditions. Is there anything else that’s strange in my code and I should look up?
Thanks again

4

See: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

1 Like