diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index 81717e62bdb9aab35849a9efb4a2f8fd835766c2..d72a2655b26885b1510b694d02b03611bc2b8444 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -479,20 +479,21 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
-        Rd momentum_fluxes       = zero;
-        double energy_fluxes     = 0;
-        Rdxd cauchy_green_fluxes = zero;
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        Rdxd gradv           = zero;
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          const NodeId r   = cell_nodes[R];
-          const Rdxd gradv = tensorProduct(ur[r], Cjr(j, R));
+          const NodeId r = cell_nodes[R];
+          gradv += tensorProduct(ur[r], Cjr(j, R));
           momentum_fluxes += Fjr(j, R);
           energy_fluxes += dot(Fjr(j, R), ur[r]);
-          cauchy_green_fluxes += gradv * CG[j] + CG[j] * transpose(gradv);
         }
-        const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+        const Rdxd cauchy_green_fluxes = gradv * CG[j] + CG[j] * transpose(gradv);
+        const double dt_over_Mj        = dt / (rho[j] * Vj[j]);
+        const double dt_over_Vj        = dt / Vj[j];
         new_u[j] += dt_over_Mj * momentum_fluxes;
         new_E[j] += dt_over_Mj * energy_fluxes;
-        new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
+        new_CG[j] += dt_over_Vj * cauchy_green_fluxes;
         new_CG[j] += transpose(new_CG[j]);
         new_CG[j] *= 0.5;
       });