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!
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, [volumes], marker) gmsh.model.setPhysicalName(volumes, 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)) values = -10.0 values = 10.0 if dim == 3: values = 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