Adding 1d physical groups to gmsh model causes errors when convert it to dolfinx mesh

mesh
1

Hello, I wrote a function to generate a 2d area by gmsh. But there are some problems in the generate_mesh function.

  1. If I only add the region group or region + parabola group, the code runs correctly. It prints 1/6 as I expected.
  2. If I add all the groups (region + parabola +x_axis) or I add the groups region + x_axis , MPI error occurs.
# run the python code with mpirun -n 2 
Invalid rank, error stack:
internal_Issend(61644): MPI_Issend(buf=0x7fffe0320db2, count=1, MPI_BYTE, 2, 1, comm=0xc400000c, request=0x7fffe02f3e38) failed
internal_Issend(61605): Invalid rank has value 2 but must be nonnegative and less than 2
Invalid rank, error stack:
internal_Issend(61644): MPI_Issend(buf=0x7fffe4ae99c2, count=1, MPI_BYTE, 2, 1, comm=0xc4000005, request=0x7fffe4b0dc08) failed
internal_Issend(61605): Invalid rank has value 2 but must be nonnegative and less than 2
  1. If I only add the group parabola or x_axis, there is no error, but the code prints strange values.

Does anyone have some ideas about this ? Thanks.
Here is the code

import gmsh
import dolfinx, ufl
import numpy as np
from mpi4py import MPI
comm = MPI.COMM_WORLD
def generate_mesh(grid_size, order:int=1):
    gmsh.initialize()
    gmsh.model.add("parabola")

    nx = 50
    x_array = np.linspace(0, 1, nx)
    points = [gmsh.model.occ.add_point(x, x-x**2, 0, grid_size) for x in x_array]
    parabola = gmsh.model.occ.add_spline(points)

    p0 = gmsh.model.occ.add_point(0, 0, 0.0, grid_size)
    p1 = gmsh.model.occ.add_point(1.0, 0, 0.0, grid_size)
    line_x_axis = gmsh.model.occ.add_line(p1, p0)

    curves = [parabola, line_x_axis]
    curve_loop = gmsh.model.occ.addCurveLoop(curves)
    surface = gmsh.model.occ.add_plane_surface([curve_loop])
    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(2, [surface], name="region")
    gmsh.model.addPhysicalGroup(1, [parabola], name="parabola")
    #gmsh.model.addPhysicalGroup(1, [line_x_axis], name="x_axis")

    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", grid_size)
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.setOrder(order)
    mesh, cell_marker, facet_marker = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
    gmsh.finalize()

    return mesh, cell_marker, facet_marker


mesh, cell_marker, facet_marker = generate_mesh(0.05, 2)
dx = ufl.Measure("dx", domain=mesh)
local_area = dolfinx.fem.assemble_scalar(dolfinx.fem.form(1.0*dx))
global_area = mesh.comm.allreduce(local_area, op=MPI.SUM)
if comm.Get_rank() == 0:
    print(global_area)
2

This is because line_x_axis is not part of the surface, as seen by calling:

    boundary = gmsh.model.getBoundary([(2, surface)])
    print(f"{boundary=}, {line_x_axis=}")

Thus, if you instead add:

    for boundary_dim, boundary_tag in boundary:
        gmsh.model.addPhysicalGroup(boundary_dim, [boundary_tag], tag=boundary_tag)

you will get all your outer boundaries tagged.

2 Likes