CellIterator causes a OutOfMemoryError

[Joroks]

This only happens if the cellset used for the FieldHandler doesn't contain every cell in the grid.

using Ferrite

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

FE = MixedDofHandler(grid)
fe_domain = Set(1:4) # with = Set(1:27) the error doesn't happen
fe_field = Field(:u, Lagrange{3, RefCube, 1}(), 3)
fe_fieldHandler = FieldHandler([fe_field], fe_domain)

add!(FE, fe_fieldHandler)
close!(FE)

begin
    cell = CellCache(FE)

    for i in fe_domain
        reinit!(cell, i)
    end

    println("works")
end

begin
    for cell in CellIterator(FE) # causes OutOfMemoryError

    end
end

I'm on Ferrite version 0.3.14

here is my versioninfo()

Julia Version 1.10.3
Commit 0b4590a550 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × AMD Ryzen 7 3700X 8-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
Threads: 8 default, 0 interactive, 4 GC (on 16 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 8

and here is the Stacktrace

LoadError: OutOfMemoryError()
Stacktrace:
  [1] _growend!
    @ .\array.jl:1072 [inlined]
  [2] resize!
    @ .\array.jl:1315 [inlined]
  [3] reinit!(cc::CellCache{Vec{3, Float64}, Grid{3, Hexahedron, Float64}, MixedDofHandler{3, Float64, Grid{3, Hexahedron, Float64}}}, i::Int64)
    @ Ferrite C:\Users\johan\.julia\packages\Ferrite\5Ozov\src\iterators.jl:96
  [4] iterate(ci::CellIterator{CellCache{Vec{3, Float64}, Grid{3, Hexahedron, Float64}, MixedDofHandler{3, Float64, Grid{3, Hexahedron, Float64}}}, UnitRange{Int64}}, state_in::Int64)
    @ Ferrite C:\Users\johan\.julia\packages\Ferrite\5Ozov\src\iterators.jl:176
  [5] top-level scope
    @ c:\Users\johan\dev\capriccio-julia\demos\OutOfMemory.jl:26
  [6] eval
    @ .\boot.jl:385 [inlined]
  [7] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)        
    @ Base .\loading.jl:2076
  [8] include_string(m::Module, txt::String, fname::String)
    @ Base .\loading.jl:2086
  [9] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base .\essentials.jl:892
 [10] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base .\essentials.jl:889
 [11] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:271
 [12] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:181
 [13] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [14] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [15] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [16] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:150
 [17] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:515
 [18] with_logger
    @ .\logging.jl:627 [inlined]
 [19] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [20] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [21] invokelatest(::Any)
    @ Base .\essentials.jl:889
 [22] (::VSCodeServer.var"#64#65")()
    @ VSCodeServer c:\Users\johan\.vscode\extensions\julialang.language-julia-1.79.2\scripts\packages\VSCodeServer\src\eval.jl:34

[termi-official]

Hi @Joroks . Thanks for the report and welcome to the community!

A possible workaround is to pass the subdomain manually into the iterator

using Ferrite

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

FE = MixedDofHandler(grid)
fe_domain = Set(1:4) # with = Set(1:27) the error doesn't happen
fe_field = Field(:u, Lagrange{3, RefCube, 1}(), 3)
fe_fieldHandler = FieldHandler([fe_field], fe_domain)

add!(FE, fe_fieldHandler)
close!(FE)

for cell in CellIterator(FE, Set(1:4) )
end

Interestingly the issue still persists in some form on master, although it manifests differently:

using Ferrite

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

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

for cell in CellIterator(sdh) # works
end

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

for cell in CellIterator(dh) # errors
end

error message

BoundsError: attempt to access 1-element Vector{SubDofHandler{DofHandler{3, Grid{3, Hexahedron, Float64}}}} at index [0]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] ndofs_per_cell
   @ ~/Repos/Ferrite.jl/src/Dofs/DofHandler.jl:186 [inlined]
 [3] reinit!(cc::CellCache{Vec{3, Float64}, Grid{3, Hexahedron, Float64}, DofHandler{3, Grid{3, Hexahedron, Float64}}}, i::Int64)
   @ Ferrite ~/Repos/Ferrite.jl/src/iterators.jl:85
 [4] iterate(iterator::CellIterator{CellCache{Vec{…}, Grid{…}, DofHandler{…}}, UnitRange{Int64}}, state_in::Int64)
   @ Ferrite ~/Repos/Ferrite.jl/src/iterators.jl:392
 [5] top-level scope
   @ ./REPL[10]:

This can be tracked down to the either
A) Celliterator constructor always using the full domain.
B) A missing check in ndofs_per_cell(dh::DofHandler, cell::Int) not checking if the subdofhandler exists.

I will fix the later first.