diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 074a14f6331ae0c2821c863fd07061706bf35f14..2827907480cf10a382ab0dd67e6fb9fed491f928 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -7,6 +7,156 @@
#include <dune/disc/shapefunctions/lagrangeshapefunctions.hh>
+
+template <class GridType, class RT>
+RT RodLocalStiffness<GridType, RT>::
+energy(const EntityPointer& element,
+ const std::vector<Configuration>& localSolution,
+ const std::vector<Configuration>& localReferenceConfiguration,
+ int k)
+{
+ RT energy = 0;
+
+ //const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet();
+
+ double elementEnergy = 0;
+
+ // Extract local solution on this element
+ const Dune::LagrangeShapeFunctionSet<double, double, 1> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, 1>::general(element->type(), k);
+ int numOfBaseFct = baseSet.size();
+
+ // ///////////////////////////////////////////////////////////////////////////////
+ // The following two loops are a reduced integration scheme. We integrate
+ // the transverse shear and extensional energy with a first-order quadrature
+ // formula, even though it should be second order. This prevents shear-locking
+ // ///////////////////////////////////////////////////////////////////////////////
+
+ const int shearingPolOrd = 2;
+ const Dune::QuadratureRule<double, 1>& shearingQuad
+ = Dune::QuadratureRules<double, 1>::rule(element->type(), shearingPolOrd);
+
+ for (size_t pt=0; pt<shearingQuad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<double,1>& quadPos = shearingQuad[pt].position();
+
+ const double integrationElement = element->geometry().integrationElement(quadPos);
+
+ double weight = shearingQuad[pt].weight() * integrationElement;
+
+ Dune::FieldVector<double,6> strain = getStrain(localSolution, element, quadPos);
+
+ // The reference strain
+ Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+
+ for (int i=0; i<3; i++) {
+ energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
+ elementEnergy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
+ }
+ }
+
+ // Get quadrature rule
+ const int polOrd = 2;
+ const Dune::QuadratureRule<double, 1>& bendingQuad = Dune::QuadratureRules<double, 1>::rule(element->type(), polOrd);
+
+ for (size_t pt=0; pt<bendingQuad.size(); pt++) {
+
+ // Local position of the quadrature point
+ const Dune::FieldVector<double,1>& quadPos = bendingQuad[pt].position();
+
+ double weight = bendingQuad[pt].weight() * element->geometry().integrationElement(quadPos);
+
+ Dune::FieldVector<double,6> strain = getStrain(localSolution, element, quadPos);
+
+ // The reference strain
+ Dune::FieldVector<double,6> referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+
+ // Part II: the bending and twisting energy
+ for (int i=0; i<3; i++) {
+ energy += weight * 0.5 * K_[i] * (strain[i+3] - referenceStrain[i+3]) * (strain[i+3] - referenceStrain[i+3]);
+ elementEnergy += weight * 0.5 * K_[i] * (strain[i+3] - referenceStrain[i+3]) * (strain[i+3] - referenceStrain[i+3]);
+ }
+
+ }
+
+ return energy;
+}
+
+template <class GridType, class RT>
+Dune::FieldVector<double, 6> RodLocalStiffness<GridType, RT>::
+getStrain(const std::vector<Configuration>& localSolution,
+ const EntityPointer& element,
+ const Dune::FieldVector<double,1>& pos) const
+{
+ if (!element->isLeaf())
+ DUNE_THROW(Dune::NotImplemented, "Only for leaf elements");
+
+ assert(localSolution.size() == 2);
+
+ // Strain defined on each element
+ Dune::FieldVector<double, 6> strain(0);
+
+ // Extract local solution on this element
+ const Dune::LagrangeShapeFunctionSet<double, double, 1> & baseSet
+ = Dune::LagrangeShapeFunctions<double, double, 1>::general(element->type(), 1);
+ int numOfBaseFct = baseSet.size();
+
+ const Dune::FieldMatrix<double,1,1>& inv = element->geometry().jacobianInverseTransposed(pos);
+
+ // ///////////////////////////////////////
+ // Compute deformation gradient
+ // ///////////////////////////////////////
+ Dune::FieldVector<double,1> shapeGrad[numOfBaseFct];
+
+ for (int dof=0; dof<numOfBaseFct; dof++) {
+
+ for (int i=0; i<1; i++)
+ shapeGrad[dof][i] = baseSet[dof].evaluateDerivative(0,i,pos);
+
+ // multiply with jacobian inverse
+ Dune::FieldVector<double,1> tmp(0);
+ inv.umv(shapeGrad[dof], tmp);
+ shapeGrad[dof] = tmp;
+
+ }
+
+ // //////////////////////////////////
+ // Interpolate
+ // //////////////////////////////////
+
+ Dune::FieldVector<double,3> r_s;
+ for (int i=0; i<3; i++)
+ r_s[i] = localSolution[0].r[i]*shapeGrad[0][0] + localSolution[1].r[i]*shapeGrad[1][0];
+
+ // Interpolate the rotation at the quadrature point
+ Quaternion<double> q = Quaternion<double>::interpolate(localSolution[0].q, localSolution[1].q, pos);
+
+ // Get the derivative of the rotation at the quadrature point by interpolating in $H$
+ Quaternion<double> q_s = Quaternion<double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
+ pos, 1/shapeGrad[1]);
+
+ // /////////////////////////////////////////////
+ // Sum it all up
+ // /////////////////////////////////////////////
+
+ // Part I: the shearing and stretching strain
+ strain[0] = r_s * q.director(0); // shear strain
+ strain[1] = r_s * q.director(1); // shear strain
+ strain[2] = r_s * q.director(2); // stretching strain
+
+ // Part II: the Darboux vector
+
+ Dune::FieldVector<double,3> u = darboux(q, q_s);
+
+ strain[3] = u[0];
+ strain[4] = u[1];
+ strain[5] = u[2];
+
+ return strain;
+}
+
+
template <class GridType>
void Dune::RodAssembler<GridType>::
getNeighborsPerVertex(MatrixIndexSet& nb) const
@@ -377,6 +527,208 @@ getLocalMatrix( EntityPointer &entity,
}
+template <class GridType>
+void Dune::RodAssembler<GridType>::
+assembleMatrixFD(const std::vector<Configuration>& sol,
+ BCRSMatrix<MatrixBlock>& matrix) const
+{
+ double eps = 1e-2;
+
+ typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >::row_type::iterator ColumnIterator;
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = sol;
+ std::vector<Configuration> backwardSolution = sol;
+
+ std::vector<Configuration> forwardForwardSolution = sol;
+ std::vector<Configuration> forwardBackwardSolution = sol;
+ std::vector<Configuration> backwardForwardSolution = sol;
+ std::vector<Configuration> backwardBackwardSolution = sol;
+
+ // ////////////////////////////////////////////////////
+ // Create local assembler
+ // ////////////////////////////////////////////////////
+
+ Dune::array<double,3> K = {K_[0], K_[1], K_[2]};
+ Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
+ RodLocalStiffness<GridType,double> localStiffness(K, A);
+
+ // //////////////////////////////////////////////////////////////////
+ // Store pointers to all elements so we can access them by index
+ // //////////////////////////////////////////////////////////////////
+
+ const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet();
+
+ std::vector<EntityPointer> elements(grid_->size(0),grid_->template leafbegin<0>());
+
+ typename GridType::template Codim<0>::LeafIterator eIt = grid_->template leafbegin<0>();
+ typename GridType::template Codim<0>::LeafIterator eEndIt = grid_->template leafend<0>();
+
+ for (; eIt!=eEndIt; ++eIt)
+ elements[indexSet.index(*eIt)] = eIt;
+
+ std::vector<Configuration> localReferenceConfiguration(2);
+ std::vector<Configuration> localSolution(2);
+
+ // ///////////////////////////////////////////////////////////////
+ // Loop over all blocks of the outer matrix
+ // ///////////////////////////////////////////////////////////////
+ for (int i=0; i<matrix.N(); i++) {
+
+ ColumnIterator cIt = matrix[i].begin();
+ ColumnIterator cEndIt = matrix[i].end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // ////////////////////////////////////////////////////////////////////////////
+ // Compute a finite-difference approximation of this hessian matrix block
+ // ////////////////////////////////////////////////////////////////////////////
+
+ for (int j=0; j<6; j++) {
+
+ for (int k=0; k<6; k++) {
+
+ if (i==cIt.index() && j==k) {
+
+ double forwardEnergy = 0;
+ double solutionEnergy = 0;
+ double backwardEnergy = 0;
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+ if (i != grid_->size(1)-1) {
+
+ localReferenceConfiguration[0] = referenceConfiguration_[i];
+ localReferenceConfiguration[1] = referenceConfiguration_[i+1];
+
+ localSolution[0] = forwardSolution[i];
+ localSolution[1] = forwardSolution[i+1];
+
+ forwardEnergy += localStiffness.energy(elements[i], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = sol[i];
+ localSolution[1] = sol[i+1];
+
+ solutionEnergy += localStiffness.energy(elements[i], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = backwardSolution[i];
+ localSolution[1] = backwardSolution[i+1];
+
+ backwardEnergy += localStiffness.energy(elements[i], localSolution, localReferenceConfiguration);
+
+ }
+
+ if (i != 0) {
+
+ localReferenceConfiguration[0] = referenceConfiguration_[i-1];
+ localReferenceConfiguration[1] = referenceConfiguration_[i];
+
+ localSolution[0] = forwardSolution[i-1];
+ localSolution[1] = forwardSolution[i];
+
+ forwardEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = sol[i-1];
+ localSolution[1] = sol[i];
+
+ solutionEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = backwardSolution[i-1];
+ localSolution[1] = backwardSolution[i];
+
+ backwardEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+
+ }
+
+ // Second derivative
+ (*cIt)[j][k] = (forwardEnergy - 2*solutionEnergy + backwardEnergy) / (eps*eps);
+
+ forwardSolution[i] = sol[i];
+ backwardSolution[i] = sol[i];
+
+ } else {
+
+ double forwardForwardEnergy = 0;
+ double forwardBackwardEnergy = 0;
+ double backwardForwardEnergy = 0;
+ double backwardBackwardEnergy = 0;
+
+ infinitesimalVariation(forwardForwardSolution[i], eps, j);
+ infinitesimalVariation(forwardForwardSolution[cIt.index()], eps, k);
+ infinitesimalVariation(forwardBackwardSolution[i], eps, j);
+ infinitesimalVariation(forwardBackwardSolution[cIt.index()], -eps, k);
+ infinitesimalVariation(backwardForwardSolution[i], -eps, j);
+ infinitesimalVariation(backwardForwardSolution[cIt.index()], eps, k);
+ infinitesimalVariation(backwardBackwardSolution[i], -eps, j);
+ infinitesimalVariation(backwardBackwardSolution[cIt.index()],-eps, k);
+
+ std::set<int> elementsInvolved;
+ if (i>0)
+ elementsInvolved.insert(i-1);
+ if (i<grid_->size(1)-1)
+ elementsInvolved.insert(i);
+ if (cIt.index()>0)
+ elementsInvolved.insert(cIt.index()-1);
+ if (cIt.index()<grid_->size(1)-1)
+ elementsInvolved.insert(cIt.index());
+
+ for (typename std::set<int>::iterator it = elementsInvolved.begin();
+ it != elementsInvolved.end();
+ ++it) {
+
+ localReferenceConfiguration[0] = referenceConfiguration_[*it];
+ localReferenceConfiguration[1] = referenceConfiguration_[*it+1];
+
+ localSolution[0] = forwardForwardSolution[*it];
+ localSolution[1] = forwardForwardSolution[*it+1];
+
+ forwardForwardEnergy += localStiffness.energy(elements[*it], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = forwardBackwardSolution[*it];
+ localSolution[1] = forwardBackwardSolution[*it+1];
+
+ forwardBackwardEnergy += localStiffness.energy(elements[*it], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = backwardForwardSolution[*it];
+ localSolution[1] = backwardForwardSolution[*it+1];
+
+ backwardForwardEnergy += localStiffness.energy(elements[*it], localSolution, localReferenceConfiguration);
+
+ localSolution[0] = backwardBackwardSolution[*it];
+ localSolution[1] = backwardBackwardSolution[*it+1];
+
+ backwardBackwardEnergy += localStiffness.energy(elements[*it], localSolution, localReferenceConfiguration);
+
+
+ }
+
+ (*cIt)[j][k] = (forwardForwardEnergy + backwardBackwardEnergy
+ - forwardBackwardEnergy - backwardForwardEnergy) / (4*eps*eps);
+
+ forwardForwardSolution[i] = sol[i];
+ forwardForwardSolution[cIt.index()] = sol[cIt.index()];
+ forwardBackwardSolution[i] = sol[i];
+ forwardBackwardSolution[cIt.index()] = sol[cIt.index()];
+ backwardForwardSolution[i] = sol[i];
+ backwardForwardSolution[cIt.index()] = sol[cIt.index()];
+ backwardBackwardSolution[i] = sol[i];
+ backwardBackwardSolution[cIt.index()] = sol[cIt.index()];
+
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+
+}
+
template <class GridType>
void Dune::RodAssembler<GridType>::
strainDerivative(const std::vector<Configuration>& localSolution,
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 106dc3153db86c94d5eaed73bf3276a90ed29e66..fe98940a3b0d2a115299f7cac133ea91806c2fee 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -5,10 +5,118 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
+#include <dune/disc/operators/localstiffness.hh>
#include "../../common/boundarypatch.hh"
#include "configuration.hh"
+template<class GridType, class RT>
+class RodLocalStiffness
+ : public Dune::LocalStiffness<RodLocalStiffness<GridType,RT>,GridType,RT,6>
+{
+ // grid types
+ typedef typename GridType::ctype DT;
+ typedef typename GridType::template Codim<0>::Entity Entity;
+ typedef typename GridType::template Codim<0>::EntityPointer EntityPointer;
+
+ // some other sizes
+ enum {dim=GridType::dimension};
+
+public:
+ // define the number of components of your system, this is used outside
+ // to allocate the correct size of (dense) blocks with a FieldMatrix
+ enum {m=6};
+
+ enum {SIZE = Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE};
+
+ // types for matrics, vectors and boundary conditions
+ typedef Dune::FieldMatrix<RT,m,m> MBlockType; // one entry in the stiffness matrix
+ typedef Dune::FieldVector<RT,m> VBlockType; // one entry in the global vectors
+ typedef Dune::array<Dune::BoundaryConditions::Flags,m> BCBlockType; // componentwise boundary conditions
+
+ // /////////////////////////////////
+ // The material parameters
+ // /////////////////////////////////
+
+ /** \brief Material constants */
+ double K_[3];
+ double A_[3];
+
+ //! Default Constructor
+ RodLocalStiffness ()
+ {
+ // this->currentsize_ = 0;
+
+ // For the time being: all boundary conditions are homogeneous Neumann
+ // This means no boundary condition handling is done at all
+ for (int i=0; i<Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE; i++)
+ for (size_t j=0; j<this->bctype[i].size(); j++)
+ this->bctype[i][j] = Dune::BoundaryConditions::neumann;
+ }
+
+ //! Default Constructor
+ RodLocalStiffness (const Dune::array<double,3>& K, const Dune::array<double,3>& A)
+ {
+ for (int i=0; i<3; i++) {
+ K_[i] = K[i];
+ A_[i] = A[i];
+ }
+ // this->currentsize_ = 0;
+
+ // For the time being: all boundary conditions are homogeneous Neumann
+ // This means no boundary condition handling is done at all
+ for (int i=0; i<Dune::LocalStiffness<RodLocalStiffness,GridType,RT,m>::SIZE; i++)
+ for (size_t j=0; j<this->bctype[i].size(); j++)
+ this->bctype[i][j] = Dune::BoundaryConditions::neumann;
+ }
+
+ //! assemble local stiffness matrix for given element and order
+ /*! On exit the following things have been done:
+ - The stiffness matrix for the given entity and polynomial degree has been assembled and is
+ accessible with the mat() method.
+ - The boundary conditions have been evaluated and are accessible with the bc() method
+ - The right hand side has been assembled. It contains either the value of the essential boundary
+ condition or the assembled source term and neumann boundary condition. It is accessible via the rhs() method.
+ @param[in] e a codim 0 entity reference
+ \param[in] localSolution Current local solution, because this is a nonlinear assembler
+ @param[in] k order of Lagrange basis
+ */
+ template<typename TypeTag>
+ void assemble (const Entity& e,
+ const Dune::BlockVector<Dune::FieldVector<double, dim> >& localSolution,
+ int k=1);
+
+
+ RT energy (const EntityPointer& e,
+ const std::vector<Configuration>& localSolution,
+ const std::vector<Configuration>& localReferenceConfiguration,
+ int k=1);
+
+ Dune::FieldVector<double, 6> getStrain(const std::vector<Configuration>& localSolution,
+ const EntityPointer& element,
+ const Dune::FieldVector<double,1>& pos) const;
+
+ template <class T>
+ static Dune::FieldVector<T,3> darboux(const Quaternion<T>& q, const Dune::FieldVector<T,4>& q_s)
+ {
+ Dune::FieldVector<double,3> u; // The Darboux vector
+
+ u[0] = 2 * (q.B(0) * q_s);
+ u[1] = 2 * (q.B(1) * q_s);
+ u[2] = 2 * (q.B(2) * q_s);
+
+ return u;
+ }
+
+
+ //! should allow to assmble boundary conditions only
+// template<typename Tag>
+// void assembleBoundaryCondition (const Entity& e, int k=1)
+// {
+// }
+
+};
+
namespace Dune
{
@@ -44,6 +152,17 @@ namespace Dune
/** \brief The stress-free configuration */
std::vector<Configuration> referenceConfiguration_;
+ /** \todo Only for the fd approximations */
+ static void infinitesimalVariation(Configuration& c, double eps, int i)
+ {
+ if (i<3)
+ c.r[i] += eps;
+ else
+ c.q = c.q.mult(Quaternion<double>::exp((i==3)*eps,
+ (i==4)*eps,
+ (i==5)*eps));
+ }
+
public:
//! ???
@@ -106,11 +225,16 @@ namespace Dune
//exit(0);
}
- /** \brief Assemble the tangent stiffness matrix and the right hand side
+ /** \brief Assemble the tangent stiffness matrix
*/
void assembleMatrix(const std::vector<Configuration>& sol,
BCRSMatrix<MatrixBlock>& matrix) const;
+ /** \brief Assemble the tangent stiffness matrix using a finite difference approximation
+ */
+ void assembleMatrixFD(const std::vector<Configuration>& sol,
+ BCRSMatrix<MatrixBlock>& matrix) const;
+
void strainDerivative(const std::vector<Configuration>& localSolution,
double pos,
FieldVector<double,1> shapeGrad[2],