diff --git a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
index 2546dd2..bb951ea 100644
--- a/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_diffusion/transport_program.cpp
@@ -1,236 +1,239 @@
 #include "transport_program.hpp"
 #include "transport_parameters.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 #include <iostream>
 
 #define EPS_J 1e-8
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace siasaturated {
 
 index_t SaturatedDiffusionProgram::get_tot_ndf() const {
     return m_mesh->nb_nodes()*get_ndf();
 }
 
 SaturatedDiffusionProgram::SaturatedDiffusionProgram(
         std::shared_ptr<mesh::Mesh1D> the_mesh,
         std::shared_ptr<database::DataContainer> the_data,
         std::shared_ptr<SaturatedDiffusionTransportParameters> the_param
         ):
     m_mesh(the_mesh),
     m_data(the_data),
     m_param(the_param),
     m_external_flow(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1))),
     m_internal_flow(Vector::Zero(the_mesh->nb_nodes()*(the_data->nb_component-1)))
 {}
 
 
 void SaturatedDiffusionProgram::number_equations(std::vector<index_t> list_fixed_nodes
         )
 {
     m_ideq = Eigen::VectorXd::Zero(get_ndf()*m_mesh->nb_nodes());
     for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
     {
         for (index_t component: m_data->range_aqueous_component())
         {
             m_ideq(get_dof_component(*it, component)) = no_equation;
         }
     }
 
     index_t neq = 0;
     for (index_t dof=0; dof<m_ideq.rows(); ++dof)
     {
         if (m_ideq(dof) != no_equation)
         {
             m_ideq(dof) = neq;
             ++neq;
         }
     }
     m_neq = neq;
 }
 
 void SaturatedDiffusionProgram::compute_residuals(
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual
         )
 {
     m_internal_flow.setZero();
     residual = Eigen::VectorXd::Zero(get_neq());
     for (index_t element: m_mesh->range_elements())
     {
         element_residual_transport(element, displacement, velocity, residual);
     }
 }
 
 
 void SaturatedDiffusionProgram::element_residual_transport(
         index_t element,
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual)
 {
     Eigen::VectorXd element_residual(2);
     for (index_t component: m_data->range_aqueous_component())
     {
         element_residual_transport_component(element, component, displacement, velocity, element_residual);
 
         for (index_t en: m_mesh->range_nen())
         {
             const index_t node = m_mesh->get_node(element, en);
             const index_t id = m_ideq(get_dof_component(node, component));
             if (id != no_equation) {residual(id) += element_residual(en);}
         }
     }
 }
 
 void SaturatedDiffusionProgram::element_residual_transport_component(
         index_t element,
         index_t component,
         const Vector &displacement,
         const Vector &velocity,
         Vector& element_residual)
 {
     Eigen::Matrix<scalar_t, 2, 2> jacob;
     Eigen::Matrix<scalar_t, 2, 1> evelocity, edisplacement;
-    scalar_t mass_coeff = -(m_param->density_water()*
-                          m_mesh->get_volume_element(element)/2.0);
+    //scalar_t mass_coeff = -(m_param->density_water()*
+    //                      m_mesh->get_volume_element(element)/2.0);
+
+    const scalar_t mass_coeff_0 = -m_param->density_water()*m_mesh->get_volume_cell_element(element, 0);
+    const scalar_t mass_coeff_1 = -m_param->density_water()*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 porosity = (m_param->porosity(node_0) +
                                m_param->porosity(node_1) )/2.0;
     const scalar_t diff_coeff = 1.0/(0.5/m_param->diffusion_coefficient(node_0) +
                                      0.5/m_param->diffusion_coefficient(node_1));
 
     const index_t dof_0 = get_dof_component(node_0, component);
     const index_t dof_1 = get_dof_component(node_1, component);
     scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
-                           * porosity
+                           //* porosity
                            * m_param->density_water()
                            * diff_coeff
                            );
     jacob << 1.0, -1.0, -1.0, 1.0;
     jacob *= flux_coeff;
     evelocity << velocity(dof_0)*
-                     mass_coeff*m_param->porosity(node_0),
+                     mass_coeff_0*m_param->porosity(node_0),
             velocity(dof_1)*
-                     mass_coeff*m_param->porosity(node_1);
+                     mass_coeff_1*m_param->porosity(node_1);
 
     edisplacement << displacement(dof_0),
                      displacement(dof_1);
 
     element_residual = jacob*edisplacement;
 
     m_internal_flow(dof_0) += element_residual(0);
     m_internal_flow(dof_1) += element_residual(1);
 
     for (index_t en: m_mesh->range_nen())
     {
         element_residual(en) += evelocity(en);
         element_residual(en) += 1e6*(m_mesh->get_volume_element(element)/2.0
                                  *external_flow(m_mesh->get_node(element, en), component));
     }
 }
 
 
 void SaturatedDiffusionProgram::compute_jacobian(
         Vector& displacement,
         Vector& velocity,
         Eigen::SparseMatrix<scalar_t>& jacobian,
         scalar_t alphadt
         )
 {
 
     list_triplet_t jacob;
     const index_t ncomp = m_data->nb_component-1;
     const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
     jacob.reserve(estimation);
     for (index_t element: m_mesh->range_elements())
     {
         element_jacobian_transport(element, displacement, velocity, jacob, alphadt);
     }
     jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
 // Transport
 // =========
 
 void SaturatedDiffusionProgram::element_jacobian_transport(
         index_t element,
         Vector& displacement,
         Vector& velocity,
         list_triplet_t& jacobian,
         scalar_t alphadt)
 {
     for (index_t component: m_data->range_aqueous_component())
     {
         Eigen::VectorXd element_residual_orig(Eigen::VectorXd::Zero(2));
         element_residual_transport_component(element, component, displacement, velocity, element_residual_orig);
 
         for (index_t en: m_mesh->range_nen())
         {
             Eigen::VectorXd element_residual(Eigen::VectorXd::Zero(2));
 
             const index_t node = m_mesh->get_node(element, en);
             const index_t dof = get_dof_component(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_J*std::abs(tmp_v);
             if (h < 1e-4*EPS_J) h = EPS_J;
             velocity(dof) = tmp_v + h;
             h = velocity(dof) - tmp_v;
 
             displacement(dof) = tmp_d + alphadt*h;
 
             element_residual_transport_component(element, component, displacement, velocity, element_residual);
             velocity(dof) = tmp_v;
             displacement(dof) = tmp_d;
 
             for (index_t enr: m_mesh->range_nen())
             {
                 const index_t noder = m_mesh->get_node(element, enr);
                 const index_t idr = m_ideq(get_dof_component(noder, component));
 
                 if (idr == no_equation) continue;
                 jacobian.push_back(triplet_t(idr, idc, (element_residual(enr) - element_residual_orig(enr))/h));
             }
         }
 
     }
 }
 
 void SaturatedDiffusionProgram::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: m_data->range_aqueous_component())
         {
             const index_t dof = get_dof_component(node, component);
             const index_t id = m_ideq(dof);
             if (id == no_equation) continue;
             velocity(dof) += lambda*update(id);
         }
     }
     displacement = predictor + alpha_dt*velocity;
 }
 
 } // end namespace siasaturated
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp