diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 070b10d481..50a7120ad3 100644 --- a/src/Dofs/DofHandler.jl +++ b/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) diff --git a/src/iterators.jl b/src/iterators.jl index 16877ba910..a643158fbb 100644 --- a/src/iterators.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index d03a4e64b4..def659697e 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/test_abstractgrid.jl b/test/test_abstractgrid.jl new file mode 100644 index 0000000000..89b3136ecf --- /dev/null +++ b/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