Import of gmsh mesh leaves nodes empty / misses basis function

When working with gmsh to create some rather basic geometries I encountered the problem that at some nodes one misses the basis function apparently. At least the interpolation does not produce any observable vectors at these nodes in paraview. However paraview also does not consider these nodes to have a zero value. In some more complicated code example which I do not share here, This also leads to unexpected behaviour, therefore I am quite certain that the basis functions at these nodes are in some sense missing in the FEM space. I tried several geometries and algorithms in gmsh so far. I also tried removing the duplicates and saving the mesh as .msh file before importing it with dolfinx.io again which leads to the same results. Any help would be appreciated!

MWE:

from mpi4py import MPI
import numpy as np
import gmsh
from dolfinx import *


def create_2D_mesh(algo = 1):
    # Using gmsh    
    gmsh.initialize()
    gmsh.clear() # should not be necessary if only called once

    gdim = 2
    c1 = gmsh.model.occ.addDisk(0, 0, 0, 2, 2)  
    gmsh.model.occ.synchronize()

    # add volume marker
    marker = 1
    volumes = gmsh.model.getEntities(dim=gdim)
    assert(len(volumes) == 1)
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], marker)
    gmsh.model.setPhysicalName(volumes[0][0], marker, "full")

    gmsh.option.setNumber("Mesh.Algorithm", algo)
    gmsh.option.set_number("Geometry.Tolerance", 1e-3)
    gmsh.model.geo.remove_all_duplicates()
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(dim = gdim )
    return gmsh.model

def ic(x: np.ndarray) -> np.ndarray:
    dim = 2
    # x hase shape (dimension, points)
    values = np.zeros((dim, x.shape[1]))
    values[0] = -10.0
    values[1] = 10.0
    if dim == 3: values[2] = 0.0
    return values
  
if __name__ == "__main__":
    dim = 2
    # create gmsh model
    model = create_2D_mesh(1)
    # extract mesh
    domain, cell_tags, meshtags = io.gmshio.model_to_mesh(model,  MPI.COMM_WORLD, 0, gdim=dim)

    # CREATE FEM FUNCTION SPACES    
    V       = fem.VectorFunctionSpace(domain, ("Lagrange", 1), dim = dim)
    u       = fem.Function(V)

    u.interpolate(ic) 

    with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
        xdmf.write_mesh(domain)
        xdmf.write_function(u)

Plot created with paraview

This Seems to be due to your visualization in Paraview, as you can choose how many glyphs to display.
Ref image below:

Thank you very much! That error source, I most certainly did not expect…