RFC: Introducing a `MixedGrid` #445

kimauth · 2022-04-10T10:04:29Z

Following up on our discussion from last week, here comes a suggestion on how to handle mixed grids in a more data-oriented manner.

The main gripes about the current implementation are:

Grids with several celltypes store their cells as Vector{Union{CellType1, CelType2,...}} or even worse as Vector{AbstractCell}. This makes extracting a cell based on its global id inherently type unstable.
Due to the above reason, the MixedDofHandler stores its own grid representation so that it can efficiently access cell nodes and cell coordinates.
DofHandler & MixedDofHandler. Mixed Celltypes in a grid are now reflected by need to use MixedDofHandler. It would be nice to reflect this in the Grid datatype instead and use a single DofHandler no matter what the grid type. (Currently the DofHandler additionally does not support fields on subsets of the domain, while the MixedDofHandler does.)
When assembling with a MixedDofHandler, it is crucial for performance to not iterate over all cells at once. Instead a function barrier should be set for each FieldHandler. This is not mentioned in our docs.

This PR attempts to solve the first two points and lay ground for solving the 3rd in the future.

MixedGrid type

It introduces a MixedGrid that stores its cells as Tuple{Vector{CellType1}, Vector{CellType2}, ...}.
The cells in the grid are indexed by

struct CellId{I}
    i::Int
end

where I is the index within the tuple and i is the index within the Vector. This can be seen as local cell id. Cells also have global cell ids (Ints), which count through the cells tuple from beginning to end. Global ids can thus easily computed from local ids.

This means that cells will likely be renumbered when importing a mixed grid from an external grid generator. I don't see this as a huge problem, but an alternative approach could be to store a mapping local id -> global id.

MixedDofHandler + MixedGrid

A grid with mixed celltypes should always be represented as MixedGrid and always be indexed by CellId. That is, the FieldHandlers then store a (concretely typed!) Set{CellId{I}} and all cellsets + facesets in the mixed grid should reference cells by CellId instead of their global id (the latter isn't implemented yet for this PR).

The MixedDofHandler does not store cell nodes and cell coordinates anymore. Methods for extracting them fall back to the implementations for the underlying grid.

The MixedDofHandler allows to use either Grid or MixedGrid with it. Therefore, the celldofs are still stored by global cell ids.

DofHandler methods + Benchmarking

Here is a list of all functions that use the DofHandlers (to get an overview). Some methods are defined for AbstractDofHandler, but will only work for the DofHandler, thus they're listed on the left.

DofHandler	MixedDofHandler
`getfielddim(dh::DofHandler, field_idx::Int)`	no exact equivalent
`getfielddim(dh::DofHandler, name::Symbol)`	`getfielddim(dh::MixedDofHandler, name::Symbol)`
`ndim(dh::AbstractDofHandler, field_name::Symbol)`, same functionality as getfielddim	doesn't exist
`nfields(dh::AbstractDofHandler)`, see #444	`nfields(dh::MixedDofHandler)`
`ndofs_per_cell(dh::AbstractDofHandler, cell::Int=1)`	`ndofs_per_cell(dh::MixedDofHandler, cell::Int=1)`
`getfieldnames(dh::AbstractDofHandler)`	`getfieldnames(dh::MixedDofHandler)`, `getfieldnames(fh::FieldHandler)`

`find_field(dh::DofHandler, field_name::Symbol)`	`find_field(fh::FieldHandler, field_name::Symbol)`
`field_offset(dh::DofHandler, field_name::Symbol)`	`field_offset(fh::FieldHandler, field_name::Symbol)`
`dof_range(dh::DofHandler, field_name::Symbol)`	`dof_range(fh::FieldHandler, field_name::Symbol)`
`getfieldinterpolation(dh::DofHandler, field_idx::Int)`	no exact equivalent
no excact equivalent	`getfieldinterpolations(fh::FieldHandler)`

`cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int)`	`cellcoords!(global_coords::Vector{Vec{dim,T}}, dh::MixedDofHandler, i::Int) where {dim,T}`
`celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)`	`celldofs!(global_dofs::Vector{Int}, dh::MixedDofHandler, i::Int)`
`celldofs(dh::DofHandler, i::Int)`	`celldofs(dh::MixedDofHandler, i::Int)`
use grid version	`cellnodes!(global_nodes::Vector{Int}, dh::MixedDofHandler, i::Int)`

Most of these are used for dof distribution, I'd expect that only celldofs! and cellcoords! are used in performance relevant code, so let's focus on these for some micro benchmarks. The benchmarks extract cellcoords and celldofs for a Quadrilateral in a mixed grid with Triangles and Quadrilaterals.

	`DofHandler`	`MixedDofHandler#master` + mixed grid	`MixedDofHandler#thisPR` + mixed grid	`MixedDofHandler#thisPR` + concrete grid
`cellcoords!`	4.100 ns (0 allocations: 0 bytes)	9.300 ns (0 allocations: 0 bytes)	4.300 ns (0 allocations: 0 bytes)	4.100 ns (0 allocations: 0 bytes)
`celldofs!`	9.500 ns (0 allocations: 0 bytes)	9.810 ns (0 allocations: 0 bytes)	9.300 ns (0 allocations: 0 bytes)	9.300 ns (0 allocations: 0 bytes)

For this case, the MixedDofHandler now performs equally well as the DofHandler, no matter if it is used with Grid or with MixedGrid.

Technically, it was possible to use the grid version of cellcoords! before this PR (and thus get rid of the grid representation in the MixedDofHandler), however that can be (very) slow:

cell type	`Quadrilateral` (concretly typed)	`Union{Triangle, Quadrilateral}`	`AbstractCell`
`cellcoords!` `(xe, grid, i)`	3.700 ns (0 allocations: 0 bytes)	6.700 ns (0 allocations: 0 bytes)	1.710 μs (22 allocations: 720 bytes)

Usage

The usage pattern remains exactly the same as before.

# e.g. for the quads and tris example above:
buffers_quads = (xe_quads, dofs_quads, cv_quads, ke_quads, re_quads)
buffers_tris = (xe_tris, dofs_tris, cv_tris, ke_tris, re_tris)
buffers = [buffers_quads, buffers_tris]

for (fh_idx, fh) in enumerate(dh.fieldhandlers)
    xe, dofs, cv, ke, re = buffers[fh_idx]
    loop_fh(dh, fh, ke, re, xe, dofs, cv)
end

function loop_fh(dh, fh, ke, re, xe, dofs, cv)
    for cellid in fh.cellset
        cellcoords!(xe, dh, cellid)
        celldofs!(dofs, dh, cellid)
        # ... do stuff like element routine
        # @views ue = u[dofs] 
        # assemble_cell!(ke, re, xe, cv, ue)
    end
end

koehlerson

I like this PR; included a small suggestion. Would be a nice working example of using the grid interface. We should take care that functionality doesn't diverge too much. DofHandler cannot use any other grid than Grid. #379 tried to address it, but seems to be stalled. So with including something like this, it would be quite unnatural for new or not so experienced users to have a DofHandler who can't take the MixedGrid (which isn't probably that useful). Therefore, a long term strategy would be nice how this problem will be resolved.

koehlerson · 2022-04-21T09:44:02Z

src/Dofs/DofHandler.jl

-    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


There is an interface function which should be used, because the functionality of cellcoords! can be recovered by the already defined interface functions. I think it would be worth considering to put this function in another file, e.g. in the grid.jl file, doesn't seem to be exclusively related the dofhandler

Suggested change

global_coords[j] = grid.nodes[nodes[j]].x

global_coords[j] = getnodes(grid,nodes[j]).x

Using the interface function is definitely the right way to go.

I also agree that this should live in the grid file. I simply didn't move it yet so that it's easier to see the actual diff caused by this PR (and not the one by shuffling around functions).

koehlerson · 2022-04-21T09:50:23Z

src/Grid/mixed_grid.jl

+    edgesets::Dict{String,Set{EdgeIndex}} 
+    vertexsets::Dict{String,Set{VertexIndex}} 
+    # Boundary matrix (faces per cell × cell)
+    boundary_matrix::SparseMatrixCSC{Bool,Int}


@termi-official is there a better way to encode the boundary information into the grid? Would be a good opportunity to include it in this struct

Yes and no. It really depends on your use case. You can get away without using the boundary_matrix and I honestly think that it should not be in any grid struct, but in a separate (boundary information) data structure.

I mean we have the facesets data structure anyway to distinguish different portions of the boundary - and a set lookup should be O(1). If we want to check if we are on the interior, then the upcoming topology PR should handle this better.

Maybe @fredrikekre does know more.

I'd be open to removing this field - I don't use it myself. The reason that it's still there is to stay as close to existing code as possible.

Note that removing the field is definitely a breaking change, so I do not think we can remove it immediately. (See

Ferrite.jl/src/iterators.jl

Line 90 in ea716bc

@inline onboundary(ci::CellIterator, face::Int) = ci.grid.boundary_matrix[face, ci.current_cellid[]]

which is exported.) This is also why I asked Frederik for a comment.

I guess the reason this exists is that it is faster than the set-lookup:

julia> grid = generate_grid(Triangle, (20, 20)); julia> ci = CellIterator(grid); reinit!(ci, 1); julia> @btime onboundary($ci, $1); 4.300 ns (0 allocations: 0 bytes) julia> faceindex = FaceIndex(1,1) FaceIndex((1, 1)) julia> @btime in($faceindex, $(grid.facesets["bottom"])) 16.032 ns (0 allocations: 0 bytes)

So from that perspective, I'd suggest leaving it as it is for now and eventually tackle if it should be in the grid or not in another PR?

It might make more sense to have it be a lookup for whether the cell has any face on a boundary. Then you can check once per cell, instead of once epr face. (This also works better with multi-celltype grids).

lijas

Seems like a good PR to get a type-stable grid with mixed celltypes, and getting rid of the grid-representation in the MixedDofhandler. The MixedGrid is quite minimalistic, and fits into the grid api without changing other parts of the code too much.

Just one thought:

It seems like just type-annotating the normal Grid (when it includes mixed celltypes) gives similar benchmarks when calling cellcords!(...), as compared to with your PR, e.g.:


function f1(coords, grid::Grid{2,Ferrite.AbstractCell,Float64}, cellset::Set{Int})
    for i in cellset
        cell::Quadrilateral = grid.cells[i]; #type annotate
        Ferrite.cellcoords!(coords, grid, cell)
        #...
    end
end

I think it is typically the case that one can type-annotate their cell-loops (that is at least the case in my code), so It is worth considering if we want want to add a completely new Grid-type, if we dont gain that much?
Personally I would probably use the type-annotation approach above, but this PR supports that as well, so that is nice.

lijas · 2022-05-01T18:25:09Z

src/Grid/mixed_grid.jl

+
+function MixedGrid(cells::C,
+        nodes::Vector{Node{dim,T}};
+        cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),


I guess cellsets should not be a vector of Ints in the MixedGrid?

No, working on that :)

lijas · 2022-05-01T18:42:45Z

src/Grid/mixed_grid.jl

+end
+
+struct CellId{I}
+    i::Int


Having the cell id as a type parameter is clever, but can give some trouble if one has a cellset with different celltypes: https://docs.julialang.org/en/v1/manual/performance-tips/#The-dangers-of-abusing-multiple-dispatch-(aka,-more-on-types-with-values-as-parameters)

I dont think that is necessarily that common though...

It's actually inherently needed here, so that we can know the celltype in a type stable manner. Of course that means that collections of cells with different I should be avoided. That is where this PR is coming from though, the purpose it to store the cells of the same type together.
It however does become a little tricky especially with facesets. I'll comment on that when I've pushed the changes on the different sets.

lijas · 2022-05-02T06:10:23Z

src/Dofs/MixedDofHandler.jl

    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)


Define > and < for a CellId instead?

KnutAM · 2022-05-03T08:01:57Z

src/iterators.jl

@@ -41,7 +41,7 @@ struct CellIterator{dim,C,T,DH<:Union{AbstractDofHandler,Nothing}}
    dh::Union{DH,Nothing}
    celldofs::Vector{Int}



Seems like the CellIterator will still have an Abstract celltype?
This PR doesn't seem necessary for fixing that but would enable lifting the type instability higher up (i.e. making the creation of CellIterator type-stable also) !

Haven't done anything on the CellIterator yet, I'd hope fixing the CellIterator for mixed grids would come naturally once we agree on a grid format.

KnutAM · 2022-06-04T16:10:37Z

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}


Is ncells_per_vector needed?
Did some quick benchmarks for getncells and couldn't see that it would be faster (actually slower but so fast that it wouldn't matter anyway)

termi-official · 2022-06-05T15:59:43Z

I see that this PR has quite some impact on the microbenchmarks, but does it really have a significant impact on overall assembly performance for non-trivial assembly loops? I mean first, we cut down the cost of calls which are usually not the most costly ones by several nanoseconds. Second, this PR feels a bit contrary to the goal of merging the DofHandler and MixedDofhandler, since we diversify the code base further. From a conceptual point, f we want to merge this code, wouldn't it be more beneficial if the MixedGrid replaces the current Grid? This also cuts maintainance cost.

fredrikekre · 2023-05-20T23:42:39Z

Superseded by the DofHandler merge.

add MixedGrid

81440fa

KnutAM mentioned this pull request Apr 14, 2022

SubDofHandler/MixedGrid/PorousMedia example #448

Merged

koehlerson reviewed Apr 21, 2022

View reviewed changes

koehlerson mentioned this pull request Apr 22, 2022

DofHandler with any <: AbstractGrid #379

Merged

lijas reviewed May 1, 2022

View reviewed changes

lijas reviewed May 2, 2022

View reviewed changes

KnutAM reviewed May 3, 2022

View reviewed changes

KnutAM reviewed Jun 4, 2022

View reviewed changes

kimauth mentioned this pull request Oct 16, 2022

Decoupling of Geometrical, Topological and Algebraical Interfaces #519

Closed

7 tasks

fredrikekre closed this May 20, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RFC: Introducing a `MixedGrid` #445

RFC: Introducing a `MixedGrid` #445

kimauth commented Apr 10, 2022

koehlerson left a comment

koehlerson Apr 21, 2022

kimauth Apr 21, 2022

koehlerson Apr 21, 2022

termi-official Apr 21, 2022

termi-official Apr 21, 2022

termi-official Apr 21, 2022

kimauth Apr 21, 2022

termi-official Apr 21, 2022

kimauth Apr 22, 2022

fredrikekre May 3, 2022

lijas left a comment

lijas May 1, 2022

kimauth Jun 3, 2022

lijas May 1, 2022

kimauth Jun 3, 2022

lijas May 2, 2022

KnutAM May 3, 2022

kimauth Jun 3, 2022

KnutAM Jun 4, 2022

termi-official commented Jun 5, 2022

fredrikekre commented May 20, 2023

	global_coords[j] = grid.nodes[nodes[j]].x
	global_coords[j] = getnodes(grid,nodes[j]).x

		@@ -41,7 +41,7 @@ struct CellIterator{dim,C,T,DH<:Union{AbstractDofHandler,Nothing}}
		dh::Union{DH,Nothing}
		celldofs::Vector{Int}

RFC: Introducing a MixedGrid #445

RFC: Introducing a MixedGrid #445

Conversation

kimauth commented Apr 10, 2022

MixedGrid type

MixedDofHandler + MixedGrid

DofHandler methods + Benchmarking

Usage

koehlerson left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lijas left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

termi-official commented Jun 5, 2022

fredrikekre commented May 20, 2023

RFC: Introducing a `MixedGrid` #445

RFC: Introducing a `MixedGrid` #445