CellIterator iterates over cells with no DOFs

[Joroks]

I've tried out the fix for #966. However, it seems like the CellIterator still iterates over every cell in the Grid instead of only the cells belonging to the DofHandler. When Iterating over a cell that doesn't belong to the DofHandler, the CellCache will be reinitialized with an empty dof vector.

using Ferrite

grid = generate_grid(Hexahedron, (3,3,3))

dh = DofHandler(grid)
sdh = SubDofHandler(dh, [1,2,3])
add!(sdh, :u, Lagrange{RefHexahedron,1}()^3)
close!(dh)

for cell in CellIterator(sdh) # works
    @assert celldofs(cell) |> !isempty
end

for cell in CellIterator(dh, [1,2,3]) # works
    @assert celldofs(cell) |> !isempty
end

for cell in CellIterator(dh) # errors
    @assert celldofs(cell) |> !isempty
end

[Joroks]

I've looked at the changes made in #967 and this behavior seems intentional:

function test_celliterator_on_true_subdomain_smoketest()
    grid = generate_grid(Hexahedron, (2,2,2))

    dh = DofHandler(grid)
    sdh = SubDofHandler(dh, [1,2,3])
    ip = Lagrange{RefHexahedron,1}()
    add!(sdh, :u, ip)
    close!(dh)

    # The following statements just check that the iterator
    # does not crash at least. Regression for #966
    for cell in CellIterator(sdh)
    end

    for cell in CellIterator(dh, [1,2,3])
    end

    for cell in CellIterator(dh)
        if cellid(cell) <= 3
            @test length(celldofs(cell)) == getnbasefunctions(ip)
        else
            @test length(celldofs(cell)) == 0
        end
    end
end

However, I think it would be better to just skip iterating over those cells entirely.

[termi-official]

Note that it is not recommended to iterate over a dof handler with multiple subdofhandlers. Do you have a specific use-case where this is necessary?

[KnutAM]

I think both the interpretations

1. "Iterate over the cells included in the DofHandler"
2. "Iterate over the cells in the DofHandler's grid"

can be motivated. I think @DRollin also was using DofHandlers which didn't include all cells, so pulling you into this discussion.

In case of interpretation 1, I would then rather say we should have "Iterate over each subdofhandler in the DofHandler" (the result is the same, just the order of iteration might differ). But as @termi-official, normally when defining subdofhandlers, I would expect one would loop over each subdofhandler. So a third option could be,

3. Error when trying to loop over all cells in a DofHandler whose subdofhandlers don't include all grids in the cell.

And then instruct users to either call CellIterator(dh, cellset) or CellIterator(sdh)

[Joroks]

I'm using multiple DofHandlers with a single SubDofHandler each, as they each field utilizes only part of the Grid. I considered using a single DofHandler for both fields, but I don't think that this is possible with how the DofHandler works at the moment.

I'm assembling each Matrix-Block of my coupled system seperately and then concatenate them to get my final Matrix. I therefore need independent DOF counting for each of my fields.

I'm also using a time-dependent non linear model for the FE part of the Simulation, so I have to reassemble the FE part of the system matrix repeatedly while all other parts of the system matrix stay unchanged throughout the entirety of the simulation.

using Ferrite

grid = generate_grid(Quadrilateral, (5,5))

fe = [1, 2, 3, 4, 5, 6, 10, 11, 15, 16, 20, 21, 22, 23, 24, 25] #cellset pure FE-domain
bd = [7, 8, 9, 12, 14, 17, 18, 19] # cellset FE-MD boundary-domain
md = [13] # cellset pure MD-domain

## DofHandler FE domain (fe+bd)
FE = DofHandler(grid)

fe_sub_handler = SubDofHandler(FE, fe ∪ bd)
add!(fe_sub_handler, :u, Lagrange{RefQuadrilateral, 1}()^2)

close!(FE)

## DofHandler for Lagrange Multipliers (bd only)
LM = DofHandler(grid)

lm_sub_handler = SubDofHandler(LM, bd)
add!(lm_sub_handler, :l, Lagrange{RefQuadrilateral, 1}()^2)

close!(LM)

## For the MD-Domain (md+bd) I'm doing something else enirely.

I don't think it's unreasonable to only want to use part of the Grid that you have loaded in in the first place - Maybe It contains more Information than you actually need. Maybe we could change the signature for the DofHandler constructor, such that the user can specify, which part of the grid he wants to use?

This could then look something like this:

## DofHandler FE domain (fe+bd)
FE = DofHandler(grid, fe ∪ bd)
add!(FE, :u, Lagrange{RefQuadrilateral, 1}()^2)

close!(FE)

## DofHandler for Lagrange Multipliers (bd only)
LM = DofHandler(grid, bd)
add!(LM, :l, Lagrange{RefQuadrilateral, 1}()^2)

close!(LM)

With this It would also be possible to support both interpretations:

## create a DofHandler that iterates over all cells in the grid but only has SubDofHandlers defined for a true subset.
FE = DofHandler(grid) # would be equivalent to DofHandler(grid, 1:getncells(grid))

fe_sub_handler = SubDofHandler(FE, fe ∪ bd)
add!(fe_sub_handler, :u, Lagrange{RefQuadrilateral, 1}()^2)

close!(FE)

[Joroks]

Full Grid:

+----+----+----+----+----+
|  1 |  2 |  3 |  4 |  5 |
+----+----+----+----+----+
|  6 |  7 |  8 |  9 | 10 |
+----+----+----+----+----+
| 11 | 12 | 13 | 14 | 15 |
+----+----+----+----+----+
| 16 | 17 | 18 | 19 | 20 |
+----+----+----+----+----+
| 21 | 22 | 23 | 24 | 25 |
+----+----+----+----+----+

FE:

+----+----+----+----+----+
|  1 |  2 |  3 |  4 |  5 |
+----+----+----+----+----+
|  6 |    |    |    | 10 |
+----+----+----+----+----+
| 11 |    |    |    | 15 |
+----+----+----+----+----+
| 16 |    |    |    | 20 |
+----+----+----+----+----+
| 21 | 22 | 23 | 24 | 25 |
+----+----+----+----+----+

BD:

+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+
|    |  7 |  8 |  9 |    |
+----+----+----+----+----+
|    | 12 |    | 14 |    |
+----+----+----+----+----+
|    | 17 | 18 | 19 |    |
+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+

MD:

+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+
|    |    | 13 |    |    |
+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+
|    |    |    |    |    |
+----+----+----+----+----+

[termi-official]

Thanks for providing the information. Sorry I was not precise in my question tho. The part I am not understanding yet is why you need

    for cell in CellIterator(FE)
         # do stuff
    end
    for cell in CellIterator(LM)
         # do stuff
    end

The intended way would be to just use the subdofhandlers

    for cell in CellIterator(fe_sub_handler)
         # do stuff
    end
    for cell in CellIterator(lm_sub_handler)
         # do stuff
    end

which iterates over the subdomain.

[Joroks]

I'd say that this is mostly a matter of convenience and ease of use. I also think that the way the CellIterator currently works is a bit unintuitive.

SubDofHandlers are intended to be used to define multiple element types and/or define different fields for seperate regions in the full domain. If the user doesn't need multiple element types or different fields for different regions, they generally don't have to bother with the SubDofHandlers as Ferrite provides convenience methods for this simplified usecase:

dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefQuadrilateral, 1}()^2)
close!()

for cell in CellIterator(dh)
    # do stuff
end

This is effectively equivalent to

dh = DofHandler(grid)
sdh = SubDofHandler(dh, 1:getncells(grid)
add!(sdh, :u, Lagrange{RefQuadrilateral, 1}()^2)
close!()

# somewhere else in the program; sdh is in a different scope

for cell in CellIterator(only(dh.subdofhandlers))
    # do stuff
end

You are right, that I could use the SubDofHandlers to achive what I want to do, but it's not obvious that that's what I have to do if I only want my DofHandler to operate on part of my grid. (With otherwise a single field and a single element type)

Furthermore It is not obvious that after I have created my DofHandler with a single SubDofHandler, that the CellIterator will still iterate over every cell In the grid - even those that don't have a SubDofHandler associated with them. It's easy to do

dh = DofHandler(grid)
sdh = SubDofHandler(dh, 1:getncells(grid)
add!(sdh, :u, Lagrange{RefQuadrilateral, 1}()^2)
close!()

for cell in CellIterator(dh)
    # do stuff
end

assuming that it should work.

I think the main problem is that up until now it was assumed that when you create a DofHandler, that you also want it to operate on the entire grid. This is simply not always the case.

DofHandler would be more versitile and easy to use if you could explicitly specify what part of the grid you want it to operate on:

function DofHandler(grid, cells=nothing)
    if cells === nothing
        cells = 1:getncells(grid)
    end
    # rest of constructor
end

It would then be possible to use the convenience methods in exactly the same way as before without running into headaches with the CellIterator

dh = DofHandler(grid, fe ∪ bd)
add!(dh, :u, Lagrange{RefQuadrilateral, 1}()^2)
close!()

for cell in CellIterator(dh)
    # do stuff
end

[KnutAM]

SubDofHandlers are intended to be used to define multiple element types and/or define different fields for seperate regions in the full domain.

In your case, this is exactly what you have, right? You have one field on part of your domain, and no fields in another part of your domain.

Ferrite provides convenience methods for this simplified usecase

Since what you describe is basically a variant of when the SubDofHandler is intended to be used, I don't think it makes sense to introduce a special syntax (and additional field in the DofHandler) for this special case.

But I do see the point that iterating over only the subdofhandlers sounds more logical, and is more useful. Another place where we assume the DofHandler to use the entire grid is during result export, e.g. when doing VTKFile(filename, dh) this is equivalent to VTKFile(filename, dh.grid), so we should have a consistent behavior there as well (which would be one export per subdofhandler in that case). This could actually solve the problem that cells not belonging to the subdofhandler gets colored, but is quite a lot more involved to solve...

[Joroks]

I agree! Changing the syntax for creation of the DofHandler was probabtly too ambitious. I would like to rework my PR such that it only changes the behavior of the CellIterator without making changes to either the DofHandler or the SubDofHandler.