Subdomain definition based on input file

mesh
1

Hello everybody,

I am a new user to fenics and have run into an issue that I can not resolve. I hope you can help/guide me.

I am trying to run a simple diffusion analysis on a 2D domain. The domain is obtained from image processing of microstructural images of steel. Hence the domain is split into different regions, where the diffusion coefficients are different. Therefore I need to import the mesh definition, as well as the subdomain definition. I have already looked into forum posts and discussions on how to do this for domains, where subdomains are defined by functions. But in my problem, the subdomain definition is random and hence I do not know how to define my subdomains.

I have stored my mesh in abaqus format (.inp). This is just for me to view the geometry in abaqus. Basically I have the node definitions, element connectivity matrix and information on elements belonging to different subdomains stored as feature vectors in the input file. I have generated this input file, so I have access to all of the nodes, elements and features matrices. All the elements are triangles.

So I have the following information(this is just an example to show you the information I have):

nodes =
[Node no., x, y, z
1, 0.0, 0.0, 0.0;
2, 1.0, 0.0, 0.0; …]
elements =
[Ele no., n1, n2, n3
1, 1, 2, 3;
2, 1, 2, 4; …]
subdomain1 = [1, 2, 3, …] (list of elements belonging to subdomain1)
subdomain2 = [4, 5, 6, …]
diffusion coefficients: D1, D2,…

I have already tried using meshio to convert the .inp file to .xdmf file. It worked and I was able to import the geometry but lost all information on the subdomain definitions. Any help is appreciated.

Mesh
Fenics file

Thank you in advance.

1 Like
2

You need to add the subdomains to the mesh (with the cell_data kwarg in meshio) to keep the data (assuming that meshio can read in the subdomain markers. You should be able to follow a strategy similar to: Converting simple 2D mesh from gmsh to dolfin - #4 by dokken
Note that I would strongly advice you to use xdmf and not xml formatting when converting your mesh.

1 Like
3

Thanks Dokken for the quick reply. I will have a look at how to implement what you have suggested. I will come back if I have any more questions.

4

Hello Dokken,

I tried using the meshio.write() as you suggested and tried the code in your comment. I ran the following line from your old comment with the mesh you have given:
meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}))
I got an error saying “TypeError: list indices must be integers or slices, not str”. I read in another forum that meshio changed the inputs for the write function from library to list. So I changed the cells to a list and was successful in writing the xdmf file without cell_data. I am not sure how I should define the cell_data. Could you please show me how I would define the cell_data.

I am using the following code to understand how the write function works.

import numpy as np
import meshio

points = [ [0, 0, 0], [0, 1, 0], [1, 0, 0], [1, 1, 0], ]
ele =[ [0, 1, 2], [1, 2, 3], ]
data = [[1], [2],]
cells = [("triangle", ele)]
cell_data = [("phi", data)]

meshio.write("check_mesh.xdmf", meshio.Mesh(points, cells, cell_data))

5

How about:

import meshio
import numpy as np
mesh = meshio.read("mesh.inp")
cells = mesh.get_cells_type("triangle")
cell_sets = mesh.cell_sets_dict
cell_data = np.zeros(len(cells))
for marker, set in cell_sets.items():
    for type, entities in set.items():
        if type == "triangle":
            # Create marker with int from integer from name of input
            cell_data[entities] = int(marker[-1])

out_mesh = meshio.Mesh(points=mesh.points, cells={
                       "triangle": cells}, cell_data={"name_to_read": [cell_data]})
meshio.write("mesh.xdmf", out_mesh)

It seems to me that you have not marked every cell in the mesh (look at the top right corner of the mesh), as well as the material parameters eems a bit oddly distributed.

6

Thank you Dokken. I really appreciate you helping me.
Yes, you are right. I just discovered a bug in my .inp file generation code a few hours ago. I have since then fixed it. Now all elements belong to a subdomain. I ran your code on the fixed input file and the generated xdmf file looks good. I will look at your code and try to understand it.
Thanks a lot again.