diff --git a/HyperElasticity/beam.pro b/HyperElasticity/beam.pro
index 30568cc67d31dbf6c712ded10a2eee5159f7352e..2d8a7016491663e702ec0052d6a357c2158ee90f 100644
--- a/HyperElasticity/beam.pro
+++ b/HyperElasticity/beam.pro
@@ -273,19 +273,22 @@ Formulation {
       EndFor
     }
     Equation {
+      // store useful intermediate computation at quadrature points
+      Galerkin { [ 0 * SetVariable[ TensorDiag[1., 1., 1.] +
+                                    TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ],
+                                    QuadraturePointIndex[] ]{$arg} , {u~{1}} ] ;
+        In Vol_Mec; Jacobian JVol; Integration GradGrad; }
+
       For i In {1:3}
         For j In {1:3}
-          Galerkin { [ C_Tot~{j}~{i}[ TensorDiag[1.0, 1.0, 1.0] +
-                                      TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ] *
-                       Dof{d u~{i}}, {d u~{j}} ] ;
+          Galerkin { [ SetVariable[ C_Tot~{j}~{i}[ GetVariable[QuadraturePointIndex[]]{$arg} ],
+                                    QuadraturePointIndex[] ]{$ctot} *
+                       Dof{d u~{i}} , {d u~{j}} ] ;
             In Vol_Mec ; Jacobian JVol ; Integration GradGrad ; }
-          Galerkin { [-C_Tot~{j}~{i}[ TensorDiag[1.0, 1.0, 1.0] +
-                                      TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ] *
-                       {d u~{i}}, {d u~{j}} ] ;
+          Galerkin { [ - GetVariable[QuadraturePointIndex[]]{$ctot} * {d u~{i}} , {d u~{j}} ] ;
             In Vol_Mec ; Jacobian JVol ; Integration GradGrad ; }
         EndFor
-        Galerkin { [ P~{i}[ TensorDiag[1., 1., 1.] +
-                            TensorV[{d u~{1}}, {d u~{2}}, {d u~{3}} ] ] , {d u~{i}} ] ;
+        Galerkin { [ P~{i}[ GetVariable[QuadraturePointIndex[]]{$arg} ] , {d u~{i}} ] ;
           In Vol_Mec; Jacobian JVol; Integration GradGrad; }
       EndFor
 
@@ -368,7 +371,7 @@ Resolution {
       Generate[Sys_Mec]; GetResidual[Sys_Mec, $res0]; Evaluate[ $res = $res0, $iter = 0 ];
       Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
       While[$res > tol_abs[] && $res / $res0 > tol_rel[] && $res / $res0 <= 1 && $iter < iter_max]{
-        Solve[Sys_Mec]; Evaluate[ $syscount = $syscount + 1 ];
+        Solve[Sys_Mec];
         Generate[Sys_Mec]; GetResidual[Sys_Mec, $res]; Evaluate[ $iter = $iter + 1 ];
         Print[{$iter, $res, $res / $res0}, Format "Residual %03g: abs %14.12e rel %14.12e"];
       }