Export to gmsh: Connectivity issues

dolfinx
1

Hello,
I am conducting explicit dynamic simulations and aim to remesh my domain and interpolate the previous solution when convergence issues arise. To this end, I write the last deformed mesh in the .msh format using Meshio. This mesh is then given as an argument to gmsh for remeshing, after which I re-export it in a format compatible with DolfinX.

However, I suspect that Meshio/gmsh disrupts the connectivity, resulting in meshes that are visually appealing yet functionally useless after remeshing.

messed_up_mesh
messed_up_mesh1542×784 76 KB

Below, you will find a Minimal Working Example (MWE) that illustrates this phenomenon:

from mpi4py import MPI
from dolfinx.fem import (FunctionSpace, Function, VectorFunctionSpace)
from dolfinx.io import gmshio, XDMFFile
import matplotlib.pyplot as plt
import gmsh
import meshio
from dolfinx.io.gmshio import model_to_mesh

#%% Mesh generation
def init_gmsh():
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
    gdim = 2
    rank = 0
    gmsh.model.add("Model")
    geom = gmsh.model.geo
    return gdim, rank, geom

def return_mesh(model, comm, rank, gdim, quad = True, degree = 1):
    gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
    if quad:
        gmsh.option.setNumber("Mesh.Algorithm", 8)
        gmsh.option.setNumber("Mesh.RecombineAll", 1)
        gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.setOrder(degree)
    gmsh.model.mesh.optimize("Netgen")

    domain, meshtags, facets = model_to_mesh(model, comm, rank, gdim=gdim)
    gmsh.finalize()
    return domain, meshtags, facets

def broken_unit_square(N_largeur, Nh, Largeur = 1, tol = None, raff = 0):
    if tol == None:
        tol = Largeur / 1e3
    gdim, model_rank, geom = init_gmsh()
    h_size = Largeur/N_largeur
    geom.synchronize()
    
    pb1 = geom.add_point(0, 0, 0)
    pb2 = geom.add_point(Largeur/2, 0, 0)
    pb3 = geom.add_point(Largeur, 0, 0)
    
    pt1 = geom.add_point(0, Largeur, 0)
    pt2 = geom.add_point(Largeur/2, Largeur, 0)
    pt3 = geom.add_point(Largeur, Largeur, 0)
    
    pm1 = geom.add_point(0, Largeur/2 - tol, 0)
    pm1bis = geom.add_point(0, Largeur/2 + tol, 0)
    pm2 = geom.add_point(Largeur/2, Largeur/2, 0)
    pm3 = geom.add_point(Largeur, Largeur/2, 0)

    left_bottom = geom.add_line(pm1, pb1)
    bottom_left = geom.add_line(pb1, pb2)
    bottom_right = geom.add_line(pb2, pb3)
    right_bottom = geom.add_line(pb3, pm3)
    right_mid = geom.add_line(pm3, pm2)
    bottom_lip = geom.add_line(pm2, pm1)
    
    left_top = geom.add_line(pm1bis, pt1)
    top_left = geom.add_line(pt1, pt2)
    top_right = geom.add_line(pt2, pt3)
    right_top = geom.add_line(pt3, pm3)
    top_lip = geom.add_line(pm2, pm1bis)
    
    boundary_1 = geom.add_curve_loop([left_bottom, bottom_left, bottom_right, right_bottom, right_mid, bottom_lip])
    surf_1 = geom.add_plane_surface([boundary_1])
    
    boundary_2 = geom.add_curve_loop([left_top, top_left, top_right, right_top, right_mid, top_lip])
    surf_2 = geom.add_plane_surface([boundary_2])
    geom.synchronize()

    gmsh.model.addPhysicalGroup(gdim, [surf_1], 1)
    gmsh.model.addPhysicalGroup(gdim, [surf_2], 2)
    gmsh.model.addPhysicalGroup(gdim - 1, [bottom_left, bottom_right], 1, name="bottom")
    gmsh.model.addPhysicalGroup(gdim - 1, [left_bottom, left_top], 2, name="left")
    gmsh.model.addPhysicalGroup(gdim - 1, [top_left, top_right], 3, name="top")
    gmsh.model.addPhysicalGroup(gdim - 2, [pb1], 4, name="corner")

    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h_size)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h_size)
        
    geom.synchronize()
    return return_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim)

#%% Mesh raff
def dolfinx_mesh_to_msh_with_meshio(input_mesh, name = "mesh.msh"):
    ufl_celltype = input_mesh.ufl_cell().cellname()
    if ufl_celltype =="triangle":
        celltype = "triangle"
        reshape_size = 3
    elif ufl_celltype =="quadrilateral":
        celltype = "quad"
        reshape_size = 4     
    cells = input_mesh.topology.connectivity(input_mesh.topology.dim, 0).array.reshape((-1, reshape_size))
    cell_type = meshio.CellBlock(celltype, cells)
    points = input_mesh.geometry.x    
    # Write .msh
    mio_mesh = meshio.Mesh(points, [cell_type])
    meshio.write(name, mio_mesh, file_format="gmsh22")
    return name

def improved_gmsh_mesh_to_dolfinx_mesh(input_mesh_name):
    gdim, model_rank, geom = init_gmsh()
    gmsh.open(input_mesh_name)
    gmsh.fltk.run()
    gmsh.model.mesh.optimize("Netgen")
    tags_surface = [entite[1] for entite in gmsh.model.getEntities(gdim)]
    tag_groupe_physique = gmsh.model.addPhysicalGroup(gdim, tags_surface)
    domain, meshtags, facets = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim)
    gmsh.finalize()
    return domain, meshtags, facets
    
def improve_mesh(domain, name = "old_mesh"):
    return improved_gmsh_mesh_to_dolfinx_mesh(dolfinx_mesh_to_msh_with_meshio(domain, name))

######## Mesh ########
domain, _, facet_tag = broken_unit_square(10, 10)
def dummy_deformation(x):
    return (x[0]-x[1]**2, x[1])

V = VectorFunctionSpace(domain, ("Lagrange", 1))
u_radial = Function(V)
u_radial.interpolate(dummy_deformation)
with XDMFFile(domain.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)

def deform_mesh(mesh, displacement_array):
    deformation_array = displacement_array.x.array.reshape((-1, mesh.geometry.dim))
    mesh.geometry.x[:, :mesh.geometry.dim] += deformation_array
    
deform_mesh(domain, u_radial)
with XDMFFile(domain.comm, "mesh_deformed.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    
new_domain, _, facet_tag = improve_mesh(domain)

with XDMFFile(domain.comm, "new_mesh_deformed.xdmf", "w") as xdmf:
    xdmf.write_mesh(new_domain)

I deform my non-trivial mesh with an arbitrary displacement field and then attempt to remesh it.
Would anyone have a suggestion for solving this problem? Thanks in advance.
Sincerely,
Paul

2

This is not correct. Here you assume a 1-1 correspondence between the topology and geometry, which is not correct.

You need to use the function entities_to_geometry to map each cell to its appropriate geometry indices (vertices in the case of a linear mesh).

Note that entities_to_geometry do not support higher order geometries, and you would have to use something like: Obtain high order mesh edge and face nodes? - #4 by August

3

Thank you very much @dokken for your quick response. If I understand correctly, entities_to_geometry returns the indices of the cells along with the coordinates of the nodes constituting them. Consequently, I have modified the function dolfinx_mesh_to_msh_with_meshio as follows:

def dolfinx_mesh_to_msh_with_meshio(input_mesh, name = "mesh.msh"):
    ufl_celltype = input_mesh.ufl_cell().cellname()
    if ufl_celltype =="triangle":
        celltype = "triangle"
    elif ufl_celltype =="quadrilateral":
        celltype = "quad"
    gdim = input_mesh.topology.dim
    num_cells_owned_by_proc = input_mesh.topology.index_map(gdim).size_local
    cell_entitites = entities_to_geometry(input_mesh._cpp_object, gdim, np.arange(num_cells_owned_by_proc, dtype=np.int32), False)
    cell_type = meshio.CellBlock(celltype, cell_entitites)
    points = input_mesh.geometry.x
    mio_mesh = meshio.Mesh(points, [cell_type])
    meshio.write(name, mio_mesh, file_format="gmsh22")
    return name

However, this modification does not seem to alter the outcome. What am I still missing there?
Thank you

4

If you work a bit with simplifying your code to a flatter structure, it is more likely that I have time to help you.
Currently, you have alot of nested function calls, that takes time for me to go through.

5

Yes, of course, my apologies @dokken . I had the impression that the error occurred only for non-trivial meshes, but even with the built-in meshes, my code yields exotic results.
EDIT: I believe my issue can be more succinctly reformulated with the following question: How can one properly export a mesh from DOLFINx in order to be readable by GMSH? Even without any mesh deformation, the following script completely disrupts the mesh connectivity. Has anyone encountered this issue before and would know which arguments to provide to Meshio to ensure a correct export?
This time, here is a much shorter code illustrating my problem:

from mpi4py import MPI
from dolfinx.io import XDMFFile
import gmsh
import meshio
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.cpp.mesh import entities_to_geometry
import numpy as np
from dolfinx.mesh import create_unit_square, CellType

######## Mesh ########
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, CellType.quadrilateral)
with XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    
name = "mesh.msh"
gdim = mesh.topology.dim
num_cells_owned_by_proc = mesh.topology.index_map(gdim).size_local
cell_entitites = entities_to_geometry(mesh._cpp_object, gdim, np.arange(num_cells_owned_by_proc, dtype=np.int32), False)
cell_type = meshio.CellBlock("quad", cell_entitites)
points = mesh.geometry.x
mio_mesh = meshio.Mesh(points, [cell_type])
meshio.write(name, mio_mesh, file_format="gmsh22")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
gmsh.model.add("Model")
gmsh.open(name)
gmsh.fltk.run()
tags_surface = [entite[1] for entite in gmsh.model.getEntities(gdim)]
tag_groupe_physique = gmsh.model.addPhysicalGroup(gdim, tags_surface)
new_mesh, meshtags, facets = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim)
gmsh.finalize()

with XDMFFile(new_mesh.comm, "new_mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(new_mesh)

Many thanks, sincerely,
Paul

6

You are missing a permutation from DOLFINx node ordering to GMSH.
I’ve added this below:

from mpi4py import MPI
from pathlib import Path
from dolfinx.io import XDMFFile
import gmsh
import meshio
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.cpp.mesh import entities_to_geometry
from dolfinx.cpp.io import perm_gmsh
import numpy as np
from dolfinx.mesh import create_unit_square, CellType

######## Mesh ########
mesh = create_unit_square(MPI.COMM_WORLD, 10, 10, CellType.quadrilateral)
with XDMFFile(mesh.comm, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

name = Path("tmp_mesh.msh")
tdim = mesh.topology.dim
num_cells_owned_by_proc = mesh.topology.index_map(tdim).size_local
cell_entitites = entities_to_geometry(mesh._cpp_object, tdim, np.arange(
    num_cells_owned_by_proc, dtype=np.int32), False)
permutation = np.asarray(perm_gmsh(
    mesh.topology.cell_types[0], cell_entitites.shape[1]), dtype=np.int32)

cell_entitites = cell_entitites[:, permutation]
cell_type = meshio.CellBlock("quad", cell_entitites)
points = mesh.geometry.x
mio_mesh = meshio.Mesh(points, [cell_type])
meshio.write(name.with_suffix(".msh"), mio_mesh, file_format="gmsh22")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # to disable meshing info
gmsh.model.add("Model")
gmsh.open(str(name))

tags_surface = [entite[1] for entite in gmsh.model.getEntities(tdim)]
tag_groupe_physique = gmsh.model.addPhysicalGroup(tdim, tags_surface)
new_mesh, meshtags, facets = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, tdim)
gmsh.finalize()

with XDMFFile(new_mesh.comm, "new_mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(new_mesh)
7

That’s exactly what i was looking for. Thanks a lot !