diff --git a/rod3d.cc b/rod3d.cc
index d58b0fa2c5329c1c2cadba056d665613ab580716..e6921ee7b89d17d2c158cc9339fde2227728fd3e 100644
--- a/rod3d.cc
+++ b/rod3d.cc
@@ -168,7 +168,7 @@ int main (int argc, char *argv[]) try
MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessian);
- rodAssembler.assembleMatrixFD(exactSolution, hessian);
+ rodAssembler.assembleMatrix(exactSolution, hessian);
double error = std::numeric_limits<double>::max();
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index ff1f4ddeece48214873650243da4200465a78087..ca3c0cc20207d1dd5a3e60a9070c4f58b92da693 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -51,11 +51,17 @@ void RodAssembler<GridType>::
assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const
{
- using namespace Dune;
+ // ////////////////////////////////////////////////////
+ // 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<typename GridType::LeafGridView,double> localStiffness(K, A);
const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
- MatrixIndexSet neighborsPerVertex;
+ Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex);
matrix = 0;
@@ -63,35 +69,35 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
- Matrix<MatrixBlock> mat;
-
for( ; it != endit; ++it ) {
- const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
+ const Dune::LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
= Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->type(), elementOrder);
const int numOfBaseFct = baseSet.size();
- mat.setSize(numOfBaseFct, numOfBaseFct);
-
// Extract local solution
std::vector<RigidBodyMotion<3> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
+ localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+
+ localStiffness.localReferenceConfiguration_.resize(numOfBaseFct);
+
+ for (int i=0; i<numOfBaseFct; i++)
+ localStiffness.localReferenceConfiguration_[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
// setup matrix
- //getLocalMatrix( it, localSolution, sol, numOfBaseFct, mat);
- DUNE_THROW(NotImplemented, "getLocalMatrix");
+ localStiffness.assemble(*it, localSolution);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
- int row = indexSet.template subIndex<gridDim>(*it,i);
+ int row = indexSet.subIndex(*it,i,gridDim);
for (int j=0; j<numOfBaseFct; j++ ) {
- int col = indexSet.template subIndex<gridDim>(*it,j);
- matrix[row][col] += mat[i][j];
+ int col = indexSet.subIndex(*it,j,gridDim);
+ matrix[row][col] += localStiffness.mat(i,j);
}
}
@@ -100,257 +106,6 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
}
-template <class GridType>
-void RodAssembler<GridType>::
-assembleMatrixFD(const std::vector<RigidBodyMotion<3> >& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix) const
-{
- using namespace Dune;
-
- double eps = 1e-4;
-
- typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >::row_type::iterator ColumnIterator;
-
- // ///////////////////////////////////////////////////////////
- // Compute gradient by finite-difference approximation
- // ///////////////////////////////////////////////////////////
- std::vector<RigidBodyMotion<3> > forwardSolution = sol;
- std::vector<RigidBodyMotion<3> > backwardSolution = sol;
-
- std::vector<RigidBodyMotion<3> > forwardForwardSolution = sol;
- std::vector<RigidBodyMotion<3> > forwardBackwardSolution = sol;
- std::vector<RigidBodyMotion<3> > backwardForwardSolution = sol;
- std::vector<RigidBodyMotion<3> > 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<typename GridType::LeafGridView,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;
-
- Dune::array<RigidBodyMotion<3>,2> localReferenceConfiguration;
- Dune::array<RigidBodyMotion<3>,2> localSolution;
-
- // ///////////////////////////////////////////////////////////////
- // 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 only the upper right triangular matrix
- if (cIt.index() < i)
- continue;
-
- // ////////////////////////////////////////////////////////////////////////////
- // Compute a finite-difference approximation of this hessian matrix block
- // ////////////////////////////////////////////////////////////////////////////
-
- for (int j=0; j<6; j++) {
-
- for (int k=0; k<6; k++) {
-
- // compute only the upper right triangular matrix
- if (i==cIt.index() && k<j)
- continue;
-
- 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()];
-
- }
-
- }
-
- }
-
- }
-
- }
-
- // ///////////////////////////////////////////////////////////////
- // Symmetrize the matrix
- // This is possible expensive, but I want to be absolute sure
- // that the matrix is symmetric.
- // ///////////////////////////////////////////////////////////////
- for (int i=0; i<matrix.N(); i++) {
-
- ColumnIterator cIt = matrix[i].begin();
- ColumnIterator cEndIt = matrix[i].end();
-
- for (; cIt!=cEndIt; ++cIt) {
-
- if (cIt.index()>i)
- continue;
-
-
- if (cIt.index()==i) {
-
- for (int j=1; j<6; j++)
- for (int k=0; k<j; k++)
- (*cIt)[j][k] = (*cIt)[k][j];
-
- } else {
-
- const FieldMatrix<double,6,6>& other = matrix[cIt.index()][i];
-
- for (int j=0; j<6; j++)
- for (int k=0; k<6; k++)
- (*cIt)[j][k] = other[k][j];
-
-
- }
-
-
- }
-
- }
-
-}
-
-
template <class GridType>
void RodAssembler<GridType>::
assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
@@ -385,13 +140,13 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
const int nDofs = 2;
// Extract local solution
- array<RigidBodyMotion<3>,nDofs> localSolution;
+ std::vector<RigidBodyMotion<3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
// Extract local reference configuration
- array<RigidBodyMotion<3>,nDofs> localReferenceConfiguration;
+ std::vector<RigidBodyMotion<3> > localReferenceConfiguration(nDofs);
for (int i=0; i<nDofs; i++)
localReferenceConfiguration[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
@@ -431,8 +186,8 @@ computeEnergy(const std::vector<RigidBodyMotion<3> >& sol) const
Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
- Dune::array<RigidBodyMotion<3>,2> localReferenceConfiguration;
- Dune::array<RigidBodyMotion<3>,2> localSolution;
+ std::vector<RigidBodyMotion<3> > localReferenceConfiguration(2);
+ std::vector<RigidBodyMotion<3> > localSolution(2);
ElementLeafIterator it = grid_->template leafbegin<0>();
ElementLeafIterator endIt = grid_->template leafend<0>();
@@ -493,7 +248,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
= Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->type(), elementOrder);
int numOfBaseFct = baseSet.size();
- array<RigidBodyMotion<3>,2> localSolution;
+ std::vector<RigidBodyMotion<3> > localSolution(2);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index a466fc18a4a5a92da638005b52b838268351a9ba..ca8346fed099ec7e9aaa307ccbca6781ac09eaa5 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -43,17 +43,6 @@
/** \brief The stress-free configuration */
std::vector<RigidBodyMotion<3> > referenceConfiguration_;
- /** \todo Only for the fd approximations */
- static void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
- {
- if (i<3)
- c.r[i] += eps;
- else
- c.q = c.q.mult(Rotation<3,double>::exp((i==3)*eps,
- (i==4)*eps,
- (i==5)*eps));
- }
-
public:
//! ???
@@ -121,11 +110,6 @@
void assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const;
- /** \brief Assemble the tangent stiffness matrix using a finite difference approximation
- */
- void assembleMatrixFD(const std::vector<RigidBodyMotion<3> >& sol,
- Dune::BCRSMatrix<MatrixBlock>& matrix) const;
-
void assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
diff --git a/src/rodlocalstiffness.hh b/src/rodlocalstiffness.hh
index 6c7c7e592183804cf1ea78016d02f0ec31dcceca..d44826f950fe956f8e609d662685125e650cb417 100644
--- a/src/rodlocalstiffness.hh
+++ b/src/rodlocalstiffness.hh
@@ -6,6 +6,7 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh>
#include <dune/disc/operators/localstiffness.hh>
+#include<dune/disc/operators/boundaryconditions.hh>
#include "rigidbodymotion.hh"
@@ -13,6 +14,8 @@ template<class GridView, class RT>
class RodLocalStiffness
: public Dune::LocalStiffness<GridView,RT,6>
{
+ typedef RigidBodyMotion<3> TargetSpace;
+
// grid types
typedef typename GridView::Grid::ctype DT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -27,6 +30,22 @@ class RodLocalStiffness
// Quadrature order used for the bending and torsion energy
enum {bendingQuadOrder = 2};
+public:
+ /** \brief For the fd approximations
+ \todo This is public because RodAssembler uses it
+ */
+ static void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
+ {
+ if (i<3)
+ c.r[i] += eps;
+ else
+ c.q = c.q.mult(Rotation<3,double>::exp((i==3)*eps,
+ (i==4)*eps,
+ (i==5)*eps));
+ }
+
+ std::vector<RigidBodyMotion<3> > localReferenceConfiguration_;
+
public:
//! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
@@ -61,17 +80,11 @@ public:
A_[i] = A[i];
}
}
+
+ void assemble(const Entity& e,
+ const std::vector<TargetSpace>& localSolution);
- //! 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
+ /** \brief assemble local stiffness matrix for given element and order
*/
void assemble (const Entity& e,
const Dune::BlockVector<Dune::FieldVector<double, 6> >& localSolution,
@@ -93,8 +106,8 @@ public:
RT energy (const Entity& e,
- const Dune::array<RigidBodyMotion<3>,2>& localSolution,
- const Dune::array<RigidBodyMotion<3>,2>& localReferenceConfiguration,
+ const std::vector<RigidBodyMotion<3> >& localSolution,
+ const std::vector<RigidBodyMotion<3> >& localReferenceConfiguration,
int k=1);
static void interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
@@ -103,14 +116,14 @@ public:
static void interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad);
- Dune::FieldVector<double, 6> getStrain(const Dune::array<RigidBodyMotion<3>,2>& localSolution,
+ Dune::FieldVector<double, 6> getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
const Entity& element,
const Dune::FieldVector<double,1>& pos) const;
/** \brief Assemble the element gradient of the energy functional */
void assembleGradient(const Entity& element,
- const Dune::array<RigidBodyMotion<3>,2>& solution,
- const Dune::array<RigidBodyMotion<3>,2>& referenceConfiguration,
+ const std::vector<RigidBodyMotion<3> >& solution,
+ const std::vector<RigidBodyMotion<3> >& referenceConfiguration,
Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const;
template <class T>
@@ -130,8 +143,8 @@ public:
template <class GridType, class RT>
RT RodLocalStiffness<GridType, RT>::
energy(const Entity& element,
- const Dune::array<RigidBodyMotion<3>,2>& localSolution,
- const Dune::array<RigidBodyMotion<3>,2>& localReferenceConfiguration,
+ const std::vector<RigidBodyMotion<3> >& localSolution,
+ const std::vector<RigidBodyMotion<3> >& localReferenceConfiguration,
int k)
{
RT energy = 0;
@@ -403,7 +416,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
template <class GridType, class RT>
Dune::FieldVector<double, 6> RodLocalStiffness<GridType, RT>::
-getStrain(const Dune::array<RigidBodyMotion<3>,2>& localSolution,
+getStrain(const std::vector<RigidBodyMotion<3> >& localSolution,
const Entity& element,
const Dune::FieldVector<double,1>& pos) const
{
@@ -479,8 +492,8 @@ getStrain(const Dune::array<RigidBodyMotion<3>,2>& localSolution,
template <class GridType, class RT>
void RodLocalStiffness<GridType, RT>::
assembleGradient(const Entity& element,
- const Dune::array<RigidBodyMotion<3>,2>& solution,
- const Dune::array<RigidBodyMotion<3>,2>& referenceConfiguration,
+ const std::vector<RigidBodyMotion<3> >& solution,
+ const std::vector<RigidBodyMotion<3> >& referenceConfiguration,
Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const
{
using namespace Dune;
@@ -666,5 +679,167 @@ assembleGradient(const Entity& element,
}
}
+
+template <class GridType, class RT>
+void RodLocalStiffness<GridType,RT>::
+assemble(const Entity& element,
+ const std::vector<TargetSpace>& localSolution)
+{
+ // 1 degree of freedom per element vertex
+ int nDofs = element.template count<dim>();
+
+ // Clear assemble data
+ this->setcurrentsize(nDofs);
+
+ this->A = 0;
+
+ for (int i=0; i<nDofs; i++) {
+ this->b[i] = 0;
+ for (int j=0; j<this->bctype[i].size(); j++)
+ this->bctype[i][j] = Dune::BoundaryConditions::neumann;
+ }
+
+ double eps = 1e-4;
+
+ typedef typename Dune::Matrix<Dune::FieldMatrix<double,6,6> >::row_type::iterator ColumnIterator;
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<RigidBodyMotion<3> > forwardSolution = localSolution;
+ std::vector<RigidBodyMotion<3> > backwardSolution = localSolution;
+
+ std::vector<RigidBodyMotion<3> > forwardForwardSolution = localSolution;
+ std::vector<RigidBodyMotion<3> > forwardBackwardSolution = localSolution;
+ std::vector<RigidBodyMotion<3> > backwardForwardSolution = localSolution;
+ std::vector<RigidBodyMotion<3> > backwardBackwardSolution = localSolution;
+
+ // ///////////////////////////////////////////////////////////////
+ // Loop over all blocks of the element matrix
+ // ///////////////////////////////////////////////////////////////
+ for (int i=0; i<this->A.N(); i++) {
+
+ ColumnIterator cIt = this->A[i].begin();
+ ColumnIterator cEndIt = this->A[i].end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // compute only the upper right triangular matrix
+ if (cIt.index() < i)
+ continue;
+
+ // ////////////////////////////////////////////////////////////////////////////
+ // Compute a finite-difference approximation of this hessian matrix block
+ // ////////////////////////////////////////////////////////////////////////////
+
+ for (int j=0; j<6; j++) {
+
+ for (int k=0; k<6; k++) {
+
+ // compute only the upper right triangular matrix
+ if (i==cIt.index() && k<j)
+ continue;
+
+ // Diagonal entries
+ if (i==cIt.index() && j==k) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+ double forwardEnergy = energy(element, forwardSolution, localReferenceConfiguration_);
+
+ double solutionEnergy = energy(element, localSolution, localReferenceConfiguration_);
+
+ double backwardEnergy = energy(element, backwardSolution, localReferenceConfiguration_);
+
+ // Second derivative
+ (*cIt)[j][k] = (forwardEnergy - 2*solutionEnergy + backwardEnergy) / (eps*eps);
+
+ forwardSolution[i] = localSolution[i];
+ backwardSolution[i] = localSolution[i];
+
+ } else {
+
+ // Off-diagonal entries
+ 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);
+
+ double forwardForwardEnergy = energy(element, forwardForwardSolution, localReferenceConfiguration_);
+
+ double forwardBackwardEnergy = energy(element, forwardBackwardSolution, localReferenceConfiguration_);
+
+ double backwardForwardEnergy = energy(element, backwardForwardSolution, localReferenceConfiguration_);
+
+ double backwardBackwardEnergy = energy(element, backwardBackwardSolution, localReferenceConfiguration_);
+
+ (*cIt)[j][k] = (forwardForwardEnergy + backwardBackwardEnergy
+ - forwardBackwardEnergy - backwardForwardEnergy) / (4*eps*eps);
+
+ forwardForwardSolution[i] = localSolution[i];
+ forwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
+ forwardBackwardSolution[i] = localSolution[i];
+ forwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
+ backwardForwardSolution[i] = localSolution[i];
+ backwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
+ backwardBackwardSolution[i] = localSolution[i];
+ backwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
+
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+
+ // ///////////////////////////////////////////////////////////////
+ // Symmetrize the matrix
+ // This is possible expensive, but I want to be absolute sure
+ // that the matrix is symmetric.
+ // ///////////////////////////////////////////////////////////////
+ for (int i=0; i<this->A.N(); i++) {
+
+ ColumnIterator cIt = this->A[i].begin();
+ ColumnIterator cEndIt = this->A[i].end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ if (cIt.index()>i)
+ continue;
+
+
+ if (cIt.index()==i) {
+
+ for (int j=1; j<6; j++)
+ for (int k=0; k<j; k++)
+ (*cIt)[j][k] = (*cIt)[k][j];
+
+ } else {
+
+ const Dune::FieldMatrix<double,6,6>& other = this->A[cIt.index()][i];
+
+ for (int j=0; j<6; j++)
+ for (int k=0; k<6; k++)
+ (*cIt)[j][k] = other[k][j];
+
+
+ }
+
+
+ }
+
+ }
+
+}
+
+
#endif
diff --git a/src/rodsolver.cc b/src/rodsolver.cc
index cd94020613badcc93883471743d7a7af00c13c9b..f048fedb9affd9e55bb121db50ef4bd6813b5b50 100644
--- a/src/rodsolver.cc
+++ b/src/rodsolver.cc
@@ -189,8 +189,8 @@ void RodSolver<GridType>::solve()
corr = 0;
rodAssembler_->assembleGradient(x_, rhs);
- //rodAssembler_->assembleMatrix(x_, *hessianMatrix_);
- rodAssembler_->assembleMatrixFD(x_, *hessianMatrix_);
+ rodAssembler_->assembleMatrix(x_, *hessianMatrix_);
+ //rodAssembler_->assembleMatrixFD(x_, *hessianMatrix_);
//gradientFDCheck(x_, rhs, *rodAssembler_);
//hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_);