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
