Assigning different physical groups to different color in gmsh and mark them as different material in fenics

Hello, I am stuck at how to assign different physical groups to different colors shown in image attached below, basically I am trying to given them different material properties in fenics to the physical tags afterwards but I am also confused on how to do that since I don’t know the geometrical locations of these different radii seeds, and I want to export only mesh part in .msh file and use it in fenics which is also attached below. Therefore my humble request to the community is to kindly guide me through this.
Thank You


Consider How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

Thanks @dokken for the reply, but I am confused that how do I define physical group which is already been defined in this forum and in my case i don’t know how to add physical tags, since they are randomly generated seeds.

@dokken I tried the code but as I was thinking about the physical tags, that came as an error when I implemented the below code by importing my .msh file.

from dolfin import *
import meshio
msh = meshio.read("/home/ratan/Documents/gmsh/mp1.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)

The error which came:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-11-1b306cf427a8> in <module>
      9         tetra_cells = cell.data
     10 
---> 11 for key in msh.cell_data_dict["gmsh:physical"].keys():
     12     if key == "triangle":
     13         triangle_data = msh.cell_data_dict["gmsh:physical"][key]

KeyError: 'gmsh:physical'

Kindly suggest me on how do solve this issue.
Thank You

If you dont have the physical tags of each seed when generating the mesh, how should one distinguish between them? Where did you get the coloring of the first image?

Basically I took the below image and then post process it in gmsh such that i got the coloring of the first image.
First Image
mp
Second Image


But I what I want to do is to replicate the process that microstructpy do which is:
Starting from:

Then the above is converted to below:

And finally mesh the above into:

along with physical tags such that I can use that cell tags to define different material properties in different subdomains in fenics.

The coloring in gmsh, what does that correspond to? Vertex tags/Cell tags?
Without the msh file, it is really hard to grasp what information you have from the msh stand-point.

I used the below code to post process the image in gmsh and basically I varied the mesh intensity using the intensity of color
The .geo file which I executed in gmsh:

If(PostProcessing.NbViews == 0)
  // Merge the image (this will create a new post-processing view, View[0]), and
  // modify the normalized pixel values (v0) to make them reasonnable
  // characteristic lengths for the mesh
  Merge "mp.jpg";
  //Plugin(ModifyComponents).Expression0 = "v0 * 10";
  Plugin(ModifyComponents).Expression0 = "1 + v0^2 * 10";
  Plugin(ModifyComponents).Run;
EndIf

// Apply the view as the current background mesh
Background Mesh View[0];

// Build a simple geometry on top of the background mesh
w = View[0].MaxX;
h = View[0].MaxY;
Point(1)={0,0,0,w};
Point(2)={w,0,0,w};
Point(3)={w,h,0,w};
Point(4)={0,h,0,w};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line Loop(5) = {3,4,1,2};
Plane Surface(6) = {5};

outfile = StrCat(CurrentDirectory, "out.png");

DefineConstant[
  algo = {Mesh.Algorithm, AutoCheck 0, GmshOption "Mesh.Algorithm",
    Choices{1="MeshAdapt", 5="Delaunay", 6="Frontal", 8="DelQuad"},
    Name "Meshing parameters/Algorithm"},
  sizeFact = {Mesh.CharacteristicLengthFactor, AutoCheck 0,
    GmshOption "Mesh.CharacteristicLengthFactor", Min 0.1, Max 10, Step 0.1,
    Name "Meshing parameters/Element size factor"},
  sizeMin = {Mesh.CharacteristicLengthMin, AutoCheck 0,
    GmshOption "Mesh.CharacteristicLengthMin", Min w/100, Max w, Step 0.1,
    Name "Meshing parameters/Minimum element size"},
  save = {StrCat("View.ShowScale=0;Print '", CurrentDirectory, "out.png';"),
    AutoCheck 0, Macro "GmshParseString",
    Name "Save PNG"}
];

Mesh.ColorCarousel = 0;
Mesh.Color.Triangles = Black;
Mesh.Color.Quadrangles = Black;
Mesh.RecombineAll = (algo == 8);
Solver.AutoMesh = 2; // always remesh//+

The problem is that the msh file does not store any of the information regarding your image, and thus you have lost all the information you started with.

You could do something like:

with some slight modifications.

@dokken Thanks for the above reference, I followed the above method and was able to define different material properties, but how do i import the .msh file which I generated through gmsh along with the material information and visualize them together, basically I want to generate an .xdmf file which I can use in fenics also and in paraview also. Thanks for the first part again.

I presume you have read your mesh into DOLFIN. Then you can write the mesh to file (a xdmf-file), along with using XDMFFile.write_checkpoint(u, "u", 0.0) to write the data to file alongside the mesh.

That is my question that how do i import my another ‘.msh’ file (which i generated using gmsh which was giving the error of absence of physical tags earlier) into dolfin and using both the data from ‘.dta’ file about material property and mesh from my ‘.msh’ file, how do i combine both and generate a ‘.xdmf’ file, I apologize for the trivial nature of the question.

Just remove the cell_data in these, i.e.

from dolfin import *
import meshio
msh = meshio.read("/home/ratan/Documents/gmsh/mp1.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif cell.type == "tetra":
        tetra_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)])
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

I want to deal with 2D mesh so I used the code given below to import .msh file into dolfin but there is an error
Code:

msh = meshio.read("/mnt/d/Research Projects/FEM for microstructures/modified_2cmp1.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)


mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("dolfinmesh.pvd").write(mesh)

Error:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_142/763113767.py in <module>
     13     return out_mesh
     14 
---> 15 line_mesh = create_mesh(msh, "line", prune_z=True)
     16 meshio.write("mf.xdmf", line_mesh)
     17 

/tmp/ipykernel_142/763113767.py in create_mesh(mesh, cell_type, prune_z)
     10     cells = mesh.get_cells_type(cell_type)
     11     if prune_z:
---> 12         out_mesh.prune_z_0()
     13     return out_mesh
     14 

NameError: name 'out_mesh' is not defined

This function makes no sense. I have no clue where you have gotten this from. Please use

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
    return out_mesh

which is adapted from Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken
by removing the cell data calls

I have modified code to -

msh = meshio.read("/mnt/d/Research Projects/FEM for microstructures/modified_2cmp1.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
    return out_mesh

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mf.xdmf", line_mesh)

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)


mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")

but now the following error comes:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/tmp/ipykernel_142/2838469791.py in <module>
     14 
     15 line_mesh = create_mesh(msh, "line", prune_z=True)
---> 16 meshio.write("mf.xdmf", line_mesh)
     17 
     18 triangle_mesh = create_mesh(msh, "triangle", prune_z=True)

~/.local/lib/python3.10/site-packages/meshio/_helpers.py in write(filename, mesh, file_format, **kwargs)
    186 
    187     # Write
--> 188     return writer(filename, mesh, **kwargs)

~/.local/lib/python3.10/site-packages/meshio/xdmf/main.py in write(*args, **kwargs)
    544 
    545 def write(*args, **kwargs):
--> 546     XdmfWriter(*args, **kwargs)
    547 
    548 

~/.local/lib/python3.10/site-packages/meshio/xdmf/main.py in __init__(self, filename, mesh, data_format, compression, compression_opts)
    352         if data_format == "HDF":
    353             self.h5_filename = self.filename.with_suffix(".h5")
--> 354             self.h5_file = h5py.File(self.h5_filename, "w")
    355 
    356         xdmf_file = ET.Element("Xdmf", Version="3.0")

~/.local/lib/python3.10/site-packages/h5py/_hl/files.py in __init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
    556 
    557             with phil:
--> 558                 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
    559                                  locking, page_buf_size, min_meta_keep, min_raw_keep,
    560                                  alignment_threshold=alignment_threshold,

~/.local/lib/python3.10/site-packages/h5py/_hl/files.py in make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, locking, page_buf_size, min_meta_keep, min_raw_keep, alignment_threshold, alignment_interval, meta_block_size, **kwds)
    129         # we default to earliest
    130         low, high = h5f.LIBVER_EARLIEST, h5f.LIBVER_LATEST
--> 131     plist.set_libver_bounds(low, high)
    132     plist.set_alignment(alignment_threshold, alignment_interval)
    133 

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/h5p.pyx in h5py.h5p.PropFAID.set_libver_bounds()

ValueError: High bound is not valid (high bound is not valid)

Just remove

As it is not needed when you dont have boundary tags.

I removed the boundary tag part -

msh = meshio.read("/mnt/d/Research Projects/FEM for microstructures/modified_2cmp1.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
    return out_mesh

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)

Now error comes from the triangle one -

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/tmp/ipykernel_142/3638126479.py in <module>
     14 
     15 triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
---> 16 meshio.write("mesh.xdmf", triangle_mesh)
     17 
     18 mesh = Mesh()

~/.local/lib/python3.10/site-packages/meshio/_helpers.py in write(filename, mesh, file_format, **kwargs)
    186 
    187     # Write
--> 188     return writer(filename, mesh, **kwargs)

~/.local/lib/python3.10/site-packages/meshio/xdmf/main.py in write(*args, **kwargs)
    544 
    545 def write(*args, **kwargs):
--> 546     XdmfWriter(*args, **kwargs)
    547 
    548 

~/.local/lib/python3.10/site-packages/meshio/xdmf/main.py in __init__(self, filename, mesh, data_format, compression, compression_opts)
    352         if data_format == "HDF":
    353             self.h5_filename = self.filename.with_suffix(".h5")
--> 354             self.h5_file = h5py.File(self.h5_filename, "w")
    355 
    356         xdmf_file = ET.Element("Xdmf", Version="3.0")

~/.local/lib/python3.10/site-packages/h5py/_hl/files.py in __init__(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
    556 
    557             with phil:
--> 558                 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
    559                                  locking, page_buf_size, min_meta_keep, min_raw_keep,
    560                                  alignment_threshold=alignment_threshold,

~/.local/lib/python3.10/site-packages/h5py/_hl/files.py in make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, locking, page_buf_size, min_meta_keep, min_raw_keep, alignment_threshold, alignment_interval, meta_block_size, **kwds)
    129         # we default to earliest
    130         low, high = h5f.LIBVER_EARLIEST, h5f.LIBVER_LATEST
--> 131     plist.set_libver_bounds(low, high)
    132     plist.set_alignment(alignment_threshold, alignment_interval)
    133 

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/h5p.pyx in h5py.h5p.PropFAID.set_libver_bounds()

ValueError: High bound is not valid (high bound is not valid)

How did you install h5py?
See

I am getting the below message from terminal