Element Connectivity of Boundary Surfaces

[Aminofa70]

Dear Ferrite Developers,

How can we find the element connectivity of the boundary surface and its corresponding global cell id in the grid, assuming that we don't have facetsets, we need Ferrite.FacetIndex(element_id, local_face_id)
for example in the following code

Nx = 5;  Ny = 5 ; Nz = 5
Lx = 2.0; Ly = 2.0; Lz = 2.0
left = Ferrite.Vec(0.0, 0.0, 0.0)
right = Ferrite.Vec(Lx, Ly, Lz)
grid = generate_grid(Ferrite.Hexahedron, (Nx, Ny, Nz), left, right)

Thank you so much in advance

[termi-official]

I probably do not understand the question. The "boundary elements" do not exist in the grid you generate in your example code. You can get full connectivity information from e.g. the ExclusiveTopology, or you just reconstruct the information you need manually.

[fredrikekre]

Might be useful if you explain a bit more what you want to achieve in case there is a better solution (since this sounds a little bit like an XY problem).

[Aminofa70]

Thanks. The aim is to convert mesh in Comodo.jl https://github.com/COMODO-research/Comodo.jl to Ferrite for the finite element

The Comodo mesh for a simple case is:

using Comodo
using Ferrite
boxDim = [2, 2, 2]
boxEl = [6, 6, 6]
E, V, F, Fb, Cb = hexbox(boxDim, boxEl)

#E: is element connectivty equal to grid.cells
#V: is vertices equal to grid.nodes
#F: Entire quad faces in the solid
#Fb: only quad face in the boundary
#Cb: is the boundary marker  (for hex here is from 1 to 6)

 No I want to convert this mesh to Ferrite grid, first I get grid.cells and grid.nodes

cells = [Ferrite.Quadrilateral((e[1], e[2], e[3], e[4])) for e in F]
nodes = [Ferrite.Node((e[1], e[2])) for e in V]

Now How to Create facetsets form Comodo
I use the following code:

function create_facetsets(face_boundary, cells)
    strs = OrderedSet{FacetIndex}()
    for (element_id, element) in enumerate(cells)
        element_faces = Ferrite.faces(element)
        for (local_face_id, element_face) in enumerate(element_faces)
            # Check if this element face matches any boundary face
            for face_vertices in face_boundary
                if Set(element_face) == Set(face_vertices)
                    push!(strs, Ferrite.FacetIndex(element_id, local_face_id))
                    break  
                end
            end
        end
    end
    
    return strs
end

# based on Ferrite.jl
Fb_bottom = Fb[Cb.==1]
Fb_front = Fb[Cb.==3]
Fb_top = Fb[Cb.==2]
Fb_back = Fb[Cb.==4]
Fb_right = Fb[Cb.==5]
Fb_left = Fb[Cb.==6]

left = create_facetsets(Fb_left, cells)
bottom = create_facetsets(Fb_bottom, cells)
right = create_facetsets(Fb_right, cells)
back = create_facetsets(Fb_back, cells,)
top = create_facetsets(Fb_top, cells,)
front = create_facetsets(Fb_front, cells,)
facetsets = Dict("left" => left, "bottom" => bottom, "right" => right, "back" => back, "top" => top, "front" => front)

The I generate the grid:

Grid(cells, nodes, facetsets=facetsets)

The main question here is that in create_facetsetsfunction I do loop for all elements and when the mesh is very dense it takes a lot of time.
Instead I want to do loop only for those element that contains the faces in the boundary instead of all element to be more efficient.
So, I don't know how to find those elements?

I hope now it is more clearer

[KristofferC]

This is just a standard programming problem where you are doing it in a very inefficient way. You could put this comment into an LLM, and it will spit out a better way of doing it.

[Aminofa70]

Of course cause I'm not a computer engineer, and unfortunately, LLM does not understand Ferrite as well.

[AbdAlazezAhmed]

Maybe I misunderstood something but doesn't hexbox generate hex not quad?
Also, using cells = [Ferrite.Quadrilateral((e[1], e[2], e[3], e[4])) for e in F] will result in duplicate elements iiuc because F contains the faces local to the elements not globally. Also note that the mesh in the example is in 3D so nodes = [Ferrite.Node((e[1], e[2])) for e in V] will likely be problematic.
A potential answer to your question which ignores the apriori knowledge of boundary facets is using Ferrite.addboundaryfacetset!.
Another solution which is probably suboptimal but that's what I could think of as for now is that you iterate over Fb once and findnext the corresponding entry in F starting from the last entry found. iiuc this should work if Fb is just a filtered F preserving order.
MWE:

function create_facetsets(Fb, F, Cb)
                 strs = [OrderedSet{FacetIndex}() for _ in 1:6]
                  for i in 1:6
                   sizehint!(strs[i], count(x->x==i, Cb))
                  end
                  let last_node = 1
                      for (i, boundary_face) in enumerate(Fb)
                          global_face_idx = findnext(x->x==boundary_face, F, last_node)
                          # Don't forget to check for nothing, I'm too lazy to do it
                          last_node = global_face_idx
                          cell_id = (global_face_idx - 1) ÷ 6 + 1
                          #assuming comodo uses same local node and face ordering as ferrite
                          local_facet_id = (global_face_idx - 1) % 6 + 1
                          push!(strs[Cb[i]], Ferrite.FacetIndex(cell_id, local_facet_id))
                      end
                  end

     return strs
 end

[Aminofa70]

@AbdAlazezAhmed Thanks a lot. Sorry I made a mistake for the elements and nodes, the correct ones are:

cells = [Ferrite.Hexahedron((e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8])) for e in E]
nodes = [Ferrite.Node((e[1], e[2], e[3])) for e in V]

I created the following code and works well. This is just for hex mesh and I will modify it for other elements
It is also fast for the dense mesh. I don't know if the code is efficient but it works quite well.

using Comodo
using Ferrite
using OrderedCollections
using Comodo
using Comodo.GLMakie
using Comodo.GLMakie.Colors
using Comodo.GeometryBasics
using Comodo.Statistics

using ComodoFerrite
using ComodoFerrite.Ferrite

GLMakie.closeall()

boxDim = [2, 2, 2]
boxEl = [10, 10, 10]
E, V, F, Fb, Cb = hexbox(boxDim, boxEl)

# Convert to Ferrite types
cells = [Ferrite.Hexahedron(Tuple(e)) for e in E]
nodes = [Ferrite.Node(Tuple(v)) for v in V]
grid = Grid(cells, nodes)

function addface!(grid, name::String, boundary_faces)
    # Create a mapping from face vertices to boundary face info
    face_to_idx = Dict{NTuple{4, Int}, Int}()
    for (i, face) in enumerate(boundary_faces)
        sorted_face = Tuple(sort(collect(face)))
        face_to_idx[sorted_face] = i
    end
    
    # Get unique node indices from boundary faces
    boundary_nodes = Set{Int}()
    for face in boundary_faces
        union!(boundary_nodes, face)
    end
    
    # Find elements that contain any boundary node
    candidate_elements = Set{Int}()
    for (element_id, cell) in enumerate(grid.cells)
        cell_nodes = Ferrite.get_node_ids(cell)
        if !isdisjoint(cell_nodes, boundary_nodes)
            push!(candidate_elements, element_id)
        end
    end
    
    # Create new facetset - only loop through candidate elements
    new_facetset = Set{FacetIndex}()
    for element_id in candidate_elements
        cell = grid.cells[element_id]
        element_faces = Ferrite.faces(cell)
        for (local_face_id, element_face) in enumerate(element_faces)
            sorted_face = Tuple(sort(collect(element_face)))
            if haskey(face_to_idx, sorted_face)
                push!(new_facetset, FacetIndex(element_id, local_face_id))
            end
        end
    end
    
    # Add to grid's facetsets
    grid.facetsets[name] = new_facetset
    
    return grid
end

# Usage example:
Fb_bottom = Fb[Cb.==1]
addface!(grid, "bottom", Fb_bottom)

GLMakie.closeall()

M = GeometryBasics.Mesh(V, F, normal = face_normals(V, F))
Mb = GeometryBasics.Mesh(V, Fb, normal = face_normals(V, Fb))

fig = Figure(size = (1200,800))

ax1 = AxisGeom(fig[1, 1], title = "Tet4 mesh", azimuth = -0.2π, elevation = 0.1π)
hp1 = meshplot!(ax1, M ,color=:gray,  strokecolor=:black, strokewidth=3.0, shading =  false, transparency = false)

ax2 = AxisGeom(fig[1, 2], title = "Boundary condition", azimuth = -0.2π, elevation = 0.1π)
hp2 = meshplot!(ax2, Mb, color=(Gray(0.95), 0.3),  strokecolor=:black, strokewidth=2.0, shading=true, transparency = true)
 
## Visualize boundary condition
### First convert the face_index and node_index to points
facesset = get_boundary_points(grid, getfacetset(grid, "bottom"), Faces, Ferrite.Tetrahedron)

scatter!(ax2, facesset, color=:blue,markersize=15.0, marker=:circle, strokecolor=:black, strokewidth=2, label = "Traction")
axislegend(ax2, position=:rb, backgroundcolor=(:white, 0.7), framecolor=:gray)

display(fig)

grid.facetsets