diff --git a/src/reactmicp/systems/saturated_react/transport_program.cpp b/src/reactmicp/systems/saturated_react/transport_program.cpp
index f28b8b8..4a5f114 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.cpp
@@ -1,288 +1,285 @@
 #include "transport_program.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 #include "variables.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 //class SaturatedDiffusion::
 
 
 SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
                    std::vector<index_t> list_fixed_nodes):
     m_mesh(variables->get_mesh()),
     m_variables(variables),
-    m_is_in_residual_computation(false)
+    m_is_in_residual_computation(false),
+    m_ndf(2*variables->nb_component()),
+    m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes())
 {
     number_equations(list_fixed_nodes, {}, {0,});
 }
 
 
 SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
                    std::vector<index_t> list_fixed_nodes,
                    std::map<index_t, index_t>  list_slave_nodes,
                    std::vector<index_t> list_immobile_components):
     m_mesh(variables->get_mesh()),
     m_variables(variables),
-    m_is_in_residual_computation(false)
+    m_is_in_residual_computation(false),
+    m_ndf(2*variables->nb_component()),
+    m_tot_ndf(2*variables->nb_component()*variables->get_mesh()->nb_nodes())
 {
     number_equations(list_fixed_nodes, list_slave_nodes, list_immobile_components);
 }
 
-//! \brief Return the number of degrees of freedom per node
-index_t SaturatedDiffusion::get_ndf() const {return 2*m_variables->nb_component();}
-
-//! \brief Return the total number of degrees of freedom
-index_t SaturatedDiffusion::get_tot_ndf() const {return 2*m_variables->nb_component()*m_variables->get_mesh()->nb_nodes();}
-
-
 void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes,
                                           std::map<index_t, index_t>  list_slave_nodes,
                                           std::vector<index_t> list_immobile_components)
 {
     m_ideq.resizeLike(m_variables->displacement());
     m_ideq.setZero();
     // flag fixed nodes
     for (index_t node: list_fixed_nodes)
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
         }
     }
     // flag slaves nodes
     // we flag them by making their ideq more negative than no_equation
     for (auto slave_pair: list_slave_nodes)
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) = no_equation-1;
         }
     }
     // set equation numbers
     index_t neq = 0;
     for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
     {
         for (index_t component: list_immobile_components)
         {
             m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
             m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
         }
         for (index_t component=0; component<m_variables->nb_component(); ++component)
         {
             index_t dof = m_variables->dof_aqueous_concentration(node, component);
             if (m_ideq(dof) > no_equation) // don't attribute an equation number if it is a slave or a fixed node
             {
                 m_ideq(dof) = neq;
                 ++neq;
             }
             dof = m_variables->dof_solid_concentration(node, component);
             m_ideq(dof) = no_equation;
         }
     }
     // slave nodes
     // attribute the correct equation number
     for (auto slave_pair: list_slave_nodes)
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) =
                     m_ideq(m_variables->dof_aqueous_concentration(slave_pair.second, component));
         }
     }
     m_neq = neq;
 }
 
 void SaturatedDiffusion::compute_residuals(
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual)
 {
     residual = Vector::Zero(get_neq());
     m_is_in_residual_computation = true;
     for (index_t element: m_mesh->range_elements())
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             Eigen::Vector2d element_residual;
             element_residual.setZero();
             residuals_element_component(element, component, displacement, velocity, element_residual);
             for (index_t en=0; en<2; ++en)
             {
                 const index_t node = m_mesh->get_node(element, en);
                 const index_t id = m_ideq(m_variables->dof_aqueous_concentration(node, component));
                 if (id != no_equation) {residual(id) += element_residual(en);}
             }
         }
     }
     m_is_in_residual_computation = false;
 }
 
 void SaturatedDiffusion::residuals_element_component(
         index_t element,
         index_t component,
         const Vector& displacement,
         const Vector& velocity,
         Eigen::Vector2d& element_residual
         )
 {
     const scalar_t mass_coeff_0 = -m_mesh->get_volume_cell_element(element, 0);
     const scalar_t mass_coeff_1 = -m_mesh->get_volume_cell_element(element, 1);
 
     const index_t node_0 = m_mesh->get_node(element, 0);
     const index_t node_1 = m_mesh->get_node(element, 1);
 
     const scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) +
                                      0.5/m_variables->diffusion_coefficient(node_1));
 
 
     scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
                            * diff_coeff
                            ) ;
 
     const index_t dof_0 = m_variables->dof_aqueous_concentration(node_0, component);
     const index_t dof_1 = m_variables->dof_aqueous_concentration(node_1, component);
 
     // diffusion
     scalar_t flux_diffusion = flux_coeff*(displacement(dof_0) - displacement(dof_1));
 
     element_residual(0) = flux_diffusion;
     element_residual(1) = - flux_diffusion;
 
     // advection
 
     if (m_variables->fluid_velocity(element) != 0.0)
     {
         scalar_t flux_advection = (m_mesh->get_face_area(element))
                 *m_variables->fluid_velocity(element);
         if (m_variables->fluid_velocity(element) > 0)
         {
             flux_advection *= (displacement(dof_0) - displacement(dof_1));
             element_residual(1) += flux_advection;
         }
         else
         {
             flux_advection *= (displacement(dof_1) - displacement(dof_0));
             element_residual(0) -= flux_advection;
         }
     }
 
     if (m_is_in_residual_computation)
     {
         m_variables->aqueous_concentration(node_0, component,
                            m_variables->transport_rate()) += element_residual(0);
         m_variables->aqueous_concentration(node_1, component,
                            m_variables->transport_rate()) += element_residual(1);
     }
     // velocity
     element_residual(0) += mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0)
                                          +m_variables->vel_porosity(node_0)*displacement(dof_0));
     element_residual(1) += mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1)
                                          + m_variables->vel_porosity(node_1)*displacement(dof_1));
     // external rate
     element_residual(0) +=  mass_coeff_0*m_variables->solid_concentration(node_0, component,
                 m_variables->chemistry_rate());
     element_residual(1) +=  mass_coeff_1*m_variables->solid_concentration(node_1, component,
                 m_variables->chemistry_rate());
 }
 
 
 void SaturatedDiffusion::compute_jacobian(
         Vector& displacement,
         Vector& velocity,
         Eigen::SparseMatrix<scalar_t>& jacobian,
         scalar_t alphadt
         )
 {
     dfpm::list_triplet_t jacob;
     const index_t ncomp = m_variables->nb_component();
     const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
     jacob.reserve(estimation);
     for (index_t element: m_mesh->range_elements())
     {
         jacobian_element(element, displacement, velocity, jacob, alphadt);
     }
     jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
 void SaturatedDiffusion::jacobian_element(
         index_t element,
         Vector& displacement,
         Vector& velocity,
         dfpm::list_triplet_t& jacobian,
         scalar_t alphadt)
 {
     for (index_t component=1; component<m_variables->nb_component(); ++component)
     {
         Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero());
         residuals_element_component(element, component, displacement, velocity, element_residual_orig);
 
         for (index_t en=0; en<2; ++en)
         {
             Eigen::Vector2d element_residual(Eigen::Vector2d::Zero());
 
             const index_t node = m_mesh->get_node(element, en);
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t idc = m_ideq(dof);
 
             if (idc == no_equation) continue;
 
             const scalar_t tmp_v = velocity(dof);
             const scalar_t tmp_d = displacement(dof);
 
             scalar_t h = eps_jacobian*std::abs(tmp_v);
             if (h < 1e-4*eps_jacobian) h = eps_jacobian;
             velocity(dof) = tmp_v + h;
             h = velocity(dof) - tmp_v;
 
             displacement(dof) = tmp_d + alphadt*h;
 
             residuals_element_component(element, component, displacement, velocity, element_residual);
             velocity(dof) = tmp_v;
             displacement(dof) = tmp_d;
 
             for (index_t enr=0; enr<2; ++enr)
             {
                 const index_t noder = m_mesh->get_node(element, enr);
                 const index_t idr = m_ideq(m_variables->dof_aqueous_concentration(noder, component));
 
                 if (idr == no_equation) continue;
                 jacobian.push_back(dfpm::triplet_t(
                               idr,
                               idc,
                               (element_residual(enr) - element_residual_orig(enr))/h
                               ));
             }
         }
     }
 }
 
 //! \brief Update the solutions
 void SaturatedDiffusion::update_solution(const Vector& update,
                      scalar_t lambda,
                      scalar_t alpha_dt,
                      Vector& predictor,
                      Vector& displacement,
                      Vector& velocity)
 {
     for (index_t node: m_mesh->range_nodes())
     {
         for (index_t component=1; component< m_variables->nb_component(); ++component)
         {
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t id = m_ideq(dof);
             if (id == no_equation) continue;
             velocity(dof) += lambda*update(id);
         }
     }
     //displacement = m_variables->predictor() + alpha_dt*velocity;
     displacement = predictor + alpha_dt*velocity;
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/transport_program.hpp b/src/reactmicp/systems/saturated_react/transport_program.hpp
index 166eb6e..c3c1b58 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.hpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.hpp
@@ -1,108 +1,111 @@
 #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP
 
 #include "dfpmsolver/parabolic_program.hpp"
 
 #include "variablesfwd.hpp"
 #include "dfpm/meshes/mesh1dfwd.hpp"
 #include "dfpm/types.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 class SaturatedDiffusion:  public dfpmsolver::ParabolicProgram<SaturatedDiffusion>
 {
 public:
 
     SaturatedDiffusion(SaturatedVariablesPtr variables,
                        std::vector<index_t> list_fixed_nodes,
                        std::map<index_t, index_t>  list_slave_nodes,
                        std::vector<index_t> list_immobile_components);
 
     SaturatedDiffusion(SaturatedVariablesPtr variables,
                        std::vector<index_t> list_fixed_nodes);
 
     //! \brief Return the number of equations
     index_t get_neq() const {return m_neq;}
     //! \brief Return the number of degrees of freedom per node
-    index_t get_ndf() const;
+    index_t get_ndf() const {return m_ndf;}
     //! \brief Return the total number of degrees of freedom
-    index_t get_tot_ndf() const;
+    index_t get_tot_ndf() const  {return m_tot_ndf;}
 
     //! \brief Method to update the variables
     void set_variables(SaturatedVariablesPtr variables) {m_variables = variables;}
 
     //! \brief Return the id of the equation corresponding to the degree of freedom 'id_dof'
     //!
     //! Return 'no_equation' if no equation exist
     index_t id_equation(index_t id_dof) const {return m_ideq(id_dof);}
 
     //! \brief Compute the residuals
     void compute_residuals(const Vector& displacement,
                            const Vector& velocity,
                            Vector& residual
                            );
     //! Compute the residuals inside 'element' for 'component'
     void residuals_element_component(
             index_t element,
             index_t component,
             const Vector& displacement,
             const Vector& velocity,
             Eigen::Vector2d& element_residual
             );
 
     //! \brief Compute the jacobian
     void compute_jacobian(Vector& displacement,
                           Vector& velocity,
                           Eigen::SparseMatrix<scalar_t>& jacobian,
                           scalar_t alphadt
                           );
     //! \brief Compute the contribution of 'element' in the jacobian
     void jacobian_element(
             index_t element,
             Vector& displacement,
             Vector& velocity,
             dfpm::list_triplet_t& jacobian,
             scalar_t alphadt);
 
     //! \brief Update the solutions
     void update_solution(const Vector& update,
                          scalar_t lambda,
                          scalar_t alpha_dt,
                          Vector& predictor,
                          Vector& displacement,
                          Vector& velocity);
 
     //! \brief Apply boundary conditions to the velocity vector
     //!
     //! by default do nothing.
     void apply_bc(scalar_t dt,
                   const Vector& displacement,
                   Vector& velocity)
     {}
 
 
 private:
     // number the equations
     void number_equations(std::vector<index_t> list_fixed_nodes,
                           std::map<index_t, index_t>  list_slave_nodes, std::vector<index_t> list_immobile_components);
 
 
     index_t m_neq;
+    index_t m_ndf;
+    index_t m_tot_ndf;
+
     Eigen::Matrix<index_t, Eigen::Dynamic, 1> m_ideq;
 
     mesh::Mesh1DPtr m_mesh;
     SaturatedVariablesPtr m_variables;
 
     bool m_is_in_residual_computation;
 
 };
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP