Visualizing mesh volume regions from GMSH labels

mesh
1

I would like to visualize the volume regions as labeled by GMSH. This is to verify that all material parameters are correctly applied. A simple toy problem illustrates what I am trying to do.

GMSH geo file:

SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 3, 3, 1};
Box(2) = {0, 0, 1, 3, 3, 1};
Box(3) = {0, 0, 2, 3, 3, 1};

BooleanFragments{Volume{1,2,3}; Delete;}{};
Physical Volume("one", 1) = {1};
Physical Volume("two", 2) = {2};
Physical Volume("three", 3) = {3};

Mesh.Algorithm = 4;
Mesh.CharacteristicLengthMin = 0.5;

Each of three regions has a unique label.

Now I use Meshio to convert msh file to xdmf.

import meshio
msh = meshio.read("TestVolumes.msh")
for cell in msh.cells:
    if  cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]

tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                           cell_data={"VolumeRegions":[tetra_data]})

meshio.write("mesh.xdmf", tetra_mesh)

So far, so good. This runs without error. I now want to import the geometry into Fenics using

from dolfin import *

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, 3)
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
    
info(mesh)

This runs without error, but when I use

File("VolSubDomains.pvd").write(cf)
File("Mesh.pvd").write(mesh)

I do not see the regions 1, 2 or 3 when I plot the resulting pvd files with Paraview. I only see the full mesh as one unified volume.

  1. What am I doing wrong?
  2. When building the variational form, is it enough to use dx(1), dx(2), dx(3), etc.?
2

Your are not loading in your cell-data to fenics. See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio, where I describe how to load in cell functions, as well as how to use them in the dx measure.

1 Like
3

That’s the bit I was missing! Thanks!