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

[violetus]

@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

[dokken]

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?

[violetus]

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

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.

[dokken]

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.

[violetus]

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//+

[dokken]

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:

How to define a subdomain using material data (from data file)?

First of all, does your mesh consist of triangles or quadrilaterals? Secondly, is your mesh uniform and structured (i.e. is it a RectangleMesh with 100 x 100 cells?) I would suggest something like: from dolfin import * import numpy as np mesh = UnitSquareMesh(100, 100) V = FunctionSpace(mesh, "DG", 0) u = Function(V) dx = 1/100 values = np.arange(100*100) local_values = np.zeros_like(u.vector().get_local()) for cell in cells(mesh): midpoint = cell.midpoint().array() i = int(midpoi…

with some slight modifications.

[violetus]

@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.

[dokken]

@violetus

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.

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.

[violetus]

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.

[dokken]

@violetus

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)

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)

[violetus]

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

[dokken]

@violetus

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

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

[violetus]

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)

[dokken]

Just remove

@violetus

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

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

[violetus]

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)

[dokken]

How did you install h5py?
See

Keras model.load_weights() error: ValueError: Invalid high library version bound (invalid high library version bound)

[violetus]

I am getting the below message from terminal

[dokken]

You need to uninstall your current installation of h5py first

[violetus]

I uninstalled and reinstalled h5py using pip install --no-binary=h5py h5py but now this error comes :

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/tmp/ipykernel_457/3234862839.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")

AttributeError: module 'h5py' has no attribute 'File'

[violetus]

@dokken is the version of h5py causing this error ?