Constraining dof > ndofs(dh) fails

[KnutAM]

With the nice new SparsityPattern (#888), we can now extend the equation system beyond ndofs(dh). (I'm seeing how #596 could work now without changing the dofhandler).

I realized now that the ConstraintHandler fails if adding constraints to these, since it assumes the maximum size is ndofs(dh). (This now happens in _sorted_setdiff! from #860 , before I think ch.free_dofs would simply be wrong but not error). In other places, we use ndofs(ch.dh) to check dimensions of the "to-be-constrained" vector, but this should really check for last(ch.prescribed_dofs) to be safe from out of bounds access.

An easy solution would be to add an optional keyword arg to ConstraintHandler(dh::DofHandler; n_dofs=ndofs(dh)), but I'm also not sure what we should call this since ndofs is taken for the function name. Anyone have some thoughts on this?

MWE

grid = generate_grid(Quadrilateral, (2,2));
dh = close!(add!(DofHandler(grid), :u, Lagrange{RefQuadrilateral,1}()));
ch = ConstraintHandler(dh);
Ferrite.add_prescribed_dof!(ch, ndofs(dh) + 1, NaN, nothing)
close!(ch)

Gives ERROR: BoundsError: attempt to access 8-element Vector{Int64} at index [9]

[KnutAM]

Posting workaround here that can be used untill we decide how to fix this

struct DHWrapper{DH<:DofHandler} <: Ferrite.AbstractDofHandler
    dh::DH
    ndofs::Int 
end
@inline Base.getproperty(dh::DHWrapper, f::Symbol) = Base.getproperty(Base.getfield(dh, :dh), f)
@inline Ferrite.ndofs(dh::DHWrapper) = Base.getfield(dh, :ndofs)

which allows

ch = ConstraintHandler(DHWrapper(dh, ndofs(dh) + num_extra))
...

and seems to work fine (which is very nice to see regarding hackability).

[lijas]

If a user wants to add global Lagrange multipliers/dofs, I think it should be possible to add these to the DofHandler directly, e.g. add_global_dofs!(dh, :lambda, 3). The sparsity-pattern can then be handled by the user. It seems reasonable that DofHandler should be aware of all the dofs of the problem.

[KnutAM]

Yes, that would be the solution here I guess: https://github.com/Ferrite-FEM/Ferrite.jl/pull/596/files (but with better names, global dofs is probably a good one) It is nice to be able to store the dofs associated with a "field"

An alternative is adding a "GlobalInterpolation" that we discussed earlier on slack. These dofs are part of the celldof for each cell in the subdofhandler (https://github.com/Ferrite-FEM/Ferrite.jl/compare/master...kam/GlobalDofs). It would also be possible to have two variants, one connected to cells and one not connected?

[termi-official]

If a user wants to add global Lagrange multipliers/dofs, I think it should be possible to add these to the DofHandler directly, e.g. add_global_dofs!(dh, :lambda, 3). The sparsity-pattern can then be handled by the user. It seems reasonable that DofHandler should be aware of all the dofs of the problem.

What is wrong about adding a

struct Point <: AbstractCell nodes::Tuple{Int} end

(and e.g. a Connector which behaves similar to a line element)?

I kinda started liking this idea, because it is also nice for the integration with visualization pipelines.

[KnutAM]

What is wrong about adding a Point <: AbstractCell

To me, the issue is that the coordinate is a dummy coordinate. In my application, the dofs could represent e.g. strain components, and I don't like if we require adding coordinates to this.

For a flexible Connector, a line element is nice and could be used. But if it is rigid, I guess it leads to numerical problems unless introduced as an AffineConstraint, and then it wouldn't really be a line element any more?

[termi-official]

What is wrong about adding a Point <: AbstractCell

To me, the issue is that the coordinate is a dummy coordinate. In my application, the dofs could represent e.g. strain components, and I don't like if we require adding coordinates to this.

Why? Isn't the strain associated to some point or region?

For a flexible Connector, a line element is nice and could be used. But if it is rigid, I guess it leads to numerical problems unless introduced as an AffineConstraint, and then it wouldn't really be a line element any more?

I was not talking about actual beams or springs here, but about some entity that connects with other entities. For example you have one additional dof, e.g. some pressure within a sphere, connected to the elements on the sphere inner surface. So you can have fine grained control on your geometry while setting up your problem and the dof handler handles everything for you (assuming you also provide some dummy interpolation to tell the dof handler how the dofs need to be connected). If you want to introduce some affine constraint with it, then you can now easily query all information from your new element instead of doing it manually.

[KnutAM]

Why? Isn't the strain associated to some point or region?

The entire domain. But the point is, the coordinate has not meaning.

I was not talking about actual beams or springs here, but about some entity that connects with other entities.

Aha, now I understand. Yes, for such cases I agree that adding an additional physical node can make sense (especially for visualization purposes). When doing that for my last conference, I needed to integrate over the face though, so I couldn't have used CellValues if I had interpreted that as cells.

I think you proposed this earlier, but to this problem I think I would find it most natural to have a custom SubDofHandler to which one could add face-regions that are connected to x number of "global" dofs. Something like

dh = DofHandler(grid)
main_sdh = SubDofHandler(dh, 1:getncells(grid))
add!(main_sdh, :p, ipp) # Regular fields
add!(main_sdh, :u, ipu) # Regular fields
for pore in pores # collection of OrderedSet{FacetIndex} or similar
    pore_fdh = FacetDofHandler(dh, pore)
    add!(pore_fdh, :pp, 1) # One pore pressure
    add!(pore_fdh, :whatever, <num_pore_dofs>)
end
close!(dh)

Since for each pore_fdh one would need to loop over all its faces?

But my conclusion is that supporting this type of extra dof is different from supporting Lagrange-multiplier dofs or AffineConstraint-only dofs.
I think it would be better to provide such different ways of doing these tasks, rather than trying to solve both cases with the same solution.

Hence, my proposed solution would be to go with @lijas suggestion to be able to add extra global dofs to the dofhandler - i.e. saying that ndofs(dh) should always return the true number of dofs. I'll try to make a PR for this with minimal changes.

I also support @termi-official suggestion that we should provide a way to solve these surface-interactions, and I think this is an interesting discussion to figure out what is the most flexible way to do!

[termi-official]

Thanks for elaborating.

I was not talking about actual beams or springs here, but about some entity that connects with other entities.

Aha, now I understand. Yes, for such cases I agree that adding an additional physical node can make sense (especially for visualization purposes). When doing that for my last conference, I needed to integrate over the face though, so I couldn't have used CellValues if I had interpreted that as cells.

I think you proposed this earlier, but to this problem I think I would find it most natural to have a custom SubDofHandler to which one could add face-regions that are connected to x number of "global" dofs. Something like

dh = DofHandler(grid)
main_sdh = SubDofHandler(dh, 1:getncells(grid))
add!(main_sdh, :p, ipp) # Regular fields
add!(main_sdh, :u, ipu) # Regular fields
for pore in pores # collection of OrderedSet{FacetIndex} or similar
    pore_fdh = FacetDofHandler(dh, pore)
    add!(pore_fdh, :pp, 1) # One pore pressure
    add!(pore_fdh, :whatever, <num_pore_dofs>)
end
close!(dh)

[...]

I also support @termi-official suggestion that we should provide a way to solve these surface-interactions, and I think this is an interesting discussion to figure out what is the most flexible way to do!

No need for this. For most practical purposes you can add some elements to your mesh which model the surface and simply add a SubDofHandler on this subdomain. This should already work out of the box on master. I will try this out in Thunderbolt soon and report back when I run into issues and we need intervention. Currently I simply manage the additional dofs manually.

[KnutAM]

No need for this. For most practical purposes you can add some elements to your mesh which model the surface and simply add a SubDofHandler on this subdomain. I will try this out in Thunderbolt soon and report back when I run into issues and we need intervention.

Cool - looking forward to seeing the solution! (My solution with FerriteAssembly is to hacky at the moment since it was pre-SparsityPattern, and I don't have time to work on that example again in the coming couple of months, so after seeing your solution I'm hoping it will be a quick fix:))