diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 57dedd6f268bf084d9570423ced7e3c9311638e0..37689816d08a599559615700b57c4055d861fca0 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -1062,6 +1062,27 @@ getStrain(const std::vector<Configuration>& sol,
 
 }
 
+template <class GridType>
+void RodAssembler<GridType>::
+getStress(const std::vector<Configuration>& sol,
+          Dune::BlockVector<Dune::FieldVector<double, blocksize> >& stress) const
+{
+    // Get the strain
+    getStrain(sol,stress);
+
+    // Get reference strain
+    Dune::BlockVector<Dune::FieldVector<double, blocksize> > referenceStrain;
+    getStrain(referenceConfiguration_, referenceStrain);
+
+    // Linear diagonal constitutive law
+    for (size_t i=0; i<stress.size(); i++) {
+        for (int j=0; j<3; j++) {
+            stress[i][j]   = (stress[i][j]   - referenceStrain[i][j])   * A_[j];
+            stress[i][j+3] = (stress[i][j+3] - referenceStrain[i][j+3]) * K_[j];
+        }
+    }
+}
+
 template <class GridType>
 Dune::FieldVector<double,3> RodAssembler<GridType>::
 getResultantForce(const BoundaryPatch<GridType>& boundary, 
@@ -1142,7 +1163,8 @@ getResultantForce(const BoundaryPatch<GridType>& boundary,
 //             assert( std::abs(localStress[2]-canonicalStress*sol[0].q.director(2)) < 1e-6 );
 
             // Multiply force times boundary normal to get the transmitted force
-            // I am not quite sure why the -1 is there, but it has to be there.
+            /** \todo The minus sign comes from the coupling conditions.  It
+                should really be in the Dirichlet-Neumann code. */
             canonicalStress *= -nIt->unitOuterNormal(FieldVector<double,0>(0))[0];
             canonicalTorque *= -nIt->unitOuterNormal(FieldVector<double,0>(0))[0];
             
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 4da849730cd2b2a41fde3acda78c9b11af3e60a8..80beadf3a228c944800899149df36ae51967e583 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -259,6 +259,9 @@ public:
         void getStrain(const std::vector<Configuration>& sol, 
                        Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const;
 
+        void getStress(const std::vector<Configuration>& sol, 
+                       Dune::BlockVector<Dune::FieldVector<double, blocksize> >& stress) const;
+
         /** \brief Return resultant force across boundary in canonical coordinates 
 
         \note Linear run-time in the size of the grid */