AffineConstraint for nonzero b and nonlinear problem + docfixes

[KnutAM]

This PR fixes an issue recently discussed on Slack. The test code and derivations are included here for reference.
In addition, it adds further docstrings relevant when using AffineConstraints.

Test code

Without this PR, the nonlinear problem doesn't converge. With this PR, the code shows how the wrong results are obtained if using the `a=-K\r` instead of `a=K\f` for the linear solution method, which this PR documents.

using Ferrite

params = (k=1.0, f=1.0, a=1.0, b=0.2, tol=1e-10, maxiter=2)

grid = generate_grid(Line, (2,))
addfaceset!(grid, "center", x->x[1]≈0.0)

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh)

eval_norm(r, ch) = sqrt(sum(i->r[i]^2, Ferrite.free_dofs(ch)))

function doassemble!(K, r, dh, a, params)
    # Spring elements 
    k = params.k
    Ke = [k -k; -k k]
    # Quick and dirty assem
    assembler = start_assemble(K, r)
    for cell in CellIterator(dh)
        ae = a[celldofs(cell)]
        re = Ke*ae
        assemble!(assembler, celldofs(cell), Ke, re)
    end
    r[3] -= params.f # Force on the main dof (right side)
end

ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfaceset(grid, "center"), (x,t)->Vec{1}((0.0,))))
add!(ch, AffineConstraint(1, [3=>params.a], params.b))
close!(ch)

K = create_sparsity_pattern(dh, ch)
r = zeros(ndofs(dh))

a = zeros(ndofs(dh))
update!(ch, 0.0)

# Nonlinear solution
apply!(a, ch) 
for niter = 0:params.maxiter
    doassemble!(K, r, dh, a, params)
    apply_zero!(K, r, ch)
    rnorm = eval_norm(r, ch)
    println("$niter: $rnorm, $r, $a")
    rnorm < params.tol && break
    Δa = -K\r 
    apply_zero!(Δa, ch)
    a .+= Δa
end
eval_norm(r, ch) > params.tol && @warn "Did not converge"
a_nonlinear = copy(a)

# Linear solution (r)
a_linear_r = zero(a)
doassemble!(K, r, dh, a_linear_r, params)
apply!(K, r, ch)
a_linear_r = -K\r 
apply!(a_linear_r, ch)

# Linear solution (f)
a_linear_f = zero(a)
doassemble!(K, r, dh, a_linear_f, params)
f = -r
apply!(K, f, ch)
a_linear_f = K\f
apply!(a_linear_f, ch)

# Analytical comparison
a3 = (params.f/params.k - params.b)/(1 + params.a)
a_analytical = [params.a*a3+params.b 0.0 a3]
@show a_linear_r
@show a_linear_f
@show a_nonlinear
@show a_analytical

Analytical solution

Spring

|------ $c$ ------|

|-- $k$ --|-- $k$ --|--> $f$

$u_1$ -> $u_2$ -> $u_3$ ->

Kinematical constraints

• Dirichlet: $u_2 = 0$
• Affine: $u_1 = au_3 + b$

Free body diagrams

<- $ku_1$ -- (1) -- $f_c$ ->

<- $(f_c+ku_3)$ -- (3) -- $f$ ->

where $f_c$ is the affine constraint force and $f$ is the external force

Equilibrium

• $0 = f_1 = -k u_1 + f_c$
• $0 = f_3 = -k u_3 - f_c + f$
• $0 = f_1+f_3 = -k (u_1 + u_3)$
• $f = k u_1 + k u_3$
• $f = k (u_3 (a+1) + b)$
• $u_3 = (f/k - b)/(1+a)$

[codecov-commenter]

Codecov Report

Base: 92.76% // Head: 92.76% // No change to project coverage 👍

Coverage data is based on head (2aa0137) compared to base (ca6bd6f).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #575   +/-   ##
=======================================
  Coverage   92.76%   92.76%           
=======================================
  Files          28       28           
  Lines        4202     4202           
=======================================
  Hits         3898     3898           
  Misses        304      304           

Impacted Files	Coverage Δ
src/Dofs/ConstraintHandler.jl	96.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

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

[fredrikekre]

LGMT, just some nits.

[KnutAM]

Windows test failure is unrelated to this PR, it also occurs on master.