Why is mass matrix singular After Converting Gmesh mesh to XDMF?

mesh
1

Hi everyone, I am writing because after creating my first using Gmesh, converting the mesh to XDMF format, and constructing a finite element space corresponding to the mesh, I am noticing that the mass matrix corresponding to the space is singular. I do not know if the issue is with Gmesh or with the conversion to XDMF file, so I am including all of my code for reference.

First, I created the mesh using the Gmesh Python API (version 4.11.1) with the following code:

import gmsh
import math
import sys
import os

gmsh.initialize()

gmsh.model.add("Second Mesh")

factory = gmsh.model.geo

factory.add_point(-5, 0, 0, .1, 1)
factory.add_point(5, 0, 0, .1, 2)
factory.add_point(0, -5, 0, .1, 3)
factory.add_point(0, 5, 0, 1, 4)
factory.add_point(0, 0, 0, 1, 5)

factory.add_circle_arc(2, 5, 3, 1)
factory.add_circle_arc(3, 5, 1, 2)
factory.add_circle_arc(1, 5, 4, 3)
factory.add_circle_arc(4, 5, 2, 4)

factory.add_curve_loop([1, 2, 3, 4], 5)

factory.add_plane_surface([5], 6)

factory.synchronize()

gmsh.model.mesh.field.add("Ball", 7)
gmsh.model.mesh.field.setNumber(7, "Radius", 4)
gmsh.model.mesh.field.setNumber(7, "VIn", .15)
gmsh.model.mesh.field.setNumber(7, "Thickness", .3)
gmsh.model.mesh.field.setNumber(7, "VOut", .08)

gmsh.model.mesh.field.setAsBackgroundMesh(7)

gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)


factory.synchronize()
gmsh.model.mesh.generate(2)

gmsh.write("BallMesh3.msh")

# Launch the GUI to see the results:
gmsh.fltk.run()

gmsh.finalize()

The following code is a minimal example that shows that after conversion to XDMF format, the mass matrices of finite element spaces corresponding to the mesh are singular.

import meshio
import numpy as np
import dolfin as dl
import ufl

# load mesh and convert to xdmf file
mesh_name = 'BallMesh3.msh'
dim = 2
directory = 'meshes'
prefix = mesh_name.split('.')[0]
msh = meshio.read("{}/{}".format(directory, mesh_name))
if dim == 2:
    cell_type = "triangle"
elif dim == 3:
    cell_type = "tetra"
data_array = [cell.data for cell in msh.cells if cell.type == cell_type]
if len(data_array) == 0:
    print("WARNING: No domain physical group found.")
else:
    data = np.concatenate(data_array)
cells = [
    meshio.CellBlock(
        cell_type=cell_type,
        data=data,
    )
]
domain = meshio.Mesh(
    points=msh.points[:, :dim],
    cells=cells)
meshio.write(
    "{}/{}_{}".format(directory, prefix, "domain.xdmf"),
    domain,
    file_format="xdmf"
)

# load xdmf file, create function space
mesh = dl.Mesh()
with dl.XDMFFile("meshes/BallMesh3_domain.xdmf") as fid:
    fid.read(mesh)
Vh = dl.FunctionSpace(mesh, "CG", 1)
trial = dl.TrialFunction(Vh)
test  = dl.TestFunction(Vh)
varfM = ufl.inner(trial,test)*ufl.dx       
M = dl.assemble(varfM).array()

matrix_rank = np.linalg.matrix_rank(M)
matrix_size = len(M)
print("The matrix rank is " + str(matrix_rank))
print("The matrix size is " + str(matrix_size))

Here I am using Python 3.8.16, Fenics 2019.1.0, and MeshIo 5.3.4.

Even if you do not know the cause of this bug but have ideas as to how to isolate the issue, please share. Thank you in advance for your help.

2

I discovered that with the code as it is written, an extra point at the center of the circle (i.e., the origin) was added to the mesh. I fixed the issue by adding a physical group for the surface in the mesh creation, after which the list of points was reduced by one. Keeping this here in case anyone else runs in to the same problem.

1 Like