Duplicates in v2d map

mesh
1

Hello. I faced with the next problem: I generate the simple circle mesh using pygmsh library and save it in xml format. After that I try to load this mesh with Mesh dolfin object and create FunctionSpace. When I try to get vertex_to_dof_map, I faced with duplicates! Moreover, when I try to project some expression on this FunctionSpace, I get nan values in the tail of the projected function vector. The example below. Is it a correct behaviour? How could I get duplicates in vertex_to_dof_map? Thank you advance!

import meshio
import numpy as np
import pygmsh
from dolfin import *

geom = pygmsh.built_in.Geometry()  # Circle mesh generation
geom.add_circle(
    [0.0, 0.0, 0.0],
    1.0,
    lcar=0.1,
    num_sections=4,
    compound=True
)
mesh_data = pygmsh.generate_mesh(geom, prune_z_0=True)
meshio.write_points_cells("mesh.xml", mesh_data.points, mesh_data.cells)

mesh = Mesh("mesh.xml")
V = FunctionSpace(mesh, 'CG', 1)
v2d = vertex_to_dof_map(V)

v2d_len = v2d.shape[0]
unique_v2d_len = np.unique(v2d).shape[0]
print(v2d_len == unique_v2d_len)  # False -> v2d contains duplicates :(

exp = Expression("x[0] + x[1]", degree=1)
f = project(exp, V)
print(f.vector().get_local())  # nan values in the tail
2

Can you try to reproduce your example with FEniCS built-in meshes? There does not seem to be an issue there.

1 Like
3

I can’t reproduce this example with FEniCS built-in meshes, but may be somebody faced with such a problem with custom meshes and can give me advice :frowning:

4

I have found the solution. The reason why some elements in vector array was “nan” is that these elements was structural, generated by gmsh and not included in any cell. Solution of this problem was add “remove_faces=True” option in generate_mesh function. The complete code below.

import meshio
import numpy as np
import pygmsh
from dolfin import *

geom = pygmsh.built_in.Geometry()  # Circle mesh generation
geom.add_circle(
    [0.0, 0.0, 0.0],
    1.0,
    lcar=0.1,
    num_sections=4,
    compound=True,
)
mesh_data = pygmsh.generate_mesh(geom, prune_z_0=True, remove_faces=True)
meshio.write_points_cells("mesh.xml", mesh_data.points, mesh_data.cells)

mesh = Mesh("mesh.xml")
V = FunctionSpace(mesh, 'CG', 1)
v2d = vertex_to_dof_map(V)

v2d_len = v2d.shape[0]
unique_v2d_len = np.unique(v2d).shape[0]
print(v2d_len == unique_v2d_len)  # True -> v2d doesn't contains duplicates!

exp = Expression("x[0] + x[1]", degree=1)
f = project(exp, V)
print(f.vector().get_local())  # no one nan value
1 Like