From 1d18d3b63465ff53cdeabd6d4bad0feb7e5bb229 Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Fri, 9 Mar 2018 16:17:26 +0100
Subject: [PATCH] Add method to compute the stress

---
 dune/gfe/rodlocalstiffness.hh | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/dune/gfe/rodlocalstiffness.hh b/dune/gfe/rodlocalstiffness.hh
index 8d268fbd..25328d86 100644
--- a/dune/gfe/rodlocalstiffness.hh
+++ b/dune/gfe/rodlocalstiffness.hh
@@ -117,6 +117,10 @@ public:
                                            const Entity& element,
                                            const Dune::FieldVector<double,1>& pos) const;
 
+    Dune::FieldVector<RT, 6> getStress(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
+                                           const Entity& element,
+                                           const Dune::FieldVector<double,1>& pos) const;
+
 protected:
 
     void getLocalReferenceConfiguration(const Entity& element,
@@ -502,6 +506,28 @@ getStrain(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
 
     return strain;
 }
+template <class GridType, class RT>
+Dune::FieldVector<RT, 6> RodLocalStiffness<GridType, RT>::
+getStress(const std::vector<RigidBodyMotion<RT,3> >& localSolution,
+              const Entity& element,
+                        const Dune::FieldVector<DT, 1>& pos) const
+{
+    const auto& indexSet = gridView_.indexSet();
+    std::vector<TargetSpace> localRefConf = {referenceConfiguration_[indexSet.subIndex(element, 0, 1)],
+                                             referenceConfiguration_[indexSet.subIndex(element, 1, 1)]};
+
+    auto&& strain = getStrain(localSolution, element, pos);
+    auto&& referenceStrain = getStrain(localRefConf, element, pos);
+
+    Dune::FieldVector<RT, 6> stress;
+    for (int i=0; i < dim; i++)
+        stress[i] = (strain[i] - referenceStrain[i]) * A_[i];
+
+    for (int i=0; i < dim; i++)
+        stress[i+3] = (strain[i+3] - referenceStrain[i+3]) * K_[i];
+    return stress;
+}
+
 
 template <class GridType, class RT>
 void RodLocalStiffness<GridType, RT>::
-- 
GitLab