def METHOD_NAME(Z, z_, k_n, g, v_D_mid, 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_)
    # Define midpoint values for brevity
    def avg(q, q_):
        return 0.5*(q + q_)
    sigma0_mid = avg(sigma0, sigma0_)
    sigma1_mid = avg(sigma1, sigma1_)
    v_mid = avg(v, v_)
    gamma_mid = avg(gamma, gamma_)
    # Define form
    n = FacetNormal(Z.mesh())
    F = (inner(inv(k_n)*A10(sigma0 - sigma0_), tau0)*dx
         + inner(A00(sigma0_mid), tau0)*dx
         + inner(inv(k_n)*A11(sigma1 - sigma1_), tau1)*dx
         + inner(div(tau0 + tau1), v_mid)*dx
         + inner(skw(tau0 + tau1), gamma_mid)*dx
         + inner(div(sigma0_mid + sigma1_mid), w)*dx
         + inner(skw(sigma0_mid + sigma1_mid), eta)*dx
         - inner(0.5*v_, (tau0 + tau1)*n)*ds(1)
         - inner(v_D_mid, (tau0 + tau1)*n)*ds(2) # Velocity on dO_D
         )
    # Tricky to enforce Dirichlet boundary conditions on varying sums
    # of components (same deal as for slip for Stokes for
    # instance). Use penalty instead
    beta = Constant(penalty_beta)
    h = tetrahedron.volume
    F_penalty = 0.5*(beta*inv(h)*inner((tau0 + tau1)*n,
                                       (sigma0 + sigma1)*n - g)*ds(1))
    F = F + F_penalty
    return F