Periodic boundary conditions cant be applied because assembly changes sparsity pattern?

[alex-fu27]

What I want to achieve is a transient simulation on a 1D grid with periodicity in the left/right boundaries (what comes out on the right goes in on the left and vice versa).

So I assemble two matrices, M and K, which I combine to A = M .+ (Δt .* K).
Then, I build a ConstraintHandler, following the "periodic boundary conditions" documentation:

ch = ConstraintHandler(dh)
mapping = collect_periodic_facets(grid, "left", "right", (x) -> x - Ferrite.Vec{2}((2.0, 0.0)))
pdbc = PeriodicDirichlet(:c, mapping, [1,])
add!(ch, pdbc)
close!(ch)

This ch displays as an empty ConstraintHandler

ConstraintHandler:
  BCs:

so I initially suspected I was creating my constraint handler wrong?

Ferrite finally errors out when I try to apply this constraint handler to my matrix:

rhs_data = get_rhs_data(ch, A)
apply!(A, ch)

SparsityError: writing to an index outside the sparsity pattern is not allowed

Stacktrace:
 [1] addindex!(A::SparseArrays.SparseMatrixCSC{Float64, Int64}, v::Float64, i::Int64, j::Int64)
   @ Ferrite ~/.julia/packages/Ferrite/Iyg5J/src/arrayutils.jl:76
 [2] _condense!(K::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::Vector{Float64}, dofcoefficients::Vector{Union{Nothing, Vector{Pair{Int64, Float64}}}}, dofmapping::Dict{Int64, Int64}, sym::Bool)
   @ Ferrite ~/.julia/packages/Ferrite/Iyg5J/src/Dofs/ConstraintHandler.jl:707
 [3] apply!(KK::SparseArrays.SparseMatrixCSC{Float64, Int64}, f::Vector{Float64}, ch::ConstraintHandler{DofHandler{2, Grid{2, Quadrilateral, Float64}}, Float64}, applyzero::Bool)
   @ Ferrite ~/.julia/packages/Ferrite/Iyg5J/src/Dofs/ConstraintHandler.jl:641
 [4] apply!(K::SparseArrays.SparseMatrixCSC{Float64, Int64}, ch::ConstraintHandler{DofHandler{2, Grid{2, Quadrilateral, Float64}}, Float64})
   @ Ferrite ~/.julia/packages/Ferrite/Iyg5J/src/Dofs/ConstraintHandler.jl:591
 [5] top-level scope
   @ ~/jl_notebook_cell_df34fa98e69747e1a8f8a730347b8e2f_X10sZmlsZQ==.jl:2

But if I run apply! on a newly created matrix, it works:

B = allocate_matrix(dh)
_ = get_rhs_data(ch, B)
apply!(B, ch)

What am I doing wrong? Any help is much appreciated, so I've a made a small non-working example 😄 for your reference.

[fredrikekre]

1. You need to pass the constraint handler to allocate_matrix becuase the constraints changes the sparsity pattern in the case of periodic constraints, e.g.

   M = allocate_matrix(dh, ch)
   K = allocate_matrix(dh, ch)

2. The operation A = M .+ (Δt .* K) results in A having different sparsity pattern because some values in M and K are zero, and these indices are not preserved in broadcasting... I don't know what that best way to fix this is, but perhaps build A like this instead (which works since we know M and K have identical patterns):

   A = copy(M)
   A.nzval .= M.nzval .+ Δt .* K.nzval

[alex-fu27]

Thanks for the tip!
I guess it's working now, I didn't expect the oscillating reflections from the border, but I'm new to FEM in general, so I can't tell if that is 'normal'.

Regarding the second problem, I've patched addindex! to allow adding to indices not in the sparsity pattern. I've created a pull request: #1169

[koehlerson]

Could be that you are plotting wrong by assuming dof order == node order

[koehlerson]

or maybe you forgot to apply the periodic BC also to your solution

[alex-fu27]

Oh no, yes! I forgot to apply the boundary conditions to the solutions. Now it's working. Thank you a lot 😄

[lijas]

so I initially suspected I was creating my constraint handler wrong?

The PeriodicDirichlet object is not stored in the ConstrainHandler, but constraint are still added.

I dont think anyone has tried combining PBC with RHSData, so I am not sure if it works. But try Fredriks solution and let us know if that solves your problem.