add MixedGrid

src/Dofs/DofHandler.jl

@@ -312,15 +312,22 @@ 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
+# shared implementation for all grids
+function cellcoords!(global_coords::Vector{Vec{dim,T}}, grid::AbstractGrid, cell::C) where {dim,C<:Ferrite.AbstractCell,T}
+    nodes = cell.nodes
     N = length(nodes)
     @assert length(global_coords) == N
     for j in 1:N
         global_coords[j] = grid.nodes[nodes[j]].x
     end
     return global_coords
 end
+
+function cellcoords!(global_coords::Vector{Vec{dim,T}}, grid::AbstractGrid, i) where {dim,T}
+    cell = getcells(grid, i)
+    cellcoords!(global_coords, grid, cell)
+end
+
 cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int) = cellcoords!(global_coords, dh.grid, i)
 
 function celldofs(dh::DofHandler, i::Int)

src/Dofs/MixedDofHandler.jl

@@ -23,9 +23,9 @@ A `FieldHandler` must fullfill the following requirements:
 
 Notice that a `FieldHandler` can hold several fields.
 """
-mutable struct FieldHandler
+mutable struct FieldHandler{CI}
     fields::Vector{Field}
-    cellset::Set{Int}
+    cellset::Set{CI}
 end
 
 struct CellVector{T}
@@ -46,18 +46,16 @@ Construct a `MixedDofHandler` 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 MixedDofHandler{dim,T,G<:AbstractGrid{dim}} <: AbstractDofHandler
+struct MixedDofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     fieldhandlers::Vector{FieldHandler}
     cell_dofs::CellVector{Int}
-    cell_nodes::CellVector{Int}
-    cell_coords::CellVector{Vec{dim,T}}
     closed::ScalarWrapper{Bool}
     grid::G
     ndofs::ScalarWrapper{Int}
 end
 
-function MixedDofHandler(grid::Grid{dim,C,T}) where {dim,C,T}
-    MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], CellVector(Int[],Int[],Int[]), CellVector(Int[],Int[],Int[]), CellVector(Vec{dim,T}[],Int[],Int[]), Ferrite.ScalarWrapper(false), grid, Ferrite.ScalarWrapper(-1))
+function MixedDofHandler(grid::G) where {dim, G<:AbstractGrid{dim}}
+    return MixedDofHandler{dim,G}(FieldHandler[], CellVector(Int[],Int[],Int[]), ScalarWrapper(false), grid, ScalarWrapper(-1))
 end
 
 getfieldnames(fh::FieldHandler) = [field.name for field in fh.fields]
@@ -73,26 +71,14 @@ function celldofs!(global_dofs::Vector{Int}, dh::MixedDofHandler, i::Int)
     unsafe_copyto!(global_dofs, 1, dh.cell_dofs.values, dh.cell_dofs.offset[i], length(global_dofs))
     return global_dofs
 end
+celldofs!(global_dofs::Vector{Int}, dh::MixedDofHandler, i::CellId) = celldofs!(global_dofs, dh, globalid(dh.grid, i))
 
 function celldofs(dh::MixedDofHandler, i::Int)
     @assert isclosed(dh)
     return dh.cell_dofs[i]
 end
 
-function cellcoords!(global_coords::Vector{Vec{dim,T}}, dh::MixedDofHandler, i::Int) where {dim,T}
-    @assert isclosed(dh)
-    @assert length(global_coords) == nnodes_per_cell(dh, i)
-    unsafe_copyto!(global_coords, 1, dh.cell_coords.values, dh.cell_coords.offset[i], length(global_coords))
-    return global_coords
-end
-
-function cellnodes!(global_nodes::Vector{Int}, dh::MixedDofHandler, i::Int)
-    @assert isclosed(dh)
-    @assert length(global_nodes) == nnodes_per_cell(dh, i)
-    unsafe_copyto!(global_nodes, 1, dh.cell_nodes.values, dh.cell_nodes.offset[i], length(global_nodes))
-    return global_nodes
-end
-
+cellcoords!(global_coords::Vector{<:Vec}, dh::MixedDofHandler, i) = cellcoords!(global_coords, dh.grid, i)
 
 """
     getfieldnames(dh::MixedDofHandler)
@@ -140,12 +126,12 @@ Add all fields of the [`FieldHandler`](@ref) `fh` to `dh`.
 """
 function Base.push!(dh::MixedDofHandler, fh::FieldHandler)
     @assert !isclosed(dh)
-    _check_same_celltype(dh.grid, collect(fh.cellset))
+    # _check_same_celltype(dh.grid, collect(fh.cellset))
     _check_cellset_intersections(dh, fh)
     # the field interpolations should have the same refshape as the cells they are applied to
     refshapes_fh = getrefshape.(getfieldinterpolations(fh))
     # extract the celltype from the first cell as the celltypes are all equal
-    cell_type = typeof(dh.grid.cells[first(fh.cellset)])
+    cell_type = typeof(getcells(dh.grid, first(fh.cellset)))
     refshape_cellset = getrefshape(default_interpolation(cell_type))
     for refshape in refshapes_fh
         refshape_cellset == refshape || error("The RefShapes of the fieldhandlers interpolations must correspond to the RefShape of the cells it is applied to.")
@@ -227,7 +213,8 @@ function __close!(dh::MixedDofHandler{dim}) where {dim}
     @debug "\n\nCreating dofs\n"
     for fh in dh.fieldhandlers
         # sort the cellset since we want to loop through the cells in a fixed order
-        cellnumbers = sort(collect(fh.cellset))
+        # now only works for MixedGrid
+        cellnumbers = get_sorted_cellset(fh)
         nextdof = _close!(
             dh,
             cellnumbers,
@@ -244,24 +231,14 @@ function __close!(dh::MixedDofHandler{dim}) where {dim}
     dh.ndofs[] = maximum(dh.cell_dofs.values)
     dh.closed[] = true
 
-    #Create cell_nodes and cell_coords (similar to cell_dofs)
-    push!(dh.cell_nodes.offset, 1)
-    push!(dh.cell_coords.offset, 1)
-    for cell in dh.grid.cells
-        for nodeid in cell.nodes
-            push!(dh.cell_nodes.values, nodeid)
-            push!(dh.cell_coords.values, dh.grid.nodes[nodeid].x)
-        end
-        push!(dh.cell_nodes.offset, length(dh.cell_nodes.values)+1)
-        push!(dh.cell_coords.offset, length(dh.cell_coords.values)+1)
-        push!(dh.cell_nodes.length, length(cell.nodes))
-        push!(dh.cell_coords.length, length(cell.nodes))
-    end
-
     return dh, vertexdicts, edgedicts, facedicts
 
 end
 
+get_sorted_cellset(fh::FieldHandler{Int}) = sort!(collect(fh.cellset))
+get_sorted_cellset(fh::FieldHandler{CellId{I}}) where I = sort!(collect(fh.cellset), by=c->c.i)
+globalid(::AbstractGrid, i::Int) = i # fallback for regular grids
+
 function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, field_names, field_dims, field_interpolations, nextdof, vertexdicts, edgedicts, facedicts, celldicts) where {dim}
 
     ip_infos = Ferrite.InterpolationInfo[]
@@ -277,9 +254,10 @@ function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fiel
 
     # loop over all the cells, and distribute dofs for all the fields
     for ci in cellnumbers
-        dh.cell_dofs.offset[ci] = length(dh.cell_dofs.values)+1
+        gi = globalid(dh.grid, ci)
+        dh.cell_dofs.offset[gi] = length(dh.cell_dofs.values)+1
 
-        cell = dh.grid.cells[ci]
+        cell = getcells(dh.grid, ci)
         cell_dofs = Int[]  # list of global dofs for each cell
         @debug "Creating dofs for cell #$ci"
 
@@ -301,13 +279,13 @@ function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fiel
             end
 
             if ip_info.ncelldofs > 0
-                nextdof = add_cell_dofs(cell_dofs, ci, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof)
+                nextdof = add_cell_dofs(cell_dofs, gi, celldicts[fi], field_dims[local_num], ip_info.ncelldofs, nextdof)
             end
 
         end
         # after done creating dofs for the cell, push them to the global list
         push!(dh.cell_dofs.values, cell_dofs...)
-        dh.cell_dofs.length[ci] = length(cell_dofs)
+        dh.cell_dofs.length[gi] = length(cell_dofs)
 
         @debug "Dofs for cell #$ci:\n\t$cell_dofs"
     end # cell loop

src/Ferrite.jl

@@ -51,6 +51,7 @@ include("FEValues/face_integrals.jl")
 include("Grid/grid.jl")
 include("Grid/grid_generators.jl")
 include("Grid/coloring.jl")
+include("Grid/mixed_grid.jl")
 
 # Dofs
 include("Dofs/DofHandler.jl")

src/Grid/mixed_grid.jl

@@ -0,0 +1,51 @@
+mutable struct MixedGrid{dim, C, T<:Real} <: Ferrite.AbstractGrid{dim}
+    cells::C # Tuple of concretely typed cell vectors
+    ncells_per_vector::Vector{Int}
+    nodes::Vector{Node{dim,T}}
+    # Sets
+    cellsets::Dict{String,Set{Int}}
+    nodesets::Dict{String,Set{Int}}
+    facesets::Dict{String,Set{FaceIndex}} 
+    edgesets::Dict{String,Set{EdgeIndex}} 
+    vertexsets::Dict{String,Set{VertexIndex}} 
+    # Boundary matrix (faces per cell × cell)
+    boundary_matrix::SparseMatrixCSC{Bool,Int}
+end
+
+function MixedGrid(cells::C,
+        nodes::Vector{Node{dim,T}};
+        cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
+        nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
+        facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
+        edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
+        vertexsets::Dict{String,Set{VertexIndex}}=Dict{String,Set{VertexIndex}}(),
+        boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
+        ncells_per_type = collect(length.(cells))
+    return MixedGrid(cells, ncells_per_type, nodes, cellsets, nodesets, facesets, edgesets, vertexsets, boundary_matrix)
+end
+
+struct CellId{I}
+    i::Int
+end
+
+function globalid(grid::MixedGrid, cellid::Ferrite.CellId{I}) where I
+    global_id = 0
+    for i=1:I-1
+        global_id += grid.ncells_per_vector[I]
+    end
+    global_id += cellid.i
+    return global_id
+end
+
+Ferrite.getcells(grid::MixedGrid, cellid::Ferrite.CellId{I}) where I = grid.cells[I][cellid.i]
+Ferrite.getncells(grid::MixedGrid) = sum(grid.ncells_per_vector)
+
+# Inherently type unstable
+function local_index(grid::MixedGrid, global_id::Int)
+    local_idx = global_id
+    for (i, ncells) in enumerate(grid.ncells_per_vector) 
+        local_idx - ncells < 1 && return NewCellIndex{i}(local_idx)
+        local_idx -= ncells
+    end
+    error("Local index corresponding to global_id=$global_id not found.")
+end

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,C,T},MixedDofHandler{dim,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)