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

[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 ?

[violetus]

Thanks @dokken for the reference, the error was solved when i used sudo apt-get install python3-h5py after uninstalling h5py. Can you please tell me how can I write my .dat file information into this XDMF file.
Thank You

[dokken]

See:

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…

Which writes the data to a function, Which you can write to xdmf

[violetus]

Thank you very much @dokken for your help, I was able to figure out how I will assign material properties because of you, just one last doubt is that I want to define boundaries too like here in the image shown below.

But unfortunately i dont have facet_region.xml file and want to do it with .xdmf, can you please tell me how I can do that ?
Thank You for the help.

[dokken]

The facet_region.xml file would extract facet tags from your mesh. As your mesh file does not have any, you dont have anything to extract. With xdmf you would create two files, one for the mesh and one for the facets.
This has for instance been discussed in detail in Gmsh, Meshio, Dolfin

@violetus

Thank you very much @dokken for your help, I was able to figure out how I will assign material properties because of you,

I would strongly suggest you post your solution on the forum, as it might help others in the future.

[violetus]

@dokken Yeah sure, I will share the code here but I am not getting the results I wanted (please see the image attached) so kindly help that in which step I making mistake, here’s the code:

materials = MeshFunction('double', mesh, 2)

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)

values = np.arange(4383*1) 
local_values = np.zeros_like(u.vector().get_local())
local_values_material = np.zeros_like(u.vector().get_local())

for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0])
    j = int(midpoint[1])
    k = int(midpoint[2])
    ## partha idea
    local_values_material[cell.index()] = material_array_final[cell.index()]
    materials[cell] = int(local_values_material[cell.index()])
    print(midpoint, i, j , k, "Material:", materials[cell], "Cell index: ", cell.index())
    
u.vector().set_local(local_values)
dolfin.XDMFFile(dolfin.MPI.comm_world, "mesh1.xdmf").write_checkpoint(u,"u",0)

class al(UserExpression):
    def __init__(self, materials, al0, al1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.k_0 = al0
        self.k_1 = al1
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 1:
            values[0] = self.k_0
        else:
            values[0] = self.k_1
            
E1 = 210e3
nu1 = 0.3
E2 = 21e3
nu2 = 0.25

E = al(materials, E1, E2, degree = 0)
nu = al(materials, nu1, nu2, degree = 0)

File('e.pvd') << project(E, FunctionSpace(mesh, 'DG', 0))
File('materials.pvd') << materials

I wanted the material properties as shown in gmsh but I got the right one as shown in paraview. I’m not able to figure it out that when I have assigned each cell with a cell tag of either 0 or 1 and calling the same in user-defined function then where am I doing mistake.
Second Question: I tried making another xdmf file as @dokken suggested for facets by referring to this post, but was not able to make it.
I tried:

import meshio

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

meshio.write("mesh.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
    )
)

meshio.write("mf.xdmf", meshio.Mesh(points=msh.points[:,:2], cells={"line": msh.cells_dict["line"]}))

from dolfin import * 

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)

plot(mesh)

got the error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_79/2720985724.py in <module>
     19 mvc = MeshValueCollection("size_t", mesh, 1)
     20 with XDMFFile("mf.xdmf") as infile:
---> 21     infile.read(mvc, "name_to_read")
     22 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
     23 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Any input will be appreciated.
Thank You

[howstheJosh]

In the following line:

infile.read(mvc, "name_to_read")

The "name_to_read" should be the marker "line" which we used in meshio.write("mf.xdmf", meshio.Mesh(points=msh.points[:,:2], cells={"line": msh.cells_dict["line"]}))

This is because we use that marker when writing that dictionary to the xdmf file. Here we use the "line" marker since we usually set physical tags or groups on the lines (facets in 2d) on which we will later set the boundary conditions. You can refer to Gmsh, Meshio, Dolfin for further clarification.

[violetus]

@howstheJosh Thanks for the reply, I referred to Gmsh, Meshio, Dolfin suggested by @dokken and therefore wrote the condition for facets, I tried replacing what you suggested but still the error persists.
Modified Code:

import meshio

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

meshio.write("mesh.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
    )
)

meshio.write("mf.xdmf", meshio.Mesh(points=msh.points[:,:2], cells={"line": msh.cells_dict["line"]}))

from dolfin import * 

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, "line")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

plot(mesh)

Error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_28/1003279362.py in <module>
     19 mvc = MeshValueCollection("size_t", mesh, 1)
     20 with XDMFFile("mf.xdmf") as infile:
---> 21     infile.read(mvc, "line")
     22 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
     23 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Can you please tell me what can be the possible source for this error which I could change ?
Thank You
PS: Since I am only using a plane surface and meshed it and via post processing I added image on this plane surface (kindly see the above post from where I have posed the question), I am just using mesh without any info stored in it about facet_tags or different material properties, I did that using defining different material properties at the midpoint of the cell via a .dat file which contains each point material propery, but the results are not showing in the paraview.

[howstheJosh]

Oh that was my bad. I had not read the previous answers. As @dokken mentioned in a previous reply, since you don’t have any facet regions with a physical tag:

mvc = MeshValueCollection("size_t", mesh, 1) 
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "line")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

This part of the code would not be required as it is basically reading the tagged boundaries into the object mf using the MeshValueCollection and MeshFunction functions.
Please remove these lines and try again.

[violetus]

Yes, it worked. Thank You for this. Now my main query remains that how do I replicate the inclusions ? (kindly see the image below) basically I will give an overview of my problem.
I wanted to define different material properties to red and blue part of my image so I post process the image in gmsh and meshed it but my material properties were not saved into .msh file so I took the image and extracted the points in .dat file as (x, y, z, material identity) and assigned them different material properties ( 0 → blue, 1 → red) and then I imported .msh file and .dat file and took the cells from the msh file and then assigned cell midpoint with either 0 or 1 by comparing it with dat file coordinates and then ran the code but still got very different result from what I wanted (kindly see the image below).

I humbly request to the community to kindly give me some suggestions on how to resolve this error.
Thank You

[dokken]

What are the xyz coordinates in the Dat file? Do They correspond to vertex values or cell midpoints?

As you have not given a minimal reproducible example, say a Dat file for a 2x2 unit square mesh, I cannot help you any further.

[violetus]

just a quick question @dokken I have posed a question regarding installation of mgis.fenics python module in this post but I didn’t get any suggestion or guidance on what to do. I also followed this post and tried all the steps in it but still it doesn’t work. I have also mailed on contact mail id given on their official site today, can you please tell me where can I get this issue resolved as soon as possible ?
Thank You

[dokken]

I have never used mgis.fenics, so I cannot comment on who should be contacted regarding this, except the people that made the software, Ref thelfer (Thomas Helfer) · GitHub