Defining subdomains for discontinuous composite material

I have a composite material having sphere-shaped particles as reinforcements. So, there will be two different materials, one for matrix and one for these spheres. How can I define the boundaries of the subdomain for these spheres? I’ve tried to get cell coordinates by cell.get_coordinates_dofs, but I think it didn’t work (maybe I’m missing something).

I’m using pygmsh to create geometry and mesh:

height = 0.020
radius = 0.020
with pygmsh.occ.Geometry() as geom:

        geom.characteristic_length_max=0.003  # smaller = more mesh elements
        cyl = geom.add_cylinder([0.0, 0.0, 0.0], [0.0, 0.0, height], radius, mesh_size=1)
        sphere_list=[]

        # Place voids with 3mm radius in BCC configuration
        for i in np.linspace(-0.015, 0.015, 4):
          for j in np.linspace(-0.015, 0.015, 4):
            for k in np.linspace(-0.005, 0.015, 3):
              sphere_list.append(geom.add_ball([i, j, k], 0.003))

        spheres=geom.boolean_union(sphere_list)
        cyl==geom.boolean_difference(cyl,spheres)
        mesh = geom.generate_mesh()

See for instance: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

I’ve tried to implement you solution but I couldn’t define subdomain for the spheres. I mean, how should I enter the coordinates of all spheres in the inside function?

Make several functions that each define one sphere, and tag them with a marker one by one. This can br scripted in Python. You could also use the Gmsh Python api to directly define the subdomains of yhr mesh, as shown in: Electromagnetics example — FEniCSx tutorial