How to use gmsh tags to assign material properties?

mesh
1

Hello, I have a square plate which has 3 physical entities, namely boundary, inclusions, and domain. I wish to assign different material properties to the domain and inclusions. The inclusions are circular and randomly placed, so I am not aware of the coordinates. How do I do this in FEniCS using the data I already have from gmsh?

import meshio
import dolfin
from dolfin import XDMFFile
import numpy as np
from dolfin import *
import matplotlib

msh = meshio.read("/mnt/c/Users/ABC/sfepy/testrve.msh")
print(msh)

The output of the above print function is the no. of points, lines and cells, and finally the following

Cell sets: boundary, domain, inclusions, gmsh:bounding_entities
Point data: gmsh:dim_tags
Cell data: gmsh:physical, gmsh:geometrical
Field data: boundary, domain, inclusions

The further I have written the mesh into the xdmf files as

for cell in msh.cells:
    if cell.type == "line":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

def create_mesh(mesh, cell_type, prune_z=True):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

line_mesh = create_mesh(msh, "line", prune_z=False)
meshio.write("mf.xdmf", line_mesh)
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
meshio.write("mesh.xdmf", triangle_mesh)

from mpi4py import MPI as mpi
from scipy.spatial.distance import cdist
import numpy as np

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
#boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dx = Measure("dx", domain=mesh, subdomain_data=mf_2d)

mesh
mesh817×786 42.7 KB

This is what my mesh looks like

EDIT

I tried to implement what was given in the following link,
How to extract and use physical labels from GMSH surfaces on FEniCS (2D mesh)

from mpi4py import MPI as mpi
from scipy.spatial.distance import cdist
import numpy as np

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf_2d)
dx1 = Measure("dx", domain=mesh, subdomain_data=mvc_2d,   subdomain_id=1)
dx2 = Measure("dx", domain=mesh, subdomain_data=mvc_2d, subdomain_id=2)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))

The result of which is

0.9999999999999996
0.0
0.0

2

What does

return?

I would strongly suggest you have a closer look at the post that you linked to above,

and try to figure out what is different in the mesh.xdmf file you have above and the one in the other question.

3 
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf_2d)
dx1 = Measure("ds", domain=mesh, subdomain_data=mvc_2d,   subdomain_id=1)
dx2 = Measure("ds", domain=mesh, subdomain_data=mvc_2d, subdomain_id=2)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))

returns

0.9999999999999964
0.0
0.0

I also tried something else, that is running dx(3),

domains = MeshFunction("size_t", mesh, mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mf_2d)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(3)))

which returns

0.9999999999999964
0.03609167866767062

I went through the gmsh script posted in the link, in gmsh it appears as

forumrve
forumrve843×823 12 KB

Here the square is completely filled with 4 different physical groups, while in my file, the square itself is one physical group, and the circular inclusions are another physical group, which are randomly placed.

4

The point here is what are the actual values you have tagged each circle with (look at your mesh.xdmf file in paraview).

Secondly, the marker you should be using for dx integrals should originate from mesh.xdmf, not

as this would read in the marked facets

5

meshview
meshview671×550 270 KB

The mesh looks like so in Paraview, with the parameter ‘name_to_read’ going from 1 to 2 (Square is 1, Circles are 2). So I do think that the xdmf file is able to make the distinction.

With regards to the 2nd part, I’m unable to understand what you are trying to say

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")
mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
print(mesh.topology().dim())
dx = Measure("dx",domain=mesh, subdomain_data=mvc_2d)
dx1 = Measure("dx", domain=mesh, subdomain_data=mvc_2d)
dx2 = Measure("dx", domain=mesh, subdomain_data=mvc_2d)
print(assemble(Constant(1)*dx))
print(assemble(Constant(1)*dx(1)))
print(assemble(Constant(1)*dx(2)))

Are you saying that I replace mf_2d with mvc_2d?

6

dx is a volume integral measure. The corresponding subdomain_data should use volume_markers, those marked inside mesh.xdmf, not facet markers (which are the once you are reading into your mesh value collection).

7

I’m sorry but I’m still not sure I understand.
You are saying that in
dx = Measure("dx",domain=mesh, subdomain_data=mf_2d)
mf_2d should be replaced with something else, which are supposed to be volume_markers . Where do I find these markers?

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())

I guess it should be similar to this?

Is there any other link or documentation I should refer?
Also, if my domain is 2D, why should I use volume markers?
Apologies in advance for my lack of comprehension.

8

When you call this, you add cell_data for the triangular (volume in 2D) elements.

These are the once you want to read in:

[quote=“Sarvesh_Sontakke, post:1, topic:9249”]

mvc_2d = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    infile.read(mvc_2d, "name_to_read")
cell_markers = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)