From gmsh to dolfin without file generation

CastriMik · January 19, 2022, 10:28am

Hi all,

For my simulations (with fenics 2019.2.0), I’m currently using gmsh for the mesh generation, and then I use meshio for transforming the .msh file into .xdmf files, and finally I load these files as MeshFunction. A typical workflow is:

from dolfin import *
import gmsh
import sys
import meshio
import numpy as np

gmsh.initialize(sys.argv)

if MPI.comm_world.rank == 0:

    gmsh.model.add("geometry")

    gmsh.model.occ.addCircle(0, 0, 0, 1, angle1 = 0, angle2 = 2*pi, tag=1)
    gmsh.model.occ.addCircle(0, 0, 0, 2, angle1 = 0, angle2 = 2*pi, tag=2)

    gmsh.model.occ.addCurveLoop([1], tag=1)
    gmsh.model.occ.addPlaneSurface([1], tag=1)

    gmsh.model.occ.addCurveLoop([2], tag=2)
    gmsh.model.occ.addCurveLoop([1], tag=3)
    gmsh.model.occ.addPlaneSurface([2, 3], tag=2)

    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(2, [1], tag=1)
    gmsh.model.addPhysicalGroup(2, [2], tag=2)
    gmsh.model.addPhysicalGroup(1, [2], tag=3)

    gmsh.model.mesh.setSize([(0, 1)], size = 0.1)
    gmsh.model.mesh.setSize([(0, 2)], size = 0.1)

    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")

    gmsh.finalize()

if MPI.comm_world.rank == 0:

    msh = meshio.read("mesh.msh")

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type == "triangle"]))

    total_mesh = meshio.Mesh(points=msh.points,
                             cells=[("triangle", cells)])
    total_mesh.prune_z_0()

    cells_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

    cells_mesh = meshio.Mesh(points=msh.points,
                                cells=[("triangle", cells)],
                                cell_data={"name_to_read": [cells_data]})

    cells_mesh.prune_z_0()
    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                      if cells.type == "line"]))

    facet_data = msh.cell_data_dict["gmsh:physical"]["line"]

    facet_mesh = meshio.Mesh(points=msh.points,
                             cells=[("line", facet_cells)],
                             cell_data={"name_to_read": [facet_data]})

    facet_mesh.prune_z_0()

    meshio.write("total.xdmf", total_mesh)
    meshio.write("facets.xdmf", facet_mesh)
    meshio.write("cells.xdmf", cells_mesh)

MPI.comm_world.barrier()
parameters["ghost_mode"] = "shared_vertex"        

comm = MPI.comm_world
mesh = Mesh(MPI.comm_world)

with XDMFFile(MPI.comm_world, "cells.xdmf") as infile:
    infile.read(mesh)

mvc_facets = MeshValueCollection("size_t", mesh, 1)

with XDMFFile(MPI.comm_world, "facets.xdmf") as infile:
    infile.read(mvc_facets, "name_to_read")

mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)

mvc_cells = MeshValueCollection("size_t", mesh, 2)

with XDMFFile(MPI.comm_world, "cells.xdmf") as infile:
    infile.read(mvc_cells, "name_to_read")

mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)

I’m wondering if someone has already worked out a way to shortcut the mesh loading into MeshFunction from gmsh without having to write .msh and .xdmf file. Maybe something can be done with MeshEditor, as suggested here, but I’m having troubles in understanding how I should translate this code into a gmsh compatible code.

dokken · January 19, 2022, 11:29am

You can extract the topology (cells) and geometry (points) of the gmsh object as done in:

github.com

FEniCS/dolfinx/blob/main/python/dolfinx/io.py#L81-L151

    
      
          def extract_gmsh_topology_and_markers(gmsh_model, model_name=None):
              """Extract all entities tagged with a physical marker in the gmsh
              model, and collects the data per cell type. Returns a nested
              dictionary where the first key is the gmsh MSH element type integer.
              Each element type present in the model contains the cell topology of
              the elements and corresponding markers.
          
          
    """
              if model_name is not None:
                  gmsh_model.setCurrent(model_name)
              # Get the physical groups from gmsh on the form
              # [(dim1, tag1),(dim1, tag2), (dim2, tag3),...]
              phys_grps = gmsh_model.getPhysicalGroups()
              topologies = {}
              for dim, tag in phys_grps:
                  # Get the entities for a given dimension:
                  # dim=0->Points, dim=1->Lines, dim=2->Triangles/Quadrilaterals,
                  # etc.
                  entities = gmsh_model.getEntitiesForPhysicalGroup(dim, tag)

This file has been truncated. show original

CastriMik · January 21, 2022, 11:21am

Thank you very much for the suggestion.

I think I’ve almost done it, but I’m having hard times in understanding how to correctly write the MeshValueCollection without the use of .xdmf file (Eventually, I’d like to have MeshFunction variables to use as subdomain_data for my Measure). My code is the following one at the moment:

from dolfin import *
import gmsh
import sys
import numpy as np
import dolfin as df
from ufl import vertex

def extract_gmsh_topology_and_markers(gmsh_model, model_name=None):
    if model_name is not None:
        gmsh_model.setCurrent(model_name)
    phys_grps = gmsh_model.getPhysicalGroups()
    topologies = {}
    for dim, tag in phys_grps:
        entities = gmsh_model.getEntitiesForPhysicalGroup(dim, tag)
        for entity in entities:
            element_data = gmsh_model.mesh.getElements(dim, tag=entity)
            element_types, element_tags, node_tags = element_data
            assert len(element_types) == 1
            element_type = element_types[0]
            num_el = len(element_tags[0])
            properties = gmsh_model.mesh.getElementProperties(element_type)
            name, dim, order, num_nodes, local_coords, _ = properties
            element_topology = node_tags[0].reshape(-1, num_nodes) - 1
            if element_type in topologies.keys():
                topologies[element_type]["topology"] = np.concatenate(
                    (topologies[element_type]["topology"], element_topology), axis=0)
                topologies[element_type]["cell_data"] = np.concatenate(
                    (topologies[element_type]["cell_data"], np.full(num_el, tag)), axis=0)
            else:
                topologies[element_type] = {"topology": element_topology, "cell_data": np.full(num_el, tag)}
    return topologies

def extract_gmsh_geometry(gmsh_model, model_name=None):
    if model_name is not None:
        gmsh_model.setCurrent(model_name)
    indices, points, _ = gmsh_model.mesh.getNodes()
    points = points.reshape(-1, 3)
    indices -= 1
    perm_sort = np.argsort(indices)
    assert np.all(indices[perm_sort] == np.arange(len(indices)))
    return points[perm_sort]

comm = MPI.comm_world

gmsh.initialize(sys.argv)

if comm.rank == 0:

    gmsh.model.add("geometry")

    gmsh.model.occ.addCircle(0, 0, 0, 1, angle1 = 0, angle2 = 2*pi, tag=1)
    gmsh.model.occ.addCircle(0, 0, 0, 2, angle1 = 0, angle2 = 2*pi, tag=2)

    gmsh.model.occ.addCurveLoop([1], tag=1)
    gmsh.model.occ.addPlaneSurface([1], tag=1)

    gmsh.model.occ.addCurveLoop([2], tag=2)
    gmsh.model.occ.addCurveLoop([1], tag=3)
    gmsh.model.occ.addPlaneSurface([2, 3], tag=2)

    gmsh.model.occ.synchronize()

    gmsh.model.addPhysicalGroup(2, [1], tag=1)
    gmsh.model.addPhysicalGroup(2, [2], tag=2)
    gmsh.model.addPhysicalGroup(1, [2], tag=3)

    gmsh.model.mesh.setSize([(0, 1)], size = 0.1)
    gmsh.model.mesh.setSize([(0, 2)], size = 0.1)

    gmsh.model.mesh.generate(2)

    topologies = extract_gmsh_topology_and_markers(gmsh.model)
    cells = topologies[2]["topology"]
    points = extract_gmsh_geometry(gmsh.model)
    points = points[:, 0:2]
    parameters["ghost_mode"] = "shared_vertex"        

else:
    points = cells = None

comm.barrier()

points = comm.bcast(points, root=0)
cells  = comm.bcast(cells, root=0)

editor = MeshEditor()
mesh = Mesh(comm)
editor.open(mesh, "triangle", 2, 2)

if comm.rank == 0:
   
    editor.init_vertices(points.shape[0])
    editor.init_cells(cells.shape[0])
    for k, point in enumerate(points):
        editor.add_vertex(k, point)
    for k, cell in enumerate(cells):
        editor.add_cell(k, cell)
    editor.close()

comm.barrier()

print(assemble(Constant(1)*dx(domain=mesh)))

#with XDMFFile(comm, "cells.xdmf") as infile:
#    infile.read(mesh)
#
#mvc_facets = MeshValueCollection("size_t", mesh, 1)
#
#with XDMFFile(comm, "facets.xdmf") as infile:
#    infile.read(mvc_facets, "name_to_read")
#
#mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)
#
#mvc_cells = MeshValueCollection("size_t", mesh, 2)
#
#with XDMFFile(comm, "cells.xdmf") as infile:
#    infile.read(mvc_cells, "name_to_read")
#
#mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)

The commented part is the one I wish to translate into an “xdmf-free” code, but, as I said, I was not able to find a way to define my MeshValueCollection w/o .xdmf files.

CastriMik · January 31, 2022, 9:12am

Hello everyone,

I’m still struggling finding a way to express MeshFunction without passing through .xdmf files.

Could anyone give me some hints on that? Which dolfin functions shoud I look at?

dokken · January 31, 2022, 10:16am

As far as I am aware, you need: Bitbucket
which is similar to dolfinx/xdmf_utils.h at fe741ee468b43bd4645202201aa983726eed59ab · FEniCS/dolfinx · GitHub in DOLFINx, which is used to map the gmsh entities to the correct local entity: dolfinx/demo_gmsh.py at 3a113a8854b14649a70127cae6455f2588e8be31 · FEniCS/dolfinx · GitHub

In dolfin, this function is called when reading in XDMFFiles:Bitbucket

CastriMik · February 10, 2022, 4:04pm

Thanks so much, I’ll have a look at it!

Topic		Replies	Views
Transitioning from Triangle -> GMSH -> XDMF I/O	8	913	April 16, 2023
Gmsh, Meshio, Dolfin mesh	5	1692	June 3, 2023
Create mesh with gmsh, convert it and import it in dolfin, problems mesh	5	10954	April 1, 2019
Using a simple 3d mesh from gmsh General	9	2333	November 19, 2022
Mesh 2D gmsh to fenics legacy (2019.2) mesh	3	111	September 27, 2024

From gmsh to dolfin without file generation

Related topics