Need tangent matrix with respect to velocities for deformable dofs

[xuchenhan-tri]

Background:

In the current FEM formulation for deformable bodies, we are solving for the equation

  Ma − f_d(x,v) − fₑ(x) = fₑₓₜ,   (1)

which then gets discretized as

 F(a) = M*a − f_d(x₀,v(a, v₀)) − fₑ(x(a, v₀, x₀)) - fₑₓₜ = 0,    (2)

with the Newmark scheme. We then use Newton's method to solve the discretized equations with a as the unknown variables. In that process, we produce the tangent matrix ∇ₐF = dF/da, which is related to, but not the same as, the tangent matrix ∇ᵥF = dF/dv consumed by the contact solver. It would be nice to be able to reuse the tangent matrix from FEM solve when Newton's method converges instead of recalculating another tangent matrix for contact.

Proposed solution:

Change variable in the discretization of (1) and solve

F(v) = dt*(M*a(v, a₀) − f_d(x₀, v) − fₑ(x(a₀, v, x₀)) - fₑₓₜ) = 0,  (3)

instead with v as the unknown variables. The tangent matrix needed for solving (3) with Newton's method is ∇ᵥF = dF/dv, so we can reuse the tangent matrix in the contact solve when the Newton's iterations converge. To enable this alternative formulation, we need a new class similar to NewmarkScheme that captures the same time discretization scheme, but uses v instead of a as the "unknown variable".

@amcastro-tri hinted at this idea in a PR review a while back.

[amcastro-tri]

Question. I see you got f_d(x₀, v). Shouldn't we generally have f_d(x(a₀, v, x₀), v)? that is, let the Scheme deal with the mapping from unknowns (in this case a) to FEM dofs (in this case x, v, a)

[xuchenhan-tri]

Technically yes, we should have f_d(x(a₀, v, x₀), v) instead. In practice, though, we force

x(a₀, v, x₀) = x₀ 

at the model level so we don't need to take the derivative of the stiffness matrix.

[xuchenhan-tri]

Closed via a combination of #16714 and #16800. The integrator provides the weights to combine stiffness, damping, and mass matrices, and FemModel simply calculates a linear combination.