Hi Leo, I’m trying to define subdomains for a problem with multiple materials. I was able to do this nicely in 2D by defining the subdomains with the constructive solid geometry functions of the mshr module of FEniCS and with the domain.set_subdomain() method. But it turns out that this is not implemented for 3D so I am going to do this manually by iterating over all cells and checking if the cell is inside a specific geometrical volume, for example:
My subdomain is:
marrow = Cylinder(df.Point(11.0, 13.0, 0.0), df.Point(11.0, 13.0, 150.0), r_marrow, r_marrow, 15)
then I define a markers MeshFunction:
markers = MeshFunction('size_t', mesh, 3)
finally I check if a given cell is inside my subdomain and I mark that cell
for cell in cells(mesh):
# Find cell center
coord_celda = cell.get_vertex_coordinates()
center_celda_x = (coord_celda[0] + coord_celda[3] + coord_celda[6] + coord_celda[9])/4
center_celda_y = (coord_celda[1] + coord_celda[4] + coord_celda[7] + coord_celda[10])/4
center_celda_z = (coord_celda[2] + coord_celda[5] + coord_celda[8] + coord_celda[11])/4
# If cell is in marrow
if ((centro_celda_x - 11)**2 + (centro_celda_y - 13)**2 <= r_marrow**2):
# It is marked as subdomain 1
markers[cell.index] = 1