Is it possible to apply periodic contraints directly to the DofHandler?

[Joroks]

I'm currently running simulations that couples FEM with molecular dynamics (MD). Here, I need to account for the periodic boundary constraints often used with MD on the FEM side as well. The simulation is purely dynamic and the coupling is enforced through Lagrange multipliers, which results in a system of equations that look something like this:

Therefore it would be necessary to account for the periodicity not only in M_FE, ü_FE and f_FE, but also in G_FE. The simplest solution I can come up with, is to simply assing the same dofs to the image and mirror nodes within the DofHandler (For this application the periodic boundarys are allways parallel to each other). With this I could avoid explicitly accounting for the periodicity in the assembling routines for M_FE, G_FE, and f_fe. Is this something that can be done with Ferrite? I also don't mind, if this requires me to implement this my self by meddleing with Ferrite internals.

P.s. The compoments of the Matrix on the left-hand side are constant throughout the entire simulation and therefore only get assembled once. I also never actually assemble the Matrix on the left-hand side, as I have manipulated the individual components of this equation in a manner that is more efficient to solve.

[KnutAM]

I think in principle yes (you would have to renumber the celldofs manually in the dof-handler), and update the ndofs field, and then I think it would work. However, the more standard way would be to use the AffineConstraints created when using the PeriodicDirichlet with the ConstraintHandler.

[Joroks]

Thanks! I took a closer look at how periodic and affine constraints work in Ferrite and came up with this solution. It's certainly not the most elegant or efficient way to do this, but it appears to do the job.

function apply_periodicity!(dh::DofHandler, face_mapping)
    ch = ConstraintHandler(dh)
    add!(ch, PeriodicDirichlet(:u, face_mapping))
    close!(ch)
    C, _ = Ferrite.create_constraint_matrix(ch)
    perm = permutedims(C).rowval
    Ferrite._renumber!(dh, perm)
    dh.ndofs = size(C, 2)
    return dh
end

[fredrikekre]

Why don't you want to use the ConstraintHandler directly?

[Joroks]

I think it comes down to the fact, that modifying the DofHandler directly instead of using the ConstraintHandler, allows me to account for the periodic boundary contitions only once when setting up my DofHandler. This then just makes the periodic boundary constraints work without any further adjustments to my simulation code. Using the ConstraintHandler, I have to account for this constraint at multiple places in my simulation code. I honestly struggled a bit, trying to figure out, at what point of my simulation I have to call apply! and on what. I also struggled reasoning about, how this would interact with my coupling setup.