Create simple mesh with Gridap

Lx, Ly, Lz = (5.0, 5.0, 1.0)
nx, ny, nz = (10, 10, 2)
domain = (0.0,Lx, 0.0,Ly, 0.0,Lz)
partition = (nx,ny,nz)
model = CartesianDiscreteModel(domain, partition)

d = num_cell_dims(model)
Γ = BoundaryTriangulation(model)
labels = get_face_labeling(model)


# 各境界要素の重心座標を取る
coords = get_cell_coordinates(Γ)
# 各セルの重心 (VectorValue{3,Float64})
centers = [ mean(c) for c in coords ]
tol = 1e-5
tol_rough = 0.3
# surf1_0: y=Ly, x=0 の稜線
surf1_1_ids = [
    i for (i,c) in enumerate(centers)
    if abs(c[2]-Ly)<tol_rough && abs(c[1]-0.0)<tol
]

# surf1_1: y=Ly, x=Lx の稜線
surf1_2_ids = [
    i for (i,c) in enumerate(centers)
    if abs(c[2]-Ly)<tol_rough && abs(c[1]-Lx)<tol
]

# surf2: y=0 の全面
surf2_ids = [
    i for (i,c) in enumerate(centers)
    if abs(c[2]-0.0)<tol
]

# タグ登録
add_tag!(labels, "surface_1_1", surf1_1_ids)
add_tag!(labels, "surface_1_2", surf1_2_ids)
add_tag!(labels, "surface_2", surf2_ids)

# 確認用：タグ名が正しく登録されているかをチェック
println("Registered Tags: ", labels.tag_to_name)

# 正しくタグが見えるなら、取得可能
ids_s1_1 = Gridap.Geometry.get_tag_entities(labels, "surface_1_1")
ids_s1_2 = Gridap.Geometry.get_tag_entities(labels, "surface_1_2")
ids_s2 = Gridap.Geometry.get_tag_entities(labels, "surface_2")

# 取得結果の長さなど確認
println("surface_1_1 entity count: ", length(ids_s1_1))
println("surface_1_2 entity count: ", length(ids_s1_2))
println("surface_2 entity count: ", length(ids_s2))


# Triangulation のセル数に合わせた Bool → Int 配列へ
ns = num_cells(Γ)
mask_s1_1 = falses(ns); mask_s1_1[ids_s1_1] .= true
mask_s1_2 = falses(ns); mask_s1_2[ids_s1_2] .= true
mask_s2 = falses(ns); mask_s2[ids_s2] .= true
tags_s1_1 = Int.(mask_s1_1)
tags_s1_2 = Int.(mask_s1_2)
tags_s2 = Int.(mask_s2)

# VTK に一括出力
writevtk(
    Γ, "boundary_surfaces",
    cellfields = [
    "surface_1_1" => tags_s1_1,
    "surface_1_2" => tags_s1_2,
    "surface_2" => tags_s2
    ]
)