Importing PhysicalGroup of line elements (GMSH) using Mesh.load

[HoseongJeong]

Dear scikit-fem users and staff,

I would like to extend my appreciation to Kinnala and the community for your continuous support.

Recently, I have been working on modeling composite materials that include both concrete (volume) and rebar (line). While I successfully created the mesh and defined the PhysicalGroup for both concrete and rebar using GMSH, I was unable to find the PhysicalGroup for the rebar (line) in the mesh._subdomain. (I can find PhysicalGroup of volume in _subdomain)

Could anyone provide guidance on how to import the PhysicalGroup for line elements?

Thank you in advance!

-------------Below is my mesh generator--------------------

import gmsh
import numpy as np
import math
import copy
import sys
import pyvista as pv

gmsh.initialize()
gmsh.model.add("DFG 3D")

Longitudinal_rebar_dimension="3D"
Transverse_rebar_dimension="1D"
L = 500
B=200
H=200
d_l=8
c=20
l_b=70
d_st=8
c_st=20
s_st=200
offset_distance=1
radius_of_fillet=3*d_st

concrete = gmsh.model.occ.addBox(0, 0, 0, L, B, H)

if Longitudinal_rebar_dimension=="3D":
    Longitudinal_rebar = gmsh.model.occ.addCylinder(L/2-l_b/2, B/2, c+d_l/2, l_b, 0, 0, d_l/2)
    gmsh.model.occ.fragment([(3, concrete)], [(3, Longitudinal_rebar)])
    gmsh.model.occ.synchronize()
    concrete = gmsh.model.getEntities(dim=3)[1][1]
    Longitudinal_rebar_unbonded1 = gmsh.model.occ.addCylinder(0, B/2, c+d_l/2, L/2-l_b/2, 0, 0, d_l/2+offset_distance)
    Longitudinal_rebar_unbonded2 = gmsh.model.occ.addCylinder(L/2+l_b/2, B/2, c+d_l/2, L/2+l_b/2, 0, 0, d_l/2+offset_distance)
    elimination_list=[(3, Longitudinal_rebar_unbonded1),(3, Longitudinal_rebar_unbonded2)]
    gmsh.model.occ.cut([(3, concrete)], elimination_list)
    gmsh.model.occ.synchronize()

def Calculate_transverse_rebar_positions(L, s_st):
    if L%s_st==0:
        The_number_of_transverse_rebar=math.floor(L/s_st)-1
    else:
        The_number_of_transverse_rebar=math.floor(L/s_st)
    remaining_length = L - (The_number_of_transverse_rebar * s_st)
    cover_distance = remaining_length / 2
    positions = [cover_distance + i * s_st for i in range(The_number_of_transverse_rebar + 1)]
    return positions

position=Calculate_transverse_rebar_positions(L, s_st)
point_list=[]
point_list_each_Transverse_rebar=[]

if Transverse_rebar_dimension=="1D":
    original_point=[[(c_st-d_st/2-offset_distance), (c_st-d_st/2-offset_distance)],
                    [B-(c_st-d_st/2-offset_distance), (c_st-d_st/2-offset_distance)],
                    [B-(c_st-d_st/2-offset_distance), H-(c_st-d_st/2-offset_distance)],
                    [(c_st-d_st/2-offset_distance), H-(c_st-d_st/2-offset_distance)]
                    ]
    tag_count=1000
    corners=["Bottom left","Bottom right","Top right","Top left"]
    
    for i in position:
        for j in range(len(corners)):
            gmsh.model.occ.addPoint(i, original_point[j][0], original_point[j][1],tag=tag_count)
            point_list_each_Transverse_rebar.append(tag_count)
            tag_count=tag_count+1       
        point_list.append(point_list_each_Transverse_rebar)
        point_list_each_Transverse_rebar=[]
            
    Transverse_rebar_list=[]
    loop_list=[]
    
    for i in range(len(position)):
        line1=gmsh.model.occ.addLine(point_list[i][0],point_list[i][1])
        Transverse_rebar=(1,line1)
        Transverse_rebar_list.append(Transverse_rebar)
        
        line2=gmsh.model.occ.addLine(point_list[i][1],point_list[i][2])
        Transverse_rebar=(1,line2)
        Transverse_rebar_list.append(Transverse_rebar)
        
        line3=gmsh.model.occ.addLine(point_list[i][2],point_list[i][3])
        Transverse_rebar=(1,line3)
        Transverse_rebar_list.append(Transverse_rebar)
        
        line4=gmsh.model.occ.addLine(point_list[i][3],point_list[i][0])
        Transverse_rebar=(1,line4)
        Transverse_rebar_list.append(Transverse_rebar)
        
        gmsh.model.occ.synchronize()
        
    for Transverse_rebar in Transverse_rebar_list:
        gmsh.model.mesh.embed(1, [Transverse_rebar[1]], 3, 3)
    
    gmsh.model.occ.removeAllDuplicates()
    volumes = gmsh.model.getEntities(dim=3)
    
    concrete_marker = 1
    gmsh.model.addPhysicalGroup(volumes[1][0], [volumes[1][1]], concrete_marker)
    gmsh.model.setPhysicalName(volumes[1][0], concrete_marker, "Concrete volume")
    
    Longitudinal_rebar_marker = 2
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], Longitudinal_rebar_marker)
    gmsh.model.setPhysicalName(volumes[0][0], Longitudinal_rebar_marker, "Longitudinal_rebar volume")
    
    Transverse_rebar_marker = 5000
    Transverse_rebar_dimension_list=[i[0] for i in Transverse_rebar_list]
    Transverse_rebar_tag_list=[i[1] for i in Transverse_rebar_list]
    gmsh.model.addPhysicalGroup(1, Transverse_rebar_tag_list, Transverse_rebar_marker)
    gmsh.model.setPhysicalName(1, Transverse_rebar_marker, "Transverse_rebar line")

    resolution = d_l/2
    Longitudinal_rebar_surface=gmsh.model.getBoundary([(3, Longitudinal_rebar_marker)], oriented=False)
    boundary=Longitudinal_rebar_surface
    
    for i in range(len(boundary)):
        boundary[i]=boundary[i][1]
    
    Longitudinal_rebar_surface_marker=1000
    gmsh.model.addPhysicalGroup(2, boundary, Longitudinal_rebar_surface_marker)
    gmsh.model.setPhysicalName(2, Longitudinal_rebar_surface_marker, "Longitudinal_rebar surface")
    
    gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 6)
    gmsh.option.setNumber("Mesh.MeshSizeMin", 1* resolution)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 10 * resolution)
    
    gmsh.model.mesh.generate(3)

gmsh.write("mesh3D.msh")

import pyvista as pv

mesh_target_for_plot = pv.read('mesh3D.msh')

labels = dict(xlabel='x', ylabel='y',zlabel='z')
plotter = pv.Plotter()
plotter.add_mesh(mesh_target_for_plot, show_edges=True, color="lightblue", opacity=0.5)
plotter.add_axes(**labels)
plotter.show()

gmsh.model.occ.synchronize()
gmsh.fltk.run()

[HoseongJeong]

The below is code using meshio, but I am not sure if the indices of edges are the same as those of scikit-fem.

mesh_1_support = meshio.read("mesh3D.msh") #스터럽의 노드번호를 찾기위해 meshio를 사용

physical_groups = mesh_1_support.field_data # Physical Group 이름과 ID 정보를 저장
curve_physical_group_id = physical_groups['Transverse_rebar line'][0] # 예를 들어, 'Transverse_rebar line'이라는 Physical Group을 찾고자 할 때
curve_edges = []
for cell_block, cell_data in zip(mesh_1_support.cells, mesh_1_support.cell_data["gmsh:physical"]):
    if cell_block.type == "line":  # 1D 요소 (Edge, Curve)
        for i, physical_group_id in enumerate(cell_data):
            if physical_group_id == curve_physical_group_id:
                curve_edges.append(cell_block.data[i]) #해당 edge들을 저장

curve_edges=np.array(curve_edges) #어레이로 변환
transverse_rebar = [] #스터럽 리스트 저장

found_indices = set()
for row in curve_edges:
    # 각 row를 set으로 변환하여 비교 (순서 무시)
    row_set = set(row)
    found = False
    for idx, edge in enumerate(mesh_1.edges.T):
        if row_set == set(edge) and idx not in found_indices:  # 중복된 인덱스 방지
            transverse_rebar.append(idx)  # 일치하는 인덱스를 저장
            found_indices.add(idx)  # 이미 찾은 인덱스는 추가하지 않음
            found = True
            break