Merge branch 'master' into do/flexible-sparsity-pattern

CHANGELOG.md

@@ -319,6 +319,10 @@ more discussion).
   **To upgrade** replace the quadrature rule passed to `FaceValues` with a
   `FaceQuadratureRule`.
 
+- Checking if a face `(ele_id, local_face_id) ∈ faceset` has been previously implemented
+  by type piracy. In order to be invariant to the underlying `Set` datatype as well as
+  omitting type piracy, ([#835][github-835]) implemented `isequal` and `hash` for `BoundaryIndex` datatypes.
+
 ### Deprecated
 
 - The rarely (if ever) used methods of `function_value`, `function_gradient`,
@@ -867,3 +871,4 @@ poking into Ferrite internals:
 [github-749]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/749
 [github-753]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/753
 [github-756]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/756
+[github-835]: https://github.com/Ferrite-FEM/Ferrite.jl/pull/835

docs/src/devdocs/interpolations.md

@@ -15,6 +15,9 @@ Ferrite.shape_gradient(::Interpolation, ::Vec, ::Int)
 Ferrite.shape_gradient_and_value
 Ferrite.boundarydof_indices
 Ferrite.dirichlet_boundarydof_indices
+Ferrite.shape_values!
+Ferrite.shape_gradients!
+Ferrite.shape_gradients_and_values!
 ```
 
 ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element

docs/src/reference/dofhandler.md

@@ -12,6 +12,7 @@ SubDofHandler
 ## Adding fields to the DofHandlers
 ```@docs
 add!(::DofHandler, ::Symbol, ::Interpolation)
+add!(::SubDofHandler, ::Symbol, ::Interpolation)
 close!(::DofHandler)
 ```
 

src/Dofs/DofHandler.jl

@@ -12,7 +12,6 @@ Access some grid representation for the dof handler.
 """
 get_grid(dh::AbstractDofHandler)
 
-
 struct SubDofHandler{DH} <: AbstractDofHandler
     # From constructor
     dh::DH
@@ -26,8 +25,36 @@ struct SubDofHandler{DH} <: AbstractDofHandler
     # const dof_ranges::Vector{UnitRange{Int}} # TODO: Why not?
 end
 
-# TODO: Should be an inner constructor.
+"""
+    SubDofHandler(dh::AbstractDofHandler, cellset::Set{Int})
+
+Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the 
+cells in `cellset`. This allows you to add fields to parts of the domain, or using 
+different interpolations or cell types (e.g. `Triangles` and `Quadrilaterals`). All 
+fields and cell types must be the same in one `SubDofHandler`.
+
+After construction any number of discrete fields can be added to the SubDofHandler using
+[`add!`](@ref). Construction is finalized by calling [`close!`](@ref) on the parent `dh`.
+
+# Examples
+We assume we have a `grid` containing "Triangle" and "Quadrilateral" cells, 
+including the cellsets "triangles" and "quadilaterals" for to these cells. 
+```julia
+dh = DofHandler(grid)
+
+sdh_tri = SubDofHandler(dh, getcellset(grid, "triangles"))
+ip_tri = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u
+add!(sdh_tri, :u, ip_tri)
+
+sdh_quad = SubDofHandler(dh, getcellset(grid, "quadilaterals"))
+ip_quad = Lagrange{RefQuadrilateral, 2}()^2 # vector interpolation for a field u
+add!(sdh_quad, :u, ip_quad)
+
+close!(dh) # Finalize by closing the parent 
+```
+"""
 function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
+    # TODO: Should be an inner constructor.
     isclosed(dh) && error("DofHandler already closed")
     # Compute the celltype and make sure all elements have the same one
     CT = getcelltype(dh.grid, first(cellset))
@@ -66,13 +93,6 @@ function _print_field_information(io::IO, mime::MIME"text/plain", sdh::SubDofHan
     end
 end
 
-"""
-    DofHandler(grid::Grid)
-
-Construct a `DofHandler` based on `grid`. Supports:
-- `Grid`s with or without concrete element type (E.g. "mixed" grids with several different element types.)
-- One or several fields, which can live on the whole domain or on subsets of the `Grid`.
-"""
 struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     subdofhandlers::Vector{SubDofHandler{DofHandler{dim, G}}}
     field_names::Vector{Symbol}
@@ -86,6 +106,28 @@ struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     ndofs::ScalarWrapper{Int}
 end
 
+"""
+    DofHandler(grid::Grid)
+
+Construct a `DofHandler` based on the grid `grid`.
+
+After construction any number of discrete fields can be added to the DofHandler using
+[`add!`](@ref). Construction is finalized by calling [`close!`](@ref).
+
+By default fields are added to all elements of the grid. Refer to [`SubDofHandler`](@ref)
+for restricting fields to subdomains of the grid.
+
+# Examples
+
+```julia
+dh = DofHandler(grid)
+ip_u = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u
+ip_p = Lagrange{RefTriangle, 1}()   # scalar interpolation for a field p
+add!(dh, :u, ip_u)
+add!(dh, :p, ip_p)
+close!(dh)
+```
+"""
 function DofHandler(grid::G) where {dim, G <: AbstractGrid{dim}}
     ncells = getncells(grid)
     sdhs = SubDofHandler{DofHandler{dim, G}}[]

src/FEValues/cell_values.jl

@@ -68,12 +68,8 @@ function CellValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, GIP}(qr::QR, ip::I
     dMdξ = fill(zero(dMdξ_t) * T(NaN), n_geom_basefuncs, n_qpoints)
 
     for (qp, ξ) in pairs(getpoints(qr))
-        for basefunc in 1:n_func_basefuncs
-            dNdξ[basefunc, qp], N[basefunc, qp] = shape_gradient_and_value(ip, ξ, basefunc)
-        end
-        for basefunc in 1:n_geom_basefuncs
-            dMdξ[basefunc, qp], M[basefunc, qp] = shape_gradient_and_value(gip, ξ, basefunc)
-        end
+        shape_gradients_and_values!(@view(dNdξ[:, qp]), @view(N[:, qp]), ip, ξ)
+        shape_gradients_and_values!(@view(dMdξ[:, qp]), @view(M[:, qp]), gip, ξ)
     end
 
     detJdV = fill(T(NaN), n_qpoints)

src/FEValues/face_values.jl

@@ -70,12 +70,8 @@ function FaceValues{IP, N_t, dNdx_t, dNdξ_t, T, dMdξ_t, QR, Normal_t, GIP}(qr:
     dMdξ = fill(zero(dMdξ_t) * T(NaN), n_geom_basefuncs, max_n_qpoints, n_faces)
 
     for face in 1:n_faces, (qp, ξ) in pairs(getpoints(qr, face))
-        for basefunc in 1:n_func_basefuncs
-            dNdξ[basefunc, qp, face], N[basefunc, qp, face] = shape_gradient_and_value(ip, ξ, basefunc)
-        end
-        for basefunc in 1:n_geom_basefuncs
-            dMdξ[basefunc, qp, face], M[basefunc, qp, face] = shape_gradient_and_value(gip, ξ, basefunc)
-        end
+        shape_gradients_and_values!(@view(dNdξ[:, qp, face]), @view(N[:, qp, face]), ip, ξ)
+        shape_gradients_and_values!(@view(dMdξ[:, qp, face]), @view(M[:, qp, face]), gip, ξ)
     end
 
     detJdV = fill(T(NaN), max_n_qpoints, n_faces)
@@ -222,8 +218,8 @@ function BCValues(::Type{T}, func_interpol::Interpolation{refshape}, geom_interp
     nqp = zeros(Int,n_boundary_entities)
 
     for n_boundary_entity in 1:n_boundary_entities
-        for (qp, ξ) in enumerate(qrs[n_boundary_entity].points), i in 1:n_geom_basefuncs
-            M[i, qp, n_boundary_entity] = shape_value(geom_interpol, ξ, i)
+        for (qp, ξ) in pairs(qrs[n_boundary_entity].points)
+            shape_values!(@view(M[:, qp, n_boundary_entity]), geom_interpol, ξ)
         end
         nqp[n_boundary_entity] = length(qrs[n_boundary_entity].points)
     end
@@ -261,4 +257,4 @@ function Base.show(io::IO, m::MIME"text/plain", fv::FaceValues)
     print(io, "- Function interpolation: "); show(io, m, fv.func_interp)
     println(io)
     print(io, "- Geometric interpolation: "); show(io, m, fv.geo_interp)
-end
+end

src/Grid/grid.jl

@@ -765,8 +765,11 @@ for INDEX in (:VertexIndex, :EdgeIndex, :FaceIndex)
         #To be able to do a,b = faceidx
         Base.iterate(I::($INDEX), state::Int=1) = (state==3) ?  nothing : (I[state], state+1)
 
-        #For (cellid, faceidx) in faceset
-        Base.in(v::Tuple{Int, Int}, s::Set{$INDEX}) = in($INDEX(v), s)
+        # Necessary to check if, e.g. `(cellid, faceidx) in faceset`
+        Base.isequal(x::$INDEX, y::$INDEX) = x.idx == y.idx
+        Base.isequal(x::Tuple{Int, Int}, y::$INDEX) = x[1] == y.idx[1] && x[2] == y.idx[2]
+        Base.isequal(y::$INDEX, x::Tuple{Int, Int}) = x[1] == y.idx[1] && x[2] == y.idx[2]
+        Base.hash(x::$INDEX, h::UInt) = hash(x.idx, h)
     end
 end
 

src/Grid/topology.jl

@@ -260,17 +260,15 @@ function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, vertexidx::
     cellid, local_vertexid = vertexidx[1], vertexidx[2]
     cell_vertices = vertices(getcells(grid,cellid))
     global_vertexid = cell_vertices[local_vertexid]
-    if include_self
-        vertex_to_cell = top.vertex_to_cell[global_vertexid]
-        self_reference_local = Vector{VertexIndex}(undef,length(vertex_to_cell))
-        for (i,cellid) in enumerate(vertex_to_cell)
-            local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))::Int)
-            self_reference_local[i] = local_vertex
-        end
-        return [top.vertex_vertex_neighbor[global_vertexid].neighbor_info; self_reference_local]
-    else
-        return top.vertex_vertex_neighbor[global_vertexid].neighbor_info
+    vertex_to_cell = top.vertex_to_cell[global_vertexid]
+    self_reference_local = Vector{VertexIndex}()
+    sizehint!(self_reference_local, length(vertex_to_cell))
+    for (i,cellid) in enumerate(vertex_to_cell)
+        local_vertex = VertexIndex(cellid,findfirst(x->x==global_vertexid,vertices(getcells(grid,cellid)))::Int)
+        !include_self && local_vertex == vertexidx && continue
+        push!(self_reference_local, local_vertex)
     end
+    return self_reference_local
 end
 
 function getneighborhood(top::ExclusiveTopology, grid::AbstractGrid{3}, edgeidx::EdgeIndex, include_self=false)

src/PointEval/PointEvalHandler.jl

@@ -125,6 +125,7 @@ function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, globa
     for _ in 1:max_iters
         global_guess = zero(V)
         J = zero(Tensor{2, dim, T})
+        # TODO batched eval after 764 is merged.
         for j in 1:n_basefuncs
             dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j)
             global_guess += N * cell_coordinates[j]

src/PointEval/point_values.jl

@@ -59,9 +59,7 @@ function_symmetric_gradient(pv::PointValues, u::AbstractVector, args...) =
 function reinit!(pv::PointValues, x::AbstractVector{<:Vec{D}}, ξ::Vec{D}) where {D}
     qp = 1 # PointValues only have a single qp
     # TODO: Does M need to be updated too?
-    for i in 1:getnbasefunctions(pv.cv.ip)
-        pv.cv.dNdξ[i, qp], pv.cv.N[i, qp] = shape_gradient_and_value(pv.cv.ip, ξ, i)
-    end
+    shape_gradients_and_values!(@view(pv.cv.dNdξ[:, qp]), @view(pv.cv.N[:, qp]), pv.cv.ip, ξ)
     reinit!(pv.cv, x)
     return nothing
 end
@@ -86,9 +84,6 @@ shape_value(pv::PointValuesInternal, qp::Int, i::Int) = (@assert qp == 1; pv.N[i
 
 # allow on-the-fly updating
 function reinit!(pv::PointValuesInternal{IP}, coord::Vec{dim}) where {dim, shape <: AbstractRefShape{dim}, IP <: Interpolation{shape}}
-    n_func_basefuncs = getnbasefunctions(pv.ip)
-    for i in 1:n_func_basefuncs
-        pv.N[i] = shape_value(pv.ip, coord, i)
-    end
+    shape_values!(pv.N, pv.ip, coord)
     return nothing
 end

src/exports.jl

@@ -17,6 +17,7 @@ export
     DiscontinuousLagrange,
     Serendipity,
     getnbasefunctions,
+    getrefshape,
 
 # Quadrature
     QuadratureRule,