ConstraintHandler for discontinuous interpolations

[AbdAlazezAhmed]

I haven't tested it with a problem, but it seems to work upon checking debugging mode.

[codecov-commenter]

Codecov Report

Patch coverage: 78.57% and project coverage change: -0.05 ⚠️

Comparison is base (185f8d6) 92.27% compared to head (bdbd430) 92.22%.

❗ Your organization is not using the GitHub App Integration. As a result you may experience degraded service beginning May 15th. Please install the Github App Integration for your organization. Read more.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #729      +/-   ##
==========================================
- Coverage   92.27%   92.22%   -0.05%     
==========================================
  Files          30       30              
  Lines        4594     4607      +13     
==========================================
+ Hits         4239     4249      +10     
- Misses        355      358       +3     

Impacted Files	Coverage Δ
src/interpolations.jl	96.53% <70.00%> (-0.72%)	⬇️
src/Dofs/ConstraintHandler.jl	95.98% <100.00%> (+<0.01%)	⬆️
src/FEValues/face_values.jl	98.13% <100.00%> (ø)

... and 4 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[AbdAlazezAhmed]

Setup:

julia> using Ferrite

julia> grid = generate_grid(Quadrilateral, (3, 3));

julia> ip = DiscontinuousLagrange{2, RefCube, 1}()
DiscontinuousLagrange{RefQuadrilateral, 1}()

julia> dh = DofHandler(grid);

julia> add!(dh, :u, 1,ip);
 (generic function with 10 methods)

julia> close!(dh);;

julia> ch = ConstraintHandler(dh);

julia> ∂Ω = ∪(
           getfaceset(grid, "left"),
           getfaceset(grid, "right"),
           getfaceset(grid, "top"),
           getfaceset(grid, "bottom"),

       );

julia> dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
Dirichlet(var"#5#6"(), Set(FaceIndex[FaceIndex((3, 2)), FaceIndex((7, 3)), FaceIndex((8, 3)), FaceIndex((1, 4)), FaceIndex((3, 1)), FaceIndex((9, 3)), FaceIndex((2, 1)), FaceIndex((4, 4)), FaceIndex((7, 4)), FaceIndex((6, 2)), FaceIndex((9, 2)), FaceIndex((1, 1))]), :u, Int64[], Int64[], Int64[])

julia> add!(ch, dbc);

julia> close!(ch)
ConstraintHandler:
  BCs:
    Field: u, Components: [1]

Master branch:

julia> ch.prescribed_dofs
Int64[]

This PR:

julia> ch.prescribed_dofs
20-element Vector{Int64}:
  1
  2
  4
  5
  ⋮
 34
 35
 36

[fredrikekre]

I wonder if we can define facedof_indices for the interpolations, instead? That is only used for boundary conditions I think.

What happens for the piecewise constant (order 0) interpolations? I suppose for those one need to impose the boundary condition by integrating the jump, just as is done for element-element couplings?

[AbdAlazezAhmed]

Indeed, the get_continuous_interpolation() approach is faulty for piecewise constant interpolations. I'll check it, thanks!

[termi-official]

I wonder if we can define facedof_indices for the interpolations, instead? That is only used for boundary conditions I think.

Note that this one is also designed for the distributed assembly stuff to sync process boundaries - I just started chery picking stuff from the PR to make the review easier.

What happens for the piecewise constant (order 0) interpolations? I suppose for those one need to impose the boundary condition by integrating the jump, just as is done for element-element couplings?

If we have no basis function on the interior of the boundary of the element, then AFAIK we cannot use the strong Dirichlet enforcement anymore and fall back to penalty approaches.

[fredrikekre]

If we have no basis function on the interior of the boundary of the element, then AFAIK we cannot use the strong Dirichlet enforcement anymore and fall back to penalty approaches.

Right, but penalty approach (integrating the jump from the element to the prescribed value(?)) would work also for higher order elements, right? This PR is for supporting strong Dirichlet when possible then, which should simply not be supported for order 0.

[AbdAlazezAhmed]

This PR is for supporting strong Dirichlet when possible then, which should simply not be supported for order 0.

Indeed. It currently prescribes no dofs for order 0.

Note that this one is also designed for the distributed assembly stuff to sync process boundaries

Maybe we can do something like:

function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=facedof_indices) where F
    @assert issorted(components)
    local_face_dofs = Int[]
    local_face_dofs_offset = Int[1]
    getorder(interpolation) == 0 && return [], [1 for _ in 1:nfaces(interpolation)+1]
    for (_, face) in enumerate(boundaryfunc(IsDiscontinuous(interpolation) ? get_continuous_interpolation(interpolation) : interpolation))
        for fdof in face, d in 1:field_dim
            if d in components
                push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
            end
        end
        push!(local_face_dofs_offset, length(local_face_dofs) + 1)
    end
    return local_face_dofs, local_face_dofs_offset
end

[termi-official]

This PR is for supporting strong Dirichlet when possible then, which should simply not be supported for order 0.

Indeed. It currently prescribes no dofs for order 0.

Yea this case is quite tricky. Users may accidentally prescribe different values on different faces in the corners of a domain which leads to weird ambiguities. It also raises the question at which point of the boundary should a constraint be evaluated. We should just throw a warning in this case and do nothing - which should be reflected in the tests, because this feels like there are possibilities to introduce unexpected regressions.

Note that this one is also designed for the distributed assembly stuff to sync process boundaries

Maybe we can do something like:

function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=facedof_indices) where F
    @assert issorted(components)
    local_face_dofs = Int[]
    local_face_dofs_offset = Int[1]
    getorder(interpolation) == 0 && return [], [1 for _ in 1:nfaces(interpolation)+1]
    for (_, face) in enumerate(boundaryfunc(IsDiscontinuous(interpolation) ? get_continuous_interpolation(interpolation) : interpolation))
        for fdof in face, d in 1:field_dim
            if d in components
                push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
            end
        end
        push!(local_face_dofs_offset, length(local_face_dofs) + 1)
    end
    return local_face_dofs, local_face_dofs_offset
end

I should also note that the get_continuous_interpolation trick will only work for the implemented discontinuous interpolation, which is not the usually used one. Most of the interpolations in DG do not have the nodes on the boundary, but for example at the Gauss-Legendre points (or other ones, depending on your needs). It is still nice (and wanted) to have strong Dirichlet enforcement for the existing interpolation! This e.g. allows us to compare approaches in a unified framework.

[AbdAlazezAhmed]

I was thinking maybe we could add something like dirichlet_*dof_indices and make it default to *dof_indices . Or maybe make use of IsDiscontinuous when dealing with distributed assembly instead and just dispatch *dof_indices. What do you think?

[fredrikekre]

dirichlet_*dof_indices sounds like a good plan (somewhat similar to

Ferrite.jl/src/interpolations.jl

Line 53 in 0ec0e49

n_dbc_components(ip::Interpolation) = n_components(ip)

)

[fredrikekre]

Looks good, perhaps add a simple test?

[AbdAlazezAhmed]

Looks good, perhaps add a simple test?

Sure, thanks!

[termi-official]

Can you add tests for the periodic case as regression tests, too? Also, I think the devdocs are missing.

[AbdAlazezAhmed]

Can you add tests for the periodic case as regression tests, too? Also, I think the devdocs are missing.

I'm not sure what regression tests can be written for the periodic case as I haven't changed its code, it throws the same mirroring error both before and after this PR.