Definition of Subdomain by Element Sets

Hello everyone,

I checked all of the previous related posts but I could not find a way to define subdomain by element sets (element IDs). Here , I have an Abaqus .inp mesh file. I am able to convert it to .xdmf file by meshio as following:

import meshio
import numpy as np
mesh = meshio.read("trial_subdomain.inp")
cells = mesh.get_cells_type("tetra")
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 == "tetra":
            # Create marker with int from integer from name of input
            cell_data[entities] = int(marker[-1])

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

My sample model contains 2 of the subdomain as shown below:

Now, I want to define subdomain for outer and inner regions. I have element IDs as following:

OUTER_ELEMENT ELEM
1 2 3 4 5 6 7 8 9 10
11 12 13 14 15 16 17 18 19 20
21 22 23 24 25 26 27 28 29 30
31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 46 47 48 49 50
51 52 53 54 55 56 57 58 59 60
61 62 63 64 65 66 67 68 69 70
71 72 73 74 75 76 77 78 79 80
81 82 83 84 85 86 87 88 89 90
91 92 93 94 95 96 97 98 99 100
101 102 103 104 105 106 107 108 109 110
111 112 113 114 115 116 117 118 119 120
121 122 123 124 125 126 127 128 129 130
131 132 133 134 135 136 137 138 139 140
141 142 143 144 145 146 147 148 149 150
151 152 153 154 155 156 157 158 159 160
161 162 163 164 165 166 167 168 169 170
171 172 173 174 175 176 177 178 179 180
181 182 183 184 185 186 187 188 189 190
191 192 193 194 195 196 197 198 199 200
201 202 203 204 205 206 207 208 209 210
211 212 213 214 215 216 217 218 219 220
221 222 223 224 225 226 227 228 229 230
231 232 233 234 235 236 237 238 239 240
241 242 243 244 245 246 247 248 249 250
251 252 253 254 255 256 257 258 259 260
261 262 263 264 265 266 267 268 269 270
271 272 273 274 275 276

AND:

INNER_ELEMENT ELEM
277 278 279 280 281 282 283 284 285 286
287 288 289

I thought that, using element IDs is better way to describe subdomains rather than node IDs. Because if I change the order of element then nodal selection will fail in the subdomain. So, how can i define two of the subdomain by using element IDs ? Which formatting should I follow in txt. or something ? Thanks

Regards,

I would suggest changing the corresponding entries of cell data to have a unique marker for this list of elements.
As you have not provided a plain-text version of the mesh, and i am not used to work with inp, i cannot give you any code example or test if this would resolve your issue.

Hello dokken,

Thank you for your response. Actually I do not have any definition of subdomain or lets say element sets are not located in .inp file. I extracted Element IDs manually from another software as above. I am trying to create the subdomains by these elements. So, you can assume that, I have a mesh (.xdmf, .xml, .inp, etc.) and I am trying to define subdomains by element numbers (which is not in the mesh file).
And just for your information, this is how .inp looks like:

** MESH NODES
**
** Number of nodes 85
*NODE
1, 3.000000000E+001, 0.000000000E+000, 3.000000000E+001
2, 0.000000000E+000, 0.000000000E+000, 3.000000000E+001
3, 2.000000000E+001, 0.000000000E+000, 3.000000000E+001
4, 1.000000000E+001, 0.000000000E+000, 3.000000000E+001
5, 3.000000000E+001, 3.000000000E+001, 3.000000000E+001
6, 3.000000000E+001, 2.000000000E+001, 3.000000000E+001
7, 3.000000000E+001, 1.000000000E+001, 3.000000000E+001
8, 0.000000000E+000, 3.000000000E+001, 3.000000000E+001
9, 1.000000000E+001, 3.000000000E+001, 3.000000000E+001
10, 2.000000000E+001, 3.000000000E+001, 3.000000000E+001
11, 0.000000000E+000, 1.000000000E+001, 3.000000000E+001
12, 0.000000000E+000, 2.000000000E+001, 3.000000000E+001
13, 3.000000000E+001, 0.000000000E+000, 0.000000000E+000
14, 3.000000000E+001, 0.000000000E+000, 1.000000000E+001
15, 3.000000000E+001, 0.000000000E+000, 2.000000000E+001
16, 0.000000000E+000, 0.000000000E+000, 0.000000000E+000
.
.
.
**
** MESH ELEMENTS
**
*ELEMENT, TYPE=C3D4
1, 38, 68, 72, 39
2, 48, 51, 72, 49
3, 33, 41, 70, 45
4, 6, 70, 45, 7
5, 38, 39, 72, 53
6, 7, 45, 47, 70
7, 33, 70, 40, 43
8, 9, 67, 52, 10
9, 26, 29, 32, 75
10, 1, 15, 3, 77
11, 40, 64, 67, 65
12, 2, 78, 11, 66
13, 12, 71, 59, 23
.
.
.
.

The point is that cell_data is an array, with one entry per row in cells. Therefore, you can create an XDMFFile with this cell data by replacing the cell_data-indices of the cells that you have in outer_element and inner_element by new marker values. Let us consider a simple 2D example:

import meshio
import numpy as np
x = np.array([[0,0],[1,0],[0,1],[1,1], [0,2]])
cells = np.array([[0,1,2],[1,2,3], [2,3,4]], dtype=np.int32)
cell_data = np.array([33, 55,11],dtype=np.int32)
mesh = meshio.Mesh(points=x, cells={"triangle":cells}, cell_data={"name_to_read":[cell_data]})
meshio.write("test_mesh.xdmf", mesh)

In this code, you have three cells, where the first cell has the marker 33, second 55 and third 11. Say you want to change the marker value of the last two (i.e. index 1,2). You can then do:

cell_data[np.array([1,2],dtype=np.int32)] = 3

This will cause cell 1 and 3 to now have the value 3. The cell values can then be loaded into dolfin through a mesh value collection (as described in many other posts) and used to define subdomains.
I would suggest reading in the cell regions as numpy arrays when you convert the mesh, and replace values in cell data.

Hello dokken,

Thank you for your guidance. I finally got the cell_data as following. I think it represents the “0” and “1” marker tags for element sets now.

array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
        1.,  1.,  1.])

Now, I am trying to figure out how to use the mesh value connection function. I checked other posts but they used a separate .xdmf files for subdomains and they import it as

with XDMFFile(mpi_comm, boundaries_filename) as xdmf_infile:
    xdmf_infile.read(mvc_boundaries)

But I need to use the information from cell_data and I dont have an additional .xdmf file for domains. So, I could not go further to define subdomains by this procedure. Thank you.

Regards

To continue illustrating this with the example above, you can do the following:

import meshio
import numpy as np
x = np.array([[0,0],[1,0],[0,1],[1,1], [0,2]])
cells = np.array([[0,1,2],[1,2,3], [2,3,4]], dtype=np.int32)
cell_data = np.array([33, 55,11],dtype=np.int32)
cell_data[np.array([1,2],dtype=np.int32)] = 3
mesh_out = meshio.Mesh(points=x, cells={"triangle":cells}, cell_data={"name_to_read":[cell_data]})
meshio.write("test_mesh.xdmf", mesh_out)


from dolfin import *
mesh = Mesh()
with XDMFFile("test_mesh.xdmf") as xdmf:
    xdmf.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
    xdmf.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

dx = Measure("dx", domain=mesh, subdomain_data=mf)
print(assemble(1*dx(33)), assemble(1*dx(3)))

as you observe from the output, the mesh function contains each of the subdomains.

1 Like

Hello dokken,

Thank you a lot, it works well. As a last question; in this model I do not need external boundary definitions on mesh but probably I will need it in the next studies. Is it defined by a similar way ? In our current example, we define the markers by cell_data array, what will be used for boundary definitions from external software data ? Is there any example about it ? To define boundaries, elements are not an option and nodes are not enough I think, so is there any procedure for that also ? Thanks.

There are plenty of other topics in the forum discussing how to import boundary markers. See for instance: Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #3 by dokken
Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken
Converting simple 2D mesh from gmsh to dolfin - #4 by dokken
Create mesh with gmsh, convert it and import it in dolfin, problems - #2 by dokken

1 Like

Thank you a lot , I will follow these topics.