asign material to cells in a mesh toughio object

[geowomen]

Dear all, Dear Kerfonluu,

I want to assign a type of material to a complex geometry. I know the position (x, y, z) where I should assign materials.

Using the parameters from the toughio reader, I can define a function to do so:

def ajout_mat_box(my_box1, my_box2, parameters, name_material):
    for elem in parameters['elements'].keys():
        center = parameters['elements'][elem]['center']
        if inside_my_domain(my_box1, my_box2, center):
            parameters['elements'][elem]['material'] = name_material

where inside_my_domain is a function defining the domain (not a box).

But I want to use the mesh object constructed to select the element where to apply the domain:

dx = [length_along_x / 150] * 150
dy = [length_along_y / 10] * 10

dz = [length_along_z / 150] * 150
mesh = toughio.meshmaker.structured_grid(dx, dy, dz, origin=[-length_along_x/2, -length_along_y/2, -length_along_z/2], material="COOXX")
mesh.set_material("CONC_", mesh.filter.box(x0=-13.26, y0=-length_along_y/2, z0=-13.26, dx=13.26*2, dy=length_along_y, dz=13.26*2))
mesh.set_material("EDZIN", mesh.filter.box(x0=-8.32, y0=-length_along_y/2, z0=-8.32, dx=8.32*2.0, dy=length_along_y, dz=8.32*2.0))
mesh.set_material("EDZOT", mesh.filter.box(x0=-5.20, y0=-length_along_y/2, z0=-5.20, dx=5.20*2.0, dy=length_along_y, dz=5.20*2.0))
mesh.set_material("BFILL", mesh.filter.box(x0=-4.15, y0=-length_along_y/2, z0=-4.15, dx=4.15*2.0, dy=length_along_y, dz=4.15*2.0))
mesh.plot(show_edges=True)
mesh.write("meshhh.vtu")

If I don't have a simple box, how can I define a complex domain filter by looping in cells and select only the cells in my complex domain?

Thank you for your help.

[keurfonluu]

Hi @geowomen,

There is unfortunately no obvious way to do that in 3D, this is usually handled by the meshing software when setting up the geometry.
In 2D, you can use shapely to determine whether an element center is inside any polygon.

From the name of the materials, I guess that you are trying to model tunnels. Is your geometry symmetric along the tunnel axis? That would be helpful as you could generate a 2D mesh, assign the materials using polygons, and then extrude the mesh along the tunnel axis.

[geowomen]

Thank you Kerfonluu for your answer. If I understand correctly, there is no way to access the center of cells from a toughio mesh structure and create a loop to select specific positions and assign specific materials?

As for my second question, is there a way to extrude a surface or geometry with a line, the line is changing direction, not just a constant direction?

[keurfonluu]

You can access cell centers using mesh.centers.
You can only extrude in one direction with mesh.extrude_to_3d, but you can easily shift the coordinates of the points afterwards (i.e., mesh.points) following your line.

[geowomen]

import meshio
dx = [longeur_suivant_x/150]*150
dy = [longeur_suivant_y/10]*10
dy1 = np.linspace(-longeur_suivant_y/2,-longeur_suivant_y/2 +40,2)
dy2 = np.linspace(-longeur_suivant_y/2 +40,-longeur_suivant_y/2 +40 +13.26,15)
dy3 = np.linspace(-longeur_suivant_y/2 +40 +13.26,longeur_suivant_y/2,15)
dz = [longeur_suivant_z/150]*150
mesh = toughio.meshmaker.structured_grid(dx, dy, dz, origin = [-longeur_suivant_x/2,-longeur_suivant_y/2,-longeur_suivant_z/2],material="COOXX")
def insid_box(my_box,centre):
    x_good = 0
    y_good = 0
    if ( centre[0] < my_box['x_max']  and  centre[0] > my_box['x_min']):
        x_good = 1.0
    if ( centre[2] < my_box['z_max']  and  centre[2] > my_box['z_min']):
        y_good = 1.0
    #print("my center =", centre, " give me ", "x_good =", x_good ,"y_good =", y_good)
    return x_good and y_good
def inside_two_box(my_box1,my_box2,blk):
    #print(" insid_box(my_box1,blk) =",insid_box(my_box1,blk))
    #print(" insid_box(my_box2,blk) =",insid_box(my_box2,blk))
    if insid_box(my_box1,blk) and (not insid_box(my_box2,blk)):
        return True
    return False

def ajout_mat_box(my_box1,my_box2,parameters, name_materiel):
    for  elem in parameters['elements'].keys() :
        centre = parameters['elements'][elem]['center']
        if inside_two_box(my_box1,my_box2,centre):
            parameters['elements'][elem]['material'] = name_materiel
#mesh.plot(show_edges=True)
mesh.write("mesh_travail.vtu")
my_box1_remblais = {'x_max' :4.15  , 'x_min':-4.15  , 'z_max' :4.15 , 'z_min':-4.15 }
my_box2_remblais = {'x_max' : 0 , 'x_min':0 , 'z_max' : 0, 'z_min': 0}
# Load your mesh from the .vtu file
mesh = meshio.read("mesh_travail.vtu")
dir(mesh.cells)
for cell in mesh.cells:
    print(cell)
print("matt =",mesh.cell_data['material'])
# Example: Assign material "BET" to cells where x > 2 and z > 1
print("mesh.cells[0].data =",mesh.cells[0].data)
for cell_id, cell in enumerate(mesh.cells[0].data):
    # Get the vertices (nodes) of the cell
    vertices = mesh.points[cell.data]
    if cell_id < 2 : print((vertices.shape))
    if cell_id < 2 :print("vertices[cell_id]",vertices[7])
    if cell_id < 2 :print("vertices",vertices)
    # Check if any vertex satisfies the conditions
    #if any(vertices[cell_id][:][0] > 2):
    #print([inside_two_box(my_box2_remblais,my_box1_remblais,vertices[k])  for k in range(8)])
    if np.all([inside_two_box(my_box1_remblais,my_box2_remblais,vertices[k])  for k in range(8)]):
        # Assign material ID "BET" to this cell
        print("yes i found someting")
        print(mesh.cell_data["material"][0][cell_id])
        mesh.cell_data["material"][0][cell_id] = np.full_like(mesh.cell_data["material"][0][cell_id],2)

my_mesh = toughio.from_meshio(mesh, material='dfalt')
# Write the modified mesh to a new file
mesh.write("modified_mesh0.vtu", file_format="vtu")
mesh.write("modified_mesh.vtu", file_format="vtu")