From 736bdbbbc060c729979bb480973590d991acc4d3 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 15 Jun 2006 13:26:21 +0000
Subject: [PATCH] separate computation of Darboux vector into separate method
[[Imported from SVN: r840]]
---
src/rodassembler.cc | 30 +++++++-----------------------
src/rodassembler.hh | 33 ++++++++++++++++++++++++++++++++-
2 files changed, 39 insertions(+), 24 deletions(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index ba4d125c..30e1a905 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -316,10 +316,8 @@ getLocalMatrix( EntityType &entity,
hatq_s[2] = localSolution[0].q[2]*shapeGrad[0] + localSolution[1].q[2]*shapeGrad[1];
hatq_s[3] = localSolution[0].q[3]*shapeGrad[0] + localSolution[1].q[3]*shapeGrad[1];
- FieldVector<double,3> u; // The Darboux vector
- u[0] = 2 * ( hatq[3]*hatq_s[0] + hatq[2]*hatq_s[1] - hatq[1]*hatq_s[2] - hatq[0]*hatq_s[3]);
- u[1] = 2 * (-hatq[2]*hatq_s[0] + hatq[3]*hatq_s[1] + hatq[0]*hatq_s[2] - hatq[1]*hatq_s[3]);
- u[2] = 2 * ( hatq[1]*hatq_s[0] - hatq[0]*hatq_s[1] + hatq[3]*hatq_s[2] - hatq[2]*hatq_s[3]);
+ // The Darboux vector
+ FieldVector<double,3> u = darboux(hatq, hatq_s);
// Contains \partial q / \partial v^i_j at v = 0
double dum_dvij[3][2][3];
@@ -503,10 +501,8 @@ assembleGradient(const std::vector<Configuration>& sol,
hatq_s[2] = localSolution[0].q[2]*shapeGrad[0] + localSolution[1].q[2]*shapeGrad[1];
hatq_s[3] = localSolution[0].q[3]*shapeGrad[0] + localSolution[1].q[3]*shapeGrad[1];
- FieldVector<double,3> u; // The Darboux vector
- u[0] = 2 * ( hatq[3]*hatq_s[0] + hatq[2]*hatq_s[1] - hatq[1]*hatq_s[2] - hatq[0]*hatq_s[3]);
- u[1] = 2 * (-hatq[2]*hatq_s[0] + hatq[3]*hatq_s[1] + hatq[0]*hatq_s[2] - hatq[1]*hatq_s[3]);
- u[2] = 2 * ( hatq[1]*hatq_s[0] - hatq[0]*hatq_s[1] + hatq[3]*hatq_s[2] - hatq[2]*hatq_s[3]);
+ // The Darboux vector
+ FieldVector<double,3> u = darboux(hatq, hatq_s);
// Contains \partial q / \partial v^i_j at v = 0
Quaternion<double> dq_dvij[2][3];
@@ -740,12 +736,7 @@ computeEnergy(const std::vector<Configuration>& sol) const
#if 0
// Part II: the bending and twisting energy
- FieldVector<double,3> u; // The Darboux vector
- u[0] = 2 * ( q[3]*q_s[0] + q[2]*q_s[1] - q[1]*q_s[2] - q[0]*q_s[3]);
- u[1] = 2 * (-q[2]*q_s[0] + q[3]*q_s[1] + q[0]*q_s[2] - q[1]*q_s[3]);
- u[2] = 2 * ( q[1]*q_s[0] - q[0]*q_s[1] + q[3]*q_s[2] - q[2]*q_s[3]);
-
- //std::cout << "Darboux vector : " << u << std::endl;
+ FieldVector<double,3> u = darboux(q, q_s); // The Darboux vector
energy += weight * 0.5 * (K1*u[0]*u[0] + K2*u[1]*u[1] + K3*u[2]*u[2]);
#endif
@@ -829,11 +820,7 @@ computeEnergy(const std::vector<Configuration>& sol) const
#endif
// Part II: the bending and twisting energy
- FieldVector<double,3> u; // The Darboux vector
- u[0] = 2 * ( q[3]*q_s[0] + q[2]*q_s[1] - q[1]*q_s[2] - q[0]*q_s[3]);
- u[1] = 2 * (-q[2]*q_s[0] + q[3]*q_s[1] + q[0]*q_s[2] - q[1]*q_s[3]);
- u[2] = 2 * ( q[1]*q_s[0] - q[0]*q_s[1] + q[3]*q_s[2] - q[2]*q_s[3]);
-
+ FieldVector<double,3> u = darboux(q, q_s); // The Darboux vector
//std::cout << "Darboux vector : " << u << std::endl;
energy += weight * 0.5 * (K1*u[0]*u[0] + K2*u[1]*u[1] + K3*u[2]*u[2]);
@@ -958,10 +945,7 @@ getStrain(const std::vector<Configuration>& sol,
// Part II: the Darboux vector
- FieldVector<double,3> u;
- u[0] = 2 * ( q[3]*q_s[0] + q[2]*q_s[1] - q[1]*q_s[2] - q[0]*q_s[3]);
- u[1] = 2 * (-q[2]*q_s[0] + q[3]*q_s[1] + q[0]*q_s[2] - q[1]*q_s[3]);
- u[2] = 2 * ( q[1]*q_s[0] - q[0]*q_s[1] + q[3]*q_s[2] - q[2]*q_s[3]);
+ FieldVector<double,3> u = darboux(q, q_s);
// Sum it all up
strain[elementIdx][0] += weight * v[0];
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 4a3316ef..2204461d 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -58,6 +58,28 @@ namespace Dune
A3 = a3;
}
+ /** \brief Set shape constants and material parameters
+ \param A The rod section area
+ \param J1, J2 The geometric moments (Flächenträgheitsmomente)
+ \param E Young's modulus
+ \param nu Poisson number
+ */
+ void setShapeAndMaterial(double A, double J1, double J2, double E, double nu)
+ {
+ // shear modulus
+ double G = E/(2+2*nu);
+
+ K1 = E * J1;
+ K2 = E * J2;
+ K3 = G * (J1 + J2);
+
+ A1 = G * A;
+ A2 = G * A;
+ A3 = E * A;
+
+ printf("%g %g %g %g %g %g\n", K1, K2, K3, A1, A2, A3);
+ //exit(0);
+ }
/** \brief Assemble the tangent stiffness matrix and the right hand side
*/
@@ -107,7 +129,16 @@ namespace Dune
return r;
}
-
+ template <class T>
+ static FieldVector<T,3> darboux(const Quaternion<T>& q, const FieldVector<T,4>& q_s)
+ {
+ FieldVector<double,3> u; // The Darboux vector
+ u[0] = 2 * ( q[3]*q_s[0] + q[2]*q_s[1] - q[1]*q_s[2] - q[0]*q_s[3]);
+ u[1] = 2 * (-q[2]*q_s[0] + q[3]*q_s[1] + q[0]*q_s[2] - q[1]*q_s[3]);
+ u[2] = 2 * ( q[1]*q_s[0] - q[0]*q_s[1] + q[3]*q_s[2] - q[2]*q_s[3]);
+
+ return u;
+ }
}; // end class
--
GitLab