Solve for the coupled variables in the mixed function space and then conduct post - processing

General
1 
# 
P2 = element("Lagrange", domain.basix_cell(), 1, dtype=np.float64)  
P1 = element("Lagrange", domain.basix_cell(), 1, dtype=np.float64)  # 

# 
mixed_elem = mixed_element([P1, P1])  
W = fem.functionspace(domain, mixed_elem)

(phi1, phi2) = ufl.TrialFunctions(W)
(v1, v2) = ufl.TestFunctions(W)

a = (
    (1 / (3 * mu_a1) * inner(grad(phi1), grad(v1)) + mu_a * (phi1 - 2 / 3 * phi2) * v1) * dx
    + (1 / (7 * mu_a3) * inner(grad(phi2), grad(v2)) + mu_a * (4 / 9 * phi2 - 2 / 3 * phi1) * v2 + mu_a2 * (5 / 9 * phi2) * v2) * dx
    + (-(1 / 8) * phi1 + (7 / 24) * phi2) * v2 * ds + (1 / 2 * phi1 - 1 / 8 * phi2 ) * v1 * ds
)
L = (Q1 * v1 + Q2 * v2) * dx

problem = LinearProblem(a, L, bcs=[], 
                        petsc_options={
    "ksp_type": "gmres",
    "pc_type": "ilu",
    "ksp_rtol": 1e-8,
    "ksp_max_it": 1000,
}
                        )

solution = problem.solve()
**phi1_h, phi2_h = solution.split()**

How to perform post - processing of the linear addition of two variables in the mixed function space so that it can be visualized in ParaView? Specifically, how to interpolate phi1_h + phi2_h into a function.
thanks

2

See for instance dolfinx/python/demo/demo_elasticity.py at v0.9.0 · FEniCS/dolfinx · GitHub, where you would replace sigma_vm with phi1_h + phi2_h, and your W would be a P1 space.