Simplify `create_formulation` for mixed domain approach

@jorgensd suggested we simplify the create_formulation method with:

```python
entity_maps = {
sd.submesh: sd.parent_to_submesh for sd in self.volume_subdomains
}

list_of_spaces = [sd.u.function_space for sd in self.volume_subdomains]
W = ufl.MixedFunctionSpace(*list_of_spaces)
v = ufl.TestFunctions(W)

all_mobile_species = [spe for spe in self.species if spe.mobile]
if len(all_mobile_species) > 1:
raise NotImplementedError("Multiple mobile species not implemented")
H = all_mobile_species[0]
us = [H.subdomain_to_solution[sd] for sd in self.volume_subdomains]

for interface in self.interfaces:
gamma = interface.penalty_term

subdomain_0, subdomain_1 = interface.subdomains
i0 = self.volume_subdomain.index(subdomain_0)
i1 = self.volume_subdomain.index(subdomain_1)
res = interface.restriction
n_0 = n(res[0])
h_0 = 2 * cr(res[0])
h_1 = 2 * cr(res[1])

v_b = v[i0](res[0])
v_t = v[i1](res[1])

u_b = H.subdomain_to_solution[subdomain_0](res[0])
u_t = H.subdomain_to_solution[subdomain_1](res[1])

K_b = subdomain_0.material.get_solubility_coefficient(
self.mesh.mesh, self.temperature_fenics(res[0]), H
)
K_t = subdomain_1.material.get_solubility_coefficient(
self.mesh.mesh, self.temperature_fenics(res[1]), H
)

F_0 = -0.5 * mixed_term((u_b + u_t), v_b, n_0) * dInterface(
interface.id
) - 0.5 * mixed_term(v_b, (u_b / K_b - u_t / K_t), n_0) * dInterface(
interface.id
)

F_1 = +0.5 * mixed_term((u_b + u_t), v_t, n_0) * dInterface(
interface.id
) - 0.5 * mixed_term(v_t, (u_b / K_b - u_t / K_t), n_0) * dInterface(
interface.id
)
F_0 += (
2
* gamma
/ (h_0 + h_1)
* (u_b / K_b - u_t / K_t)
* v_b
* dInterface(interface.id)
)
F_1 += (
-2
* gamma
/ (h_0 + h_1)
* (u_b / K_b - u_t / K_t)
* v_t
* dInterface(interface.id)
)

subdomain_0.F += F_0
subdomain_1.F += F_1

F = sum([subdomain.F for subdomain in self.volume_subdomains])
dus = ufl.TrialFunctions(W)
J = ufl.derivative(F, us, dus)

jit_options={
"cffi_extra_compile_args": ["-O3", "-march=native"],
"cffi_libraries": ["m"],
}
# compile jacobian (J) and residual (F)
self.forms = dolfinx.fem.form(F,
entity_maps=entity_maps,,jit_options=jit_options
)
self.J = dolfinx.fem.form(
J,
entity_maps=entity_maps,
jit_options=jit_options
)
```