DofHandler with any <: AbstractGrid

src/Dofs/DofHandler.jl

@@ -9,7 +9,7 @@ Operates slightly faster than [`MixedDofHandler`](@docs). Supports:
 - `Grid`s with a single concrete cell type.
 - One or several fields on the whole domaine.
 """
-struct DofHandler{dim,C,T} <: AbstractDofHandler
+struct DofHandler{dim,T,G<:AbstractGrid{dim}} <: AbstractDofHandler
     field_names::Vector{Symbol}
     field_dims::Vector{Int}
     # TODO: field_interpolations can probably be better typed: We should at least require
@@ -19,11 +19,11 @@ struct DofHandler{dim,C,T} <: AbstractDofHandler
     cell_dofs::Vector{Int}
     cell_dofs_offset::Vector{Int}
     closed::ScalarWrapper{Bool}
-    grid::Grid{dim,C,T}
+    grid::G
     ndofs::ScalarWrapper{Int}
 end
 
-function DofHandler(grid::Grid)
+function DofHandler(grid::AbstractGrid)
     isconcretetype(getcelltype(grid)) || error("Grid includes different celltypes. Use MixedDofHandler instead of DofHandler")
     DofHandler(Symbol[], Int[], Interpolation[], BCValues{Float64}[], Int[], Int[], ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
 end
@@ -302,8 +302,8 @@ function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
     return global_dofs
 end
 
-function cellnodes!(global_nodes::Vector{Int}, grid::Grid{dim,C}, i::Int) where {dim,C}
-    nodes = grid.cells[i].nodes
+function cellnodes!(global_nodes::Vector{Int}, grid::AbstractGrid{dim}, i::Int) where {dim,C}
+    nodes = getcells(grid,i).nodes
     N = length(nodes)
     @assert length(global_nodes) == N
     for j in 1:N
@@ -312,15 +312,16 @@ function cellnodes!(global_nodes::Vector{Int}, grid::Grid{dim,C}, i::Int) where
     return global_nodes
 end
 
-function cellcoords!(global_coords::Vector{Vec{dim,T}}, grid::Grid{dim,C}, i::Int) where {dim,C,T}
-    nodes = grid.cells[i].nodes
+function cellcoords!(global_coords::Vector{Vec{dim,T}}, grid::AbstractGrid{dim}, i::Int) where {dim,C,T}
+    nodes = getcells(grid,i).nodes
     N = length(nodes)
     @assert length(global_coords) == N
     for j in 1:N
-        global_coords[j] = grid.nodes[nodes[j]].x
+        global_coords[j] = getcoordinates(getnodes(grid,nodes[j]))
     end
     return global_coords
 end
+
 cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int) = cellcoords!(global_coords, dh.grid, i)
 
 function celldofs(dh::DofHandler, i::Int)

src/iterators.jl

@@ -41,7 +41,7 @@ struct CellIterator{dim,C,T,DH<:Union{AbstractDofHandler,Nothing}}
     dh::Union{DH,Nothing}
     celldofs::Vector{Int}
 
-    function CellIterator{dim,C,T}(dh::Union{DofHandler{dim,C,T},MixedDofHandler{dim,T,G},Nothing}, cellset::Union{AbstractVector{Int},Nothing}, flags::UpdateFlags) where {dim,C,T,G}
+    function CellIterator{dim,C,T}(dh::Union{DofHandler{dim,T,G},MixedDofHandler{dim,T,G},Nothing}, cellset::Union{AbstractVector{Int},Nothing}, flags::UpdateFlags) where {dim,C,T,G}
         isconcretetype(C) || _check_same_celltype(dh.grid, cellset)
         N = nnodes_per_cell(dh.grid, cellset === nothing ? 1 : first(cellset))
         cell = ScalarWrapper(0)
@@ -64,8 +64,8 @@ end
 
 CellIterator(grid::Grid{dim,C,T}, cellset::Union{AbstractVector{Int},Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) where {dim,C,T} =
     CellIterator{dim,C,T}(grid, cellset, flags)
-CellIterator(dh::DofHandler{dim,C,T}, cellset::Union{AbstractVector{Int},Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) where {dim,C,T} =
-    CellIterator{dim,C,T}(dh, cellset, flags)
+CellIterator(dh::DofHandler{dim,T}, cellset::Union{AbstractVector{Int},Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) where {dim,C,T} =
+    CellIterator{dim,getcelltype(dh.grid),T}(dh, cellset, flags)
 CellIterator(dh::MixedDofHandler{dim,T}, cellset::Union{AbstractVector{Int},Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) where {dim,T} =
     CellIterator{dim,getcelltype(dh.grid),T}(dh, cellset, flags)
 

test/runtests.jl

@@ -16,6 +16,7 @@ include("test_assemble.jl")
 include("test_dofs.jl")
 include("test_constraints.jl")
 include("test_grid_dofhandler_vtk.jl")
+include("test_abstractgrid.jl")
 include("test_mixeddofhandler.jl")
 include("test_l2_projection.jl")
 include("test_pointevaluation.jl")

test/test_abstractgrid.jl

@@ -0,0 +1,91 @@
+@testset "AbstractGrid" begin
+
+    struct SmallGrid{dim,N,C<:Ferrite.AbstractCell} <: Ferrite.AbstractGrid{dim}
+        nodes_test::Vector{NTuple{dim,Float64}}
+        cells_test::NTuple{N,C}
+    end
+
+    Ferrite.getcells(grid::SmallGrid) = grid.cells_test
+    Ferrite.getcells(grid::SmallGrid, v::Union{Int, Vector{Int}}) = grid.cells_test[v]
+    Ferrite.getncells(grid::SmallGrid{dim,N}) where {dim,N} = N
+    Ferrite.getcelltype(grid::SmallGrid) = eltype(grid.cells_test)
+    Ferrite.getcelltype(grid::SmallGrid, i::Int) = typeof(grid.cells_test[i])
+    Ferrite.getcoordinates(x::NTuple{dim,Float64}) where dim = Vec{dim,Float64}(x)
+
+    Ferrite.getnodes(grid::SmallGrid) = grid.nodes_test
+    Ferrite.getnodes(grid::SmallGrid, v::Union{Int, Vector{Int}}) = grid.nodes_test[v]
+    Ferrite.getnnodes(grid::SmallGrid) = length(grid.nodes_test)
+    Ferrite.nnodes_per_cell(grid::SmallGrid, i::Int=1) = Ferrite.nnodes(grid.cells_test[i])
+    Ferrite.n_faces_per_cell(grid::SmallGrid) = nfaces(eltype(grid.cells_test))
+    function Ferrite.getcoordinates!(x::Vector{Vec{dim,T}}, grid::SmallGrid, cell::Int) where {dim,T}
+        for i in 1:length(x)
+            x[i] = Vec{dim,T}(grid.nodes_test[grid.cells_test[cell].nodes[i]])
+        end
+    end
+    function Ferrite.getcoordinates(grid::SmallGrid{dim}, cell::Int) where dim
+        nodeidx = grid.cells_test[cell].nodes
+        return [Vec{dim,Float64}(grid.nodes_test[i]) for i in nodeidx]::Vector{Vec{dim,Float64}}
+    end
+
+    nodes = [(-1.0,-1.0); (0.0,-1.0); (1.0,-1.0); (-1.0,0.0); (0.0,0.0); (1.0,0.0); (-1.0,1.0); (0.0,1.0); (1.0,1.0)]
+    cells = (Quadrilateral((1,2,5,4)), Quadrilateral((2,3,6,5)), Quadrilateral((4,5,8,7)), Quadrilateral((5,6,9,8)))
+    subtype_grid = SmallGrid(nodes,cells)
+    reference_grid = generate_grid(Quadrilateral, (2,2))
+
+    ip = Lagrange{2, RefCube, 1}()
+    qr = QuadratureRule{2, RefCube}(2)
+    cellvalues = CellScalarValues(qr, ip);
+
+    dhs = [DofHandler(grid) for grid in (subtype_grid, reference_grid)]
+    u1 = Vector{Float64}(undef, 9)
+    u2 = Vector{Float64}(undef, 9)
+    ∂Ω = union(getfaceset.((reference_grid, ), ["left", "right", "top", "bottom"])...)
+    dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
+
+    function doassemble!(cellvalues::CellScalarValues{dim}, K::SparseMatrixCSC, dh::DofHandler) where {dim}
+        n_basefuncs = getnbasefunctions(cellvalues)
+        Ke = zeros(n_basefuncs, n_basefuncs)
+        fe = zeros(n_basefuncs)
+        f = zeros(ndofs(dh))
+        assembler = start_assemble(K, f)
+        for cellid in 1:getncells(dh.grid)
+            fill!(Ke, 0)
+            fill!(fe, 0)
+            coords = getcoordinates(dh.grid, cellid)
+            reinit!(cellvalues, coords)
+            for q_point in 1:getnquadpoints(cellvalues)
+                dΩ = getdetJdV(cellvalues, q_point)
+                for i in 1:n_basefuncs
+                    v  = shape_value(cellvalues, q_point, i)
+                    ∇v = shape_gradient(cellvalues, q_point, i)
+                    fe[i] += v * dΩ
+                    for j in 1:n_basefuncs
+                        ∇u = shape_gradient(cellvalues, q_point, j)
+                        Ke[i, j] += (∇v ⋅ ∇u) * dΩ
+                    end
+                end
+            end
+            assemble!(assembler, celldofs(dh,cellid), fe, Ke)
+        end
+        return K, f
+    end
+
+    for (dh,u) in zip(dhs,(u1,u2))
+        push!(dh, :u, 1)
+        close!(dh)
+        ch = ConstraintHandler(dh)
+        add!(ch, dbc)
+        close!(ch)
+        update!(ch, 0.0)
+        K = create_sparsity_pattern(dh);
+        K, f = doassemble!(cellvalues, K, dh);
+        apply!(K, f, ch)
+        sol = K \ f
+        u .= sol
+    end
+
+    @test Ferrite.ndofs_per_cell(dhs[1]) == Ferrite.ndofs_per_cell(dhs[2])
+    @test Ferrite.celldofs(dhs[1],3) == Ferrite.celldofs(dhs[2],3)
+    @test Ferrite.ndofs(dhs[1]) == Ferrite.ndofs(dhs[2])
+    @test isapprox(u1,u2,atol=1e-8)
+end