Deprecate onboundary, remove boundary_matrix property from Grid (…

…#924)

This patch deprecates the `onboundary` function which checks whether an
element has any facet alone an external boundary. This functionality
relied on the `boundary_matrix` property in `Grid`, but this was only
populated correctly by the builtin grid generators (and not e.g. for
meshes parsed from other formats).

Checking `onboundary` was just a shortcut to avoid some extra work in
the facetset containment check which you always had to do anyways, but
since this is just a tuple hash it isn't expensive anyway.

In addition, nowadays there are utilities such as `FacetIterator` which
let's you explicitly loop over the facet sets you want directly.

CHANGELOG.md

@@ -385,6 +385,17 @@ more discussion).
   now you can still use them by prefixing `Ferrite.`, e.g. `Ferrite.getweights`.)
   ([#754][github-754])
 
+- The `onboundary` function (and the associated `boundary_matrix` property of the `Grid`
+  datastructure) have been removed ([#924][github-924]). Instead of first checking
+  `onboundary` and then check whether a facet belong to a specific facetset, check the
+  facetset directly. For example:
+  ```diff
+  - if onboundary(cell, local_face_id) && (cell_id, local_face_id) in getfacesets(grid, "traction_boundary")
+  + if (cell_id, local_face_id) in getfacesets(grid, "traction_boundary")
+       # integrate the "traction_boundary" boundary
+    end
+  ```
+
 ### Fixed
 
 - Benchmarks now work with master branch. ([#751][github-#751], [#855][github-#855])
@@ -902,3 +913,4 @@ poking into Ferrite internals:
 [github-835]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/835
 [github-855]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/855
 [github-880]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/880
+[github-924]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/924

docs/src/literate-gallery/helmholtz.jl

@@ -135,9 +135,8 @@ function doassemble(cellvalues::CellValues, facetvalues::FacetValues,
         # ```
         #+
         for facet in 1:nfacets(cell)
-            if onboundary(cell, facet) && 
-                   ((cellcount, facet) ∈ getfacetset(grid, "left") || 
-                    (cellcount, facet) ∈ getfacetset(grid, "bottom"))
+            if (cellcount, facet) ∈ getfacetset(grid, "left") ||
+               (cellcount, facet) ∈ getfacetset(grid, "bottom")
                 reinit!(facetvalues, cell, facet)
                 for q_point in 1:getnquadpoints(facetvalues)
                     coords_qp = spatial_coordinate(facetvalues, q_point, coords)

docs/src/literate-gallery/topology_optimization.jl

@@ -351,7 +351,7 @@ function elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
     symmetrize_lower!(Ke)
 
     @inbounds for facet in 1:nfacets(getcells(grid, cellid(element))) 
-        if onboundary(element, facet) && (cellid(element), facet) ∈ getfacetset(grid, "traction")
+        if (cellid(element), facet) ∈ getfacetset(grid, "traction")
             reinit!(facetvalues, element, facet)
             t = Vec((0.0, -1.0)) # force pointing downwards
             for q_point in 1:getnquadpoints(facetvalues)

docs/src/literate-tutorials/incompressible_elasticity.jl

@@ -158,7 +158,7 @@ function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, g
     ## We loop over all the facets in the cell, then check if the facet
     ## is in our `"traction"` facetset.
     for facet in 1:nfacets(cell)
-        if onboundary(cell, facet) && (cellid(cell), facet) ∈ getfacetset(grid, "traction")
+        if (cellid(cell), facet) ∈ getfacetset(grid, "traction")
             reinit!(facetvalues_u, cell, facet)
             for q_point in 1:getnquadpoints(facetvalues_u)
                 dΓ = getdetJdV(facetvalues_u, q_point)

src/Grid/grid.jl

@@ -292,7 +292,6 @@ There are multiple helper structures to apply boundary conditions or define subd
 - `nodesets::Dict{String,Set{Int}}`: maps a `String` key to a `Set` of global node ids
 - `facetsets::Dict{String,Set{FacetIndex}}`: maps a `String` to a `Set` of `Set{FacetIndex} (global_cell_id, local_facet_id)`
 - `vertexsets::Dict{String,Set{VertexIndex}}`: maps a `String` key to a `Set` of local vertex ids
-- `boundary_matrix::SparseMatrixCSC{Bool,Int}`: optional, only needed by `onboundary` to check if a cell is on the boundary, see, e.g. Helmholtz example
 """
 mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
     cells::Vector{C}
@@ -302,8 +301,6 @@ mutable struct Grid{dim,C<:AbstractCell,T<:Real} <: AbstractGrid{dim}
     nodesets::Dict{String,Set{Int}}
     facetsets::Dict{String,Set{FacetIndex}}
     vertexsets::Dict{String,Set{VertexIndex}}
-    # Boundary matrix (faces per cell × cell)
-    boundary_matrix::SparseMatrixCSC{Bool,Int} # TODO: Deprecate!
 end
 
 function Grid(cells::Vector{C},
@@ -313,7 +310,7 @@ function Grid(cells::Vector{C},
               facetsets::Dict{String,Set{FacetIndex}}=Dict{String,Set{FacetIndex}}(),
               facesets = nothing,
               vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
-              boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
+              boundary_matrix = nothing) where {dim,C,T}
     if facesets !== nothing 
         if isempty(facetsets)
             @warn "facesets in Grid is deprecated, use facetsets instead" maxlog=1
@@ -324,7 +321,10 @@ function Grid(cells::Vector{C},
             error("facesets are deprecated, use only facetsets")
         end
     end
-    return Grid(cells, nodes, cellsets, nodesets, facetsets, vertexsets, boundary_matrix)
+    if boundary_matrix !== nothing
+        error("`boundary_matrix` is not part of the Grid anymore and thus not a supported keyword argument.")
+    end
+    return Grid(cells, nodes, cellsets, nodesets, facetsets, vertexsets)
 end
 
 ##########################

src/Grid/grid_generators.jl

@@ -1,17 +1,3 @@
-function boundaries_to_sparse(boundary)
-    n = length(boundary)
-    I = Vector{Int}(undef, n)
-    J = Vector{Int}(undef, n)
-    V = Vector{Bool}(undef, n)
-    for (idx, el) in enumerate(boundary)
-        cell, face = el.idx
-        I[idx] = face
-        J[idx] = cell
-        V[idx] = true
-    end
-    return sparse(I, J, V)
-end
-
 """
     generate_grid(celltype::Cell, nel::NTuple, [left::Vec, right::Vec)
 
@@ -45,12 +31,10 @@ function generate_grid(::Type{Line}, nel::NTuple{1,Int}, left::Vec{1,T}=Vec{1}((
     boundary = Vector([FacetIndex(1, 1),
                        FacetIndex(nel_x, 2)])
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     facetsets = Dict("left"  => Set{FacetIndex}([boundary[1]]),
                     "right" => Set{FacetIndex}([boundary[2]]))
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 # QuadraticLine
@@ -75,12 +59,10 @@ function generate_grid(::Type{QuadraticLine}, nel::NTuple{1,Int}, left::Vec{1,T}
     boundary = FacetIndex[FacetIndex(1, 1),
                           FacetIndex(nel_x, 2)]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     facetsets = Dict("left"  => Set{FacetIndex}([boundary[1]]),
                     "right" => Set{FacetIndex}([boundary[2]]))
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 function _generate_2d_nodes!(nodes, nx, ny, LL, LR, UR, UL)
@@ -139,8 +121,6 @@ function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T},
                               [FacetIndex(cl, 3) for cl in cell_array[:,end]];
                               [FacetIndex(cl, 4) for cl in cell_array[1,:]]]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String, Set{FacetIndex}}()
@@ -149,7 +129,7 @@ function generate_grid(C::Type{Quadrilateral}, nel::NTuple{2,Int}, LL::Vec{2,T},
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 # QuadraticQuadrilateral
@@ -178,8 +158,6 @@ function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2,Int}, LL::V
                               [FacetIndex(cl, 3) for cl in cell_array[:,end]];
                               [FacetIndex(cl, 4) for cl in cell_array[1,:]]]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String, Set{FacetIndex}}()
@@ -188,7 +166,7 @@ function generate_grid(::Type{QuadraticQuadrilateral}, nel::NTuple{2,Int}, LL::V
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[:,end])) .+ offset]); offset += length(cell_array[:,end])
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,:]))   .+ offset]); offset += length(cell_array[1,:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 # Hexahedron
@@ -223,8 +201,6 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve
                               [FacetIndex(cl, 5) for cl in cell_array[1,:,:][:]];
                               [FacetIndex(cl, 6) for cl in cell_array[:,:,end][:]]]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String,Set{FacetIndex}}()
@@ -235,7 +211,7 @@ function generate_grid(::Type{Hexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Ve
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end 
 
 # Wedge
@@ -273,8 +249,6 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}(
     @views bo = [map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> FacetIndex(x,1), c_nxyz[2, :, :, 1][:])]
     @views to = [map(x -> FacetIndex(x,5), c_nxyz[1, :, :, end][:]) ; map(x -> FacetIndex(x,5), c_nxyz[2, :, :, end][:])]
 
-    boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])
-
     facetsets = Dict(
         "left" => Set(le),
         "right" => Set(ri),
@@ -284,7 +258,7 @@ function generate_grid(::Type{Wedge}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}(
         "top" => Set(to),
     )
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end 
 
 #Pyramid
@@ -336,8 +310,6 @@ function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3
     @views bo = map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:]) 
     @views to = map(x -> FacetIndex(x,1), c_nxyz[6, :, :, end][:])
 
-    boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])
-
     facetsets = Dict(
         "left" => Set(le),
         "right" => Set(ri),
@@ -347,7 +319,7 @@ function generate_grid(::Type{Pyramid}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3
         "top" => Set(to),
     )
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end 
 
 function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTuple{3,Int}, left::Vec{3,T}=Vec{3}((-1.0,-1.0,-1.0)), right::Vec{3,T}=Vec{3}((1.0,1.0,1.0))) where {T}
@@ -393,8 +365,6 @@ function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTup
                               [FacetIndex(cl, 5) for cl in cell_array[1,:,:][:]];
                               [FacetIndex(cl, 6) for cl in cell_array[:,:,end][:]]]
 
-    boundary_matrix = Ferrite.boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String,Set{FacetIndex}}()
@@ -405,7 +375,7 @@ function Ferrite.generate_grid(::Type{SerendipityQuadraticHexahedron}, nel::NTup
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,:,:][:]))   .+ offset]); offset += length(cell_array[1,:,:][:])
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[:,:,end][:])) .+ offset]); offset += length(cell_array[:,:,end][:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)   
+    return Grid(cells, nodes, facetsets=facetsets)   
 end
 
 # Triangle
@@ -433,8 +403,6 @@ function generate_grid(::Type{Triangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::V
                                [FacetIndex(cl, 2) for cl in cell_array[2,:,end]];
                                [FacetIndex(cl, 3) for cl in cell_array[1,1,:]]]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String,Set{FacetIndex}}()
@@ -443,7 +411,7 @@ function generate_grid(::Type{Triangle}, nel::NTuple{2,Int}, LL::Vec{2,T}, LR::V
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 # QuadraticTriangle
@@ -473,8 +441,6 @@ function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2,Int}, LL::Vec{2,
                               [FacetIndex(cl, 2) for cl in cell_array[2,:,end]];
                               [FacetIndex(cl, 3) for cl in cell_array[1,1,:]]]
 
-    boundary_matrix = boundaries_to_sparse(boundary)
-
     # Cell face sets
     offset = 0
     facetsets = Dict{String,Set{FacetIndex}}()
@@ -483,7 +449,7 @@ function generate_grid(::Type{QuadraticTriangle}, nel::NTuple{2,Int}, LL::Vec{2,
     facetsets["top"]    = Set{FacetIndex}(boundary[(1:length(cell_array[2,:,end])) .+ offset]); offset += length(cell_array[2,:,end])
     facetsets["left"]   = Set{FacetIndex}(boundary[(1:length(cell_array[1,1,:]))   .+ offset]); offset += length(cell_array[1,1,:])
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end
 
 # Tetrahedron
@@ -552,8 +518,6 @@ function generate_grid(::Type{Tetrahedron}, cells_per_dim::NTuple{3,Int}, left::
     @views bo = [map(x -> FacetIndex(x,1), c_nxyz[1, :, :, 1][:])   ; map(x -> FacetIndex(x,1), c_nxyz[3, :, :, 1][:])]
     @views to = [map(x -> FacetIndex(x,3), c_nxyz[5, :, :, end][:]) ; map(x -> FacetIndex(x,3), c_nxyz[6, :, :, end][:])]
 
-    boundary_matrix = boundaries_to_sparse([le; ri; bo; to; fr; ba])
-
     facetsets = Dict(
         "left" => Set(le),
         "right" => Set(ri),
@@ -563,5 +527,5 @@ function generate_grid(::Type{Tetrahedron}, cells_per_dim::NTuple{3,Int}, left::
         "top" => Set(to),
     )
 
-    return Grid(cells, nodes, facetsets=facetsets, boundary_matrix=boundary_matrix)
+    return Grid(cells, nodes, facetsets=facetsets)
 end

src/deprecations.jl

@@ -364,3 +364,8 @@ function addfaceset!(grid, name, f::Function; kwargs...)
     @warn "addfaceset! is deprecated, using addfacetset! instead"
     return addfacetset!(grid, name, f; kwargs...)
 end
+
+export onboundary
+function onboundary(::CellCache, ::Int)
+    error("`onboundary` is deprecated, check just the facetset instead of first checking `onboundary`.")
+end

src/exports.jl

@@ -92,7 +92,6 @@ export
     get_node_coordinate,
     getcoordinates,
     getcoordinates!,
-    onboundary,
     nfacets,
     addnodeset!,
     addfacetset!,

src/iterators.jl

@@ -103,7 +103,6 @@ celldofs!(v::Vector, cc::CellCache) = copyto!(v, cc.dofs) # celldofs!(v, cc.dh,
 
 # TODO: These should really be replaced with something better...
 nfacets(cc::CellCache) = nfacets(cc.grid.cells[cc.cellid[]])
-onboundary(cc::CellCache, face::Int) = cc.grid.boundary_matrix[face, cc.cellid[]]
 
 
 # TODO: Currently excluded from the docstring below. Should they be public?
@@ -155,7 +154,6 @@ for op = (:getnodes, :getcoordinates, :cellid, :celldofs)
 end
 # @inline faceid(fc::FacetCache) = fc.current_faceid[]
 @inline celldofs!(v::Vector, fc::FacetCache) = celldofs!(v, fc.cc)
-# @inline onboundary(fc::FacetCache) = onboundary(fc.cc, faceid(fc))
 # @inline faceindex(fc::FacetCache) = FaceIndex(cellid(fc), faceid(fc))
 @inline function reinit!(fv::FacetValues, fc::FacetCache)
     reinit!(fv, fc.cc, fc.current_facet_id[])

test/test_deprecations.jl

@@ -101,4 +101,11 @@ end
     @test_throws ErrorException vtk_grid("old", generate_grid(Line, (1,)))
 end
 
+@testset "onboundary" begin
+    msg = "`onboundary` is deprecated, check just the facetset instead of first checking `onboundary`."
+    @test_throws ErrorException(msg) onboundary(first(CellIterator(generate_grid(Line, (2,)))), 1)
+    msg = "`boundary_matrix` is not part of the Grid anymore and thus not a supported keyword argument."
+    @test_throws ErrorException(msg) Grid(Triangle[], Node{2,Float64}[]; boundary_matrix = something)
+end
+
 end # testset deprecations