update hyperelasticity docs to show how to use NonlinearSolve

[oscardssmith]

I've also tweaked the manual version a little to use a slightly more standard newton iteration.

[termi-official]

Just to be sure that I totally understand this correctly, in NonlinearSolve.jl there is currently no support to evaluate both at the same time, right? If so, then we can separate the functions above into jacobian and residual only:

function assemble_element_jac!(ke, cell, cv, fv, mp, ue, ΓN)
    ## Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    b = Vec{3}((0.0, -0.5, 0.0)) # Body force
    tn = 0.1 # Traction (to be scaled with surface normal)
    ndofs = getnbasefunctions(cv)

    for qp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, qp)
        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇u = function_gradient(cv, qp, ue)
        F = one(∇u) + ∇u
        C = tdot(F) # F' ⋅ F
        ## Compute stress and tangent
        S, ∂S∂C = constitutive_driver(C, mp)
        P = F ⋅ S
        I = one(S)
        ∂P∂F =  otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I)

        ## Loop over test functions
        for i in 1:ndofs
            ## Test function and gradient
            δui = shape_value(cv, qp, i)
            ∇δui = shape_gradient(cv, qp, i)
            ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation
            for j in 1:ndofs
                ∇δuj = shape_gradient(cv, qp, j)
                ## Add contribution to the tangent
                ke[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ
            end
        end
    end
end;

function assemble_element_residual!(ge, cell, cv, fv, mp, ue, ΓN)
    ## Reinitialize cell values, and reset output arrays
    reinit!(cv, cell)
    fill!(ke, 0.0)
    fill!(ge, 0.0)

    b = Vec{3}((0.0, -0.5, 0.0)) # Body force
    tn = 0.1 # Traction (to be scaled with surface normal)
    ndofs = getnbasefunctions(cv)

    for qp in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, qp)
        ## Compute deformation gradient F and right Cauchy-Green tensor C
        ∇u = function_gradient(cv, qp, ue)
        F = one(∇u) + ∇u
        C = tdot(F) # F' ⋅ F
        ## Compute stress and tangent
        S, ∂S∂C = constitutive_driver(C, mp)
        P = F ⋅ S

        ## Loop over test functions
        for i in 1:ndofs
            ## Test function and gradient
            δui = shape_value(cv, qp, i)
            ∇δui = shape_gradient(cv, qp, i)
            ## Add contribution to the residual from this test function
            ge[i] += ( ∇δui ⊡ P - δui ⋅ b ) * dΩ
        end
    end

    ## Surface integral for the traction
    for face in 1:nfaces(cell)
        if (cellid(cell), face) in ΓN
            reinit!(fv, cell, face)
            for q_point in 1:getnquadpoints(fv)
                t = tn * getnormal(fv, q_point)
                dΓ = getdetJdV(fv, q_point)
                for i in 1:ndofs
                    δui = shape_value(fv, q_point, i)
                    ge[i] -= (δui ⋅ t) * dΓ
                end
            end
        end
    end
end;

+- typos and we can do the same with with assemble_global!

I will benchmark this after I finish the prioritized tasks.

[KristofferC]

Wouldn't it be better to keep them in the same function with a flag to chose what to evaluate. That would make it easier to evaluate both for solvers that support it, to share some computations?

[termi-official]

I am not sure about the performance hit associated with splitting up the Jacobian and residual computation, especially for element formulations which come with some kind of condensation for the Jacobian.

[KnutAM]

I think this is really nice to include - thanks!

Would it make sense to put this in a how-to? E.g. "Using NonlinearSolve.jl" to (a) increase the visibility and (b) simplify the present tutorial. In that case, I would still add a link to that how-to in the hyperelasticity tutorial to point users to that option.

[KristofferC]

What's the reason for passing ΓN etc through p instead of just directly using them from the scope above (since this is a closure closing over it)?

[termi-official]

Would it make sense to put this in a how-to? E.g. "Using NonlinearSolve.jl" to (a) increase the visibility and (b) simplify the present tutorial. In that case, I would still add a link to that how-to in the hyperelasticity tutorial to point users to that option.

I am also struggling with this one. It is, like the Navier-Stokes example, somehow both. I mean on the one side it is learning-oriented and shows learning how which makes it a tutorial (https://docs.divio.com/documentation-system/tutorials/). But you can also frame the same example as a How do I ...?, which makes it also a how-to (https://docs.divio.com/documentation-system/how-to-guides/). Any thoughts on this @fredrikekre ?

We also should put the results of our discussion into the CONTRIBUTING.md @KnutAM because I think the distinction between how-to and tutorial is quite narrow for our framework and it can help to guide new contributors.

[KnutAM]

I am also struggling with this one. It is, like the Navier-Stokes example, somehow both.

For cases when the main purpose is to show how to interact with different packages it could be nice to have a separate category for that when we have enough examples.

In this case, I find this an "add-on" to the current tutorial. So instead of showing both ways, I think new users will find it easier if only one method is shown in the tutorial, and users are referred to how-to's for different extensions.

To me, the difference in Ferrite between how-to and tutorial is that a tutorial should include all theory and (typically) a complete solution.

My opinion is that how-to's are nice for building on top of existing examples (such as the multithreading how-to, except that the linear elasticity tutorial isn't merged yet...).

We also should put the results of our discussion into the CONTRIBUTING.md @KnutAM because I think the distinction between how-to and tutorial is quite narrow for our framework and it can help to guide new contributors.

👍