Extract boundary entities

[termi-official]

Following the discussion in #565 it is helpful to add utilities to extract boundary entities (faces, edges, ...). Boundary faces can be identified as a function of the grid and a topology. So maybe something like addboundaryfaceset!(grid, topology, "bottom" x -> norm(x[3]) ≈ 0.0); makes sense.

[graphitical]

Hi Dennis, the topology you list here, would that be an ExclusiveTopology or something else?

[termi-official]

Yes. Currently only the ExclusiveTopology is implemented in Ferrite. Note that the topology interface is still rather experimental. The topology is decoupled from the mesh, because this allows to have different neighborhood construction and search algorithms, depending on the specific use-case of a user.

Edit: You also might have noted that the AbstractTopology interface is not well-defined yet. We want the ExclusiveTopology to stabilize first and then derive an abstract interface from it.

[graphitical]

Thanks. Would the topology argument you list originally be necessary? I'm a bit lost on how all the pieces work together. It seems like grid and the defining boundary function should contain all the information necessary to create a boundary.

[termi-official]

The topology here is just ment to help to identify which entities are not in the interior of the grid. Note that the grid has basically all the information that is also in the topology. We precompute topological information from the grid to speed up queries like for example "Who is my neighbor on face 1", which is stored in such topology data structures. This makes boundary queries easy, because a boundary face can be defined as a face without neighbor face in a suitable topology. A solution which does not use any topology data structure is also fine and possible.

[graphitical]

I think I see. So you were originally suggesting creating a new type of topology (i.e., ExteriorTopology) which would contain all the topological information for all the exterior entities.

The intent of this would be to reduce the search space for checking if a node satisfies the query function.

Is that correct?

This is my understanding because ExclusiveTopology seems to contain topological information for the entire grid so that obviously contains a lot of points we don't want to check against.

[koehlerson]

I don't think this is exactly what Dennis meant. You can just use the ExclusiveTopology and query that by

julia> grid = generate_grid(Quadrilateral,(100,100))
Grid{2, Quadrilateral, Float64} with 10000 Quadrilateral cells and 10201 nodes

julia> topology = ExclusiveTopology(grid);

julia> findall(isempty,topology.face_neighbor)
400-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 1)
 CartesianIndex(3, 1)
 CartesianIndex(4, 1)
 CartesianIndex(5, 1)
 CartesianIndex(6, 1)
 CartesianIndex(7, 1)
 CartesianIndex(8, 1)
 CartesianIndex(9, 1)
 CartesianIndex(10, 1)
 CartesianIndex(11, 1)
 CartesianIndex(12, 1)
 CartesianIndex(13, 1)
 CartesianIndex(14, 1)
 CartesianIndex(15, 1)
 CartesianIndex(16, 1)
 CartesianIndex(17, 1)
 CartesianIndex(18, 1)
 CartesianIndex(19, 1)
 CartesianIndex(20, 1)
 CartesianIndex(21, 1)
 CartesianIndex(22, 1)
 CartesianIndex(23, 1)
 CartesianIndex(24, 1)
 ⋮
 CartesianIndex(7701, 4)
 CartesianIndex(7801, 4)
 CartesianIndex(7901, 4)
 CartesianIndex(8001, 4)
 CartesianIndex(8101, 4)
 CartesianIndex(8201, 4)
 CartesianIndex(8301, 4)
 CartesianIndex(8401, 4)
 CartesianIndex(8501, 4)
 CartesianIndex(8601, 4)
 CartesianIndex(8701, 4)
 CartesianIndex(8801, 4)
 CartesianIndex(8901, 4)
 CartesianIndex(9001, 4)
 CartesianIndex(9101, 4)
 CartesianIndex(9201, 4)
 CartesianIndex(9301, 4)
 CartesianIndex(9401, 4)
 CartesianIndex(9501, 4)
 CartesianIndex(9601, 4)
 CartesianIndex(9701, 4)
 CartesianIndex(9801, 4)
 CartesianIndex(9901, 4)

The CartesianIndex corresponds to the element id and the second entry is the local face id and can be directly translated to a FaceIndex

[AbdAlazezAhmed]

Hi, I am working on an implementation (shown below)

addboundaryfaceset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function; all::Bool=true) = 
    _addset!(grid, topology, name, f, Ferrite.faces, grid.facesets, FaceIndex; all=all)

function _addset!(grid::AbstractGrid, topology::ExclusiveTopology, name::String, f::Function, _ftype::Function, dict::Dict, _indextype::Type; all::Bool=true)
    _check_setname(dict, name)
    _set = Set{_indextype}()
    cells = getcells(grid)
    pass = all
    for boundary in findall(isempty,topology.face_neighbor)
        cell_idx = boundary[1]
        face_idx = boundary[2]
        face = _ftype(cells[cell_idx])[face_idx]
        for node_idx in face
            v = f(grid.nodes[node_idx].x)
            all ? (!v && (pass = false; break)) : (v && (pass = true; break))
        end
        pass && push!(_set, _indextype(cell_idx, face_idx))
    end
    _warn_emptyset(_set, name)
    dict[name] = _set
    grid
end

Am I going in the right direction? Will other functions be needed (addboundaryedgeset!, addboundaryvertexset) ?

Thanks.

[termi-official]

@graphitical : Sorry. As Maxi wrote, my idea was to just use the ExclusiveTopology to query which faces are on the boundary. If you have a more efficient way to extract the boundary with the information contained in the grid, then this is even better!

@AbdAlazezAhmed : Indeed, this is the right direction. If you want, then you could also try to implement your suggestions regarding boundary edges and boundary vertices. It might be interesting to implement some helpers for this to efficiently extract edges/vertices on a specific face/edge and add it to the public API.

[AbdAlazezAhmed]

Thank you, I'll implement addboundaryedgeset! and addboundaryvertexset! with the helpers.

[graphitical]

@termi-official thanks! I'm still learning the API and didn't realize ExclusiveTopology would do that. Looks like @AbdAlazezAhmed is already well on the case so maybe I'll go digging around for a different thing to contribute on.

[termi-official]

Also thanks for taking time here @graphitical and @AbdAlazezAhmed . A bit unfortunate that two people grabbed the issue at the same time, but I hope it is okay.

The topology API is still very experimental and I understand that it might be hard to use at the moment due to the lack of examples and scarce docs.

Do you have something specific in mind? You can ping me on Slack or Zulip if you need help.

[graphitical]

I don't have anything specific in mind I suppose, but things that will help me with my computational homogenization task are on my radar.

Everyone's been very generous with help on Slack so I was looking for a way to contribute. I've also never contributed to a project before so I'm still figuring out all the logistics.

[AbdAlazezAhmed]

Hi, I made a pull request with the changes. I'm not sure but I think the helpers getfacevertices and getedgevertices are redundant as the face/edge is a set of vertices (or I misunderstood what the helper should be). I added documentation similar to the existing functions and made a simple test with quadrilateral grid. (tested Tetrahedron manually) Please let me know if I got something wrong as it's my first time contributing to a project. Thanks!