Labeling Boundary with Gridap and export as VTU. What is the feature of Julia?

Data Structure in FaceLabeling

I struggled with ChatGPT to resolve bugs in my code, but finally gave up. It just tells us codes without thinking so much. What I understood is that seeing original code is fastest way. Before seeing the code, I was so reluctant to do that, seeing codes, which is implemented by others, is kind of headache. But What I found is that data structure and functions are so simple. Probably because this is a feature of Programing Language, Julia.

Lets see key data structure to label boundary, this is actual implementation that original repository uses.

struct FaceLabeling <: GridapType
    d_to_dface_to_entity :: Vector{Vector{Int32}}
    tag_to_entities      :: Vector{Vector{Int32}}
    tag_to_name          :: Vector{String}

I do not write detail here, but you can easily imagine what data is stored here.

Fundamental Usage to access mesh data

Get Entities

labels = get_face_labeling(model)
tag_A   = Gridap.Geometry.get_tag_from_name(labels, "surface_A")
ids_A   = Gridap.Geometry.get_tag_entities(labels, tag_A)

Some useful functions

d = num_cell_dims(model) - 1
topo  = Gridap.Geometry.get_grid_topology(model)
is_bdry = Gridap.Geometry.get_isboundary_face(topo, d)
verticies = Gridap.ReferenceFEs.get_face_vertices(topo, d)
cell_faces = get_cell_faces(topo)
face_entity = get_face_entity(labels, d)
nc = num_cells(model)

When I see original repo, the implementation of function get_face_entity is like this.

function get_face_entity(lab::FaceLabeling,d::Integer)
  lab.d_to_dface_to_entity[d+1]
end

What I found is, Julia is a language that combines small functions and structures to create larger processes.

Functions I coded

Export Boundary with Tag Information

To verify what tags are assigned on the boundary is important.

function write_boundary_tag_ids(
  model;
  tag_names::Vector{String},
  outname::String = "boundary_tags",
  background::Int = 0,
)
  D     = num_cell_dims(model)
  d     = D - 1 # facet
  topo  = Gridap.Geometry.get_grid_topology(model)
  labels = Gridap.Geometry.get_face_labeling(model)

  is_bdry = Gridap.Geometry.get_isboundary_face(topo, d)
  bdry_gface_ids = findall(is_bdry)

  tag_ids  = [Gridap.Geometry.get_tag_from_name(labels, nm) for nm in tag_names]
  masks    = [Gridap.Geometry.get_face_mask(labels, tid, d) for tid in tag_ids]  # ::Vector{Vector{Bool}}

  out = fill(Int(background), length(bdry_gface_ids))
  for (i, m) in enumerate(masks)
    m_bdry = m[bdry_gface_ids]
    for (k, flag) in pairs(m_bdry)
      if flag && out[k] == background
        out[k] = i
      end
    end
  end

  Γall = Gridap.Geometry.BoundaryTriangulation(model)
  Gridap.writevtk(
    Γall, outname;
    cellfields = ["tag_id" => Gridap.CellField(out, Γall)],
  )

  return nothing
end

Export Info with Independent Tag

function write_boundary_tag_ids2(
  model;
  tag_names::Vector{String},
  outname::String = "boundary_tags",
  background::Int = 0,
)
  D      = Gridap.Geometry.num_cell_dims(model)
  d      = D - 1
  topo   = Gridap.Geometry.get_grid_topology(model)
  labels = Gridap.Geometry.get_face_labeling(model)

  is_bdry         = Gridap.Geometry.get_isboundary_face(topo, d)
  bdry_gface_ids  = findall(is_bdry)  # ::Vector{Int}

  tag_ids = [Gridap.Geometry.get_tag_from_name(labels, nm) for nm in tag_names]
  masks   = [Gridap.Geometry.get_face_mask(labels, tid, d) for tid in tag_ids]  # Vector{Vector{Bool}}

  Γall = Gridap.Geometry.BoundaryTriangulation(model)

  cellfields_int = [
    begin
      m_bdry = m[bdry_gface_ids]
      vals   = ifelse.(m_bdry, Int(tid), background)  # Vector{Int}
      tag => Gridap.CellField(vals, Γall)
    end
    for (m, tid, tag) in zip(masks, tag_ids, tag_names)
  ]

  Gridap.writevtk(Γall, outname; cellfields = cellfields_int)
  return nothing
end