def METHOD_NAME(Z, z_, z__, k_n, g, v_D, ds):
    # Define trial and test functions
    (sigma0, sigma1, v, gamma) = TrialFunctions(Z)
    (tau0, tau1, w, eta) = TestFunctions(Z)
    # Extract previous components
    (sigma0_, sigma1_, v_, gamma_) = split(z_)
    (sigma0__, sigma1__, v__, gamma__) = split(z__)
    # Define complete form
    n = FacetNormal(Z.mesh())
    F = (inner(inv(k_n)*A10(1.5*sigma0 - 2.*sigma0_ + 0.5*sigma0__), tau0)*dx
         + inner(A00(sigma0), tau0)*dx
         + inner(inv(k_n)*A11(1.5*sigma1 - 2.*sigma1_ + 0.5*sigma1__), tau1)*dx
         + inner(div(tau0 + tau1), v)*dx
         + inner(skw(tau0 + tau1), gamma)*dx
         + inner(div(sigma0 + sigma1), w)*dx
         + inner(skw(sigma0 + sigma1), eta)*dx
         - inner(v_D, (tau0 + tau1)*n)*ds(2)
         )
    # Enforce essential bc on stress by penalty
    beta = Constant(penalty_beta)
    h = tetrahedron.volume
    F_penalty = beta*inv(h)*inner((tau0 + tau1)*n,
                                  (sigma0 + sigma1)*n - g)*ds(1)
    F = F + F_penalty
    return F