From 3f4037c885f3f1d2c130ffe95206fa096218c529 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 19 May 2008 09:24:42 +0000
Subject: [PATCH] add method getStress
[[Imported from SVN: r2303]]
---
src/rodassembler.cc | 24 +++++++++++++++++++++++-
src/rodassembler.hh | 3 +++
2 files changed, 26 insertions(+), 1 deletion(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 57dedd6f..37689816 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 4da84973..80beadf3 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 */
--
GitLab