BoundsError when creating OrderedFESpace on a periodic CartesianDiscreteModel

[bitsarebits]

Description

I encountered a BoundsError when attempting to instantiate an OrderedFESpace on a mesh with periodic boundary conditions.

The issue seems to be caused by the topology modification introduced by isperiodic=(true,) (where boundary vertices are topologically glued together but geometrically distinct). This likely breaks the continuous indexing assumptions made by the ordering algorithm in GridapROMs.DofMaps.OrderedFESpaces.

Minimal Working Example

Here is a self-contained script that reproduces the crash:

using Gridap
using GridapROMs

# 1D domain definition and partitioning
L = 1.0
Ω = (-L, L)
partition = (10,)

# Cartesian model creation with periodicity (the bug trigger)
model = CartesianDiscreteModel(Ω, partition, isperiodic=(true,))

# Reference element definition (1st-order Lagrangian)
reffe = ReferenceFE(lagrangian, Float64, 1)

# OrderedFESpace construction
println("Attempting to construct OrderedFESpace...")
V_ordered = OrderedFESpace(model, reffe)

Error / Stacktrace

Attempting to construct OrderedFESpace...
ERROR: LoadError: BoundsError: attempt to access 10-element LinearIndices{1, Tuple{Base.OneTo{Int64}}} at index [10:11]
Stacktrace:
  [1] throw_boundserror(A::LinearIndices{1, Tuple{Base.OneTo{Int64}}}, I::Tuple{UnitRange{Int64}})
    @ Base ./essentials.jl:15
  [2] checkbounds
    @ ./abstractarray.jl:699 [inlined]
  [3] view(A::LinearIndices{1, Tuple{Base.OneTo{Int64}}}, I::UnitRange{Int64})
    @ Base ./subarray.jl:214
  [4] get_dof_to_odof(fe_dof_basis::Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{1, Float64}, Int64}, cell_dofs_ids::Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, cells::SubArray{CartesianIndex{1}, 1, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, Tuple{Gridap.Arrays.IdentityVector{Int64}}, false}, onodes::LinearIndices{1, Tuple{Base.OneTo{Int64}}}, orders::Tuple{Int64})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:195
  [5] get_dof_to_odof
    @ ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:170 [inlined]
  [6] _get_cell_odof_info(model::CartesianDiscreteModel{1, Float64, typeof(identity)}, fe_dof_basis::FillArrays.Fill{Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{1, Float64}, Int64}, 1, Tuple{Base.OneTo{Int64}}}, cell_dofs_ids::Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, cell_to_parent_cell::Gridap.Arrays.IdentityVector{Int64}, orders::Tuple{Int64})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:146
  [7] get_cell_odof_ids(space::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.CompressedCellConformity{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:122
  [8] OrderedFESpace(f::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.CompressedCellConformity{FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:31
  [9] OrderedFESpace(model::CartesianDiscreteModel{1, Float64, typeof(identity)}, args::Tuple{Lagrangian, Tuple{DataType, Int64}, @Kwargs{}}; kwargs::@Kwargs{})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:43
 [10] OrderedFESpace(model::CartesianDiscreteModel{1, Float64, typeof(identity)}, args::Tuple{Lagrangian, Tuple{DataType, Int64}, @Kwargs{}})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/N03YZ/src/FEM/DofMaps/OrderedFESpaces.jl:42

Expected Behavior

OrderedFESpace should correctly map the degrees of freedom on a periodic topology without throwing bounds errors.

Environment

• Julia version: 1.12.6
• Gridap version: 0.20.8
• GridapROMs version: 0.1.4

[nichomueller]

Thank you for opening the issue! I have worked on the bug, and the tests should pass. I will close this issue

[bitsarebits]

Hi, thanks for the fast fix!

I switched to the main branch to test it, and the BoundsError during the OrderedFESpace instantiation is gone.

However, it seems to have exposed a downstream problem with the mapping. When triggering an interpolation on the space, it hits a new error: dof ids either positive or negative, not zero. Some periodic DoFs might have been left unmapped (ID 0) after the ordering fix.

Here is the updated MWE:

using Gridap
using GridapROMs
using GridapROMs.ParamDataStructures

# 1. Setup periodic model
domain = (0.0, 1.0)
partition = (10,)
model = CartesianDiscreteModel(domain, partition, isperiodic=(true,))
reffe = ReferenceFE(lagrangian, Float64, 1)

# 2. Construct OrderedFESpace (This now works correctly)
println("Step 1: Constructing OrderedFESpace...")
V = OrderedFESpace(model, reffe)

# 3. Setup Parametric Space
println("Step 2: Constructing TransientTrialParamFESpace...")
u₀(σ) = x -> 0.0
u₀ₚ(σ) = parameterise(u₀, σ)
U = TransientTrialParamFESpace(V, u₀ₚ)

# 4. Trigger the new error
println("Step 3: Triggering interpolation...")
σ_val = [[1.0]]
t0 = 0.0
uh = interpolate_everywhere(u₀ₚ(σ_val), U(σ_val, t0))

And the new stacktrace:

Step 1: Constructing OrderedFESpace...
Step 2: Constructing TransientTrialParamFESpace...
Step 3: Triggering interpolation...
ERROR: LoadError: dof ids either positive or negative, not zero
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] macro expansion
    @ ~/.julia/packages/Gridap/2PZdB/src/Helpers/Macros.jl:47 [inlined]
  [3] _free_and_dirichlet_values_fill!(free_vals::Vector{…}, dirichlet_vals::Vector{…}, cache_vals::Tuple{…}, cache_dofs::Tuple{…}, cell_vals::Gridap.Arrays.LazyArray{…}, cell_dofs::GridapROMs.DofMaps.OTable{…}, cells::UnitRange{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/UnconstrainedFESpaces.jl:135
  [4] gather_free_and_dirichlet_values!(fv::Vector{…}, dv::Vector{…}, f::OrderedFESpace{…}, cv::Gridap.Arrays.LazyArray{…})
    @ GridapROMs.DofMaps ~/.julia/packages/GridapROMs/WBAB0/src/FEM/DofMaps/OrderedFESpaces.jl:101
  [5] gather_free_and_dirichlet_values!(free_values::Vector{…}, dirichlet_values::Vector{…}, f::GridapROMs.ParamFESpaces.UnEvalTrialFESpace{…}, cell_vals::Gridap.Arrays.LazyArray{…})
    @ GridapROMs.ParamFESpaces ~/.julia/packages/GridapROMs/WBAB0/src/FEM/ParamFESpaces/UnEvalTrialFESpaces.jl:78
  [6] interpolate_everywhere!(object::Function, free_values::Vector{…}, dirichlet_values::Vector{…}, fs::GridapROMs.ParamFESpaces.UnEvalTrialFESpace{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/SingleFieldFESpaces.jl:217
  [7] interpolate_everywhere(object::Function, fs::GridapROMs.ParamFESpaces.UnEvalTrialFESpace{OrderedFESpace{Gridap.FESpaces.UnconstrainedFESpace{…}}, TrialFESpace{OrderedFESpace{…}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/2PZdB/src/FESpaces/SingleFieldFESpaces.jl:210