BUG: c=0 trick for steady state problems

In `HydrogenTransportProblem` we want apply `c = 0` to the formulation in steady state in subdomains where immobile species are not involved in a reaction. We do this here by looking in which subdomains each immobile species isn't defined:

https://github.com/festim-dev/FESTIM/blob/013937a860d72252436118dabcbf75b80ebc5cce/src/festim/hydrogen_transport_problem.py#L683-L699

However, this will apply $c = 0$ in all the subdomains where there is _a_ reaction, even if this reaction doesn't involve this specific species.

These lines are the issue:
https://github.com/festim-dev/FESTIM/blob/013937a860d72252436118dabcbf75b80ebc5cce/src/festim/hydrogen_transport_problem.py#L689-L693

It should check if the species is in the product of the reaction:

```python
for vol in self.volume_subdomains:
# check reactions
for reaction in self.reactions:
if spe in reaction.product: # <-- check if the species is in the reaction
if vol == reaction.volume:
not_defined_in_volume.remove(vol)
```

To reproduce the bug, run this on `fenicsx`:

```python
import festim as F
import numpy as np

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh1D(
np.concatenate([np.linspace(0, 0.5, 10), np.linspace(0.5, 1, 10)])
)

mat = F.Material(D_0=1, E_D=0.1)

vol1 = F.VolumeSubdomain1D(id=1, material=mat, borders=[0, 0.5])
vol2 = F.VolumeSubdomain1D(id=2, material=mat, borders=[0.5, 1])

surface1 = F.SurfaceSubdomain1D(id=4, x=0)
surface2 = F.SurfaceSubdomain1D(id=5, x=1)

my_model.subdomains = [vol1, vol2, surface1, surface2]

mobile = F.Species(name="H", mobile=True)
trapped = F.Species(name="H_trapped", mobile=False)
trapped2 = F.Species(name="H_trapped", mobile=False)
empty = F.ImplicitSpecies(n=0.5, others=[trapped])
empty2 = F.ImplicitSpecies(n=0.5, others=[trapped2])

my_model.species = [mobile, trapped, trapped2]

my_model.reactions = [
F.Reaction(
reactant=[mobile, empty],
product=[trapped],
k_0=1,
E_k=mat.E_D,
p_0=0.1,
E_p=0.87,
volume=vol1,
),
F.Reaction(
reactant=[mobile, empty2],
product=[trapped2],
k_0=1,
E_k=mat.E_D,
p_0=0.1,
E_p=0.87,
volume=vol2,
),
]

my_model.temperature = 600

my_model.boundary_conditions = [
F.FixedConcentrationBC(subdomain=surface1, species=mobile, value=0),
F.FixedConcentrationBC(subdomain=surface2, species=mobile, value=1),
]

my_model.settings = F.Settings(
atol=1e-10,
rtol=1e-10,
transient=False,
)

my_model.initialise()
my_model.run()
```
Produces:

```
Traceback (most recent call last):
File "/home/remidm/FESTIM/mwe_trick.py", line 76, in
my_model.run()
File "/home/remidm/FESTIM/src/festim/problem.py", line 150, in run
self.solver.solve(self.u)
File "/home/remidm/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/nls/petsc.py", line 51, in solve
n, converged = super().solve(u.x.petsc_vec)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 76, Error in external library
```