Skip to content
Snippets Groups Projects
Commit 94ecd801 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

The code to assemble the stiffness matrices by finite-difference approximations 

now works on a per-element basis and has been moved to the local rod assembler.
That way, it has lost most of its 1d-ness.  In the future, it will be usable
for other problems as well.

[[Imported from SVN: r4019]]
parent a4dbdf25
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
rod3d.cc
+ 1
− 1
View file @ 94ecd801
...@@ -168,7 +168,7 @@ int main (int argc, char *argv[]) try ...@@ -168,7 +168,7 @@ int main (int argc, char *argv[]) try
MatrixIndexSet indices(exactSolution.size(), exactSolution.size()); MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
rodAssembler.getNeighborsPerVertex(indices); rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessian); indices.exportIdx(hessian);
rodAssembler.assembleMatrixFD(exactSolution, hessian); rodAssembler.assembleMatrix(exactSolution, hessian);
double error = std::numeric_limits<double>::max(); double error = std::numeric_limits<double>::max();
......
src/rodassembler.cc
+ 24
− 269
View file @ 94ecd801
...@@ -51,11 +51,17 @@ void RodAssembler<GridType>:: ...@@ -51,11 +51,17 @@ void RodAssembler<GridType>::
assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol, assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const 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()); const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel());
MatrixIndexSet neighborsPerVertex; Dune::MatrixIndexSet neighborsPerVertex;
getNeighborsPerVertex(neighborsPerVertex); getNeighborsPerVertex(neighborsPerVertex);
matrix = 0; matrix = 0;
...@@ -63,35 +69,35 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol, ...@@ -63,35 +69,35 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() ); ElementIterator it = grid_->template lbegin<0>( grid_->maxLevel() );
ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() ); ElementIterator endit = grid_->template lend<0> ( grid_->maxLevel() );
Matrix<MatrixBlock> mat;
for( ; it != endit; ++it ) { 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); = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->type(), elementOrder);
const int numOfBaseFct = baseSet.size(); const int numOfBaseFct = baseSet.size();
mat.setSize(numOfBaseFct, numOfBaseFct);
// Extract local solution // Extract local solution
std::vector<RigidBodyMotion<3> > localSolution(numOfBaseFct); std::vector<RigidBodyMotion<3> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++) 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 // setup matrix
//getLocalMatrix( it, localSolution, sol, numOfBaseFct, mat); localStiffness.assemble(*it, localSolution);
DUNE_THROW(NotImplemented, "getLocalMatrix");
// Add element matrix to global stiffness matrix // Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) { 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++ ) { for (int j=0; j<numOfBaseFct; j++ ) {
int col = indexSet.template subIndex<gridDim>(*it,j); int col = indexSet.subIndex(*it,j,gridDim);
matrix[row][col] += mat[i][j]; matrix[row][col] += localStiffness.mat(i,j);
} }
} }
...@@ -100,257 +106,6 @@ assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol, ...@@ -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> template <class GridType>
void RodAssembler<GridType>:: void RodAssembler<GridType>::
assembleGradient(const std::vector<RigidBodyMotion<3> >& sol, assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
...@@ -385,13 +140,13 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol, ...@@ -385,13 +140,13 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
const int nDofs = 2; const int nDofs = 2;
// Extract local solution // Extract local solution
array<RigidBodyMotion<3>,nDofs> localSolution; std::vector<RigidBodyMotion<3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++) for (int i=0; i<nDofs; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)]; localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
// Extract local reference configuration // Extract local reference configuration
array<RigidBodyMotion<3>,nDofs> localReferenceConfiguration; std::vector<RigidBodyMotion<3> > localReferenceConfiguration(nDofs);
for (int i=0; i<nDofs; i++) for (int i=0; i<nDofs; i++)
localReferenceConfiguration[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)]; localReferenceConfiguration[i] = referenceConfiguration_[indexSet.subIndex(*it,i,gridDim)];
...@@ -431,8 +186,8 @@ computeEnergy(const std::vector<RigidBodyMotion<3> >& sol) const ...@@ -431,8 +186,8 @@ computeEnergy(const std::vector<RigidBodyMotion<3> >& sol) const
Dune::array<double,3> A = {A_[0], A_[1], A_[2]}; Dune::array<double,3> A = {A_[0], A_[1], A_[2]};
RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A); RodLocalStiffness<typename GridType::LeafGridView,double> localStiffness(K, A);
Dune::array<RigidBodyMotion<3>,2> localReferenceConfiguration; std::vector<RigidBodyMotion<3> > localReferenceConfiguration(2);
Dune::array<RigidBodyMotion<3>,2> localSolution; std::vector<RigidBodyMotion<3> > localSolution(2);
ElementLeafIterator it = grid_->template leafbegin<0>(); ElementLeafIterator it = grid_->template leafbegin<0>();
ElementLeafIterator endIt = grid_->template leafend<0>(); ElementLeafIterator endIt = grid_->template leafend<0>();
...@@ -493,7 +248,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol, ...@@ -493,7 +248,7 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
= Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->type(), elementOrder); = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(it->type(), elementOrder);
int numOfBaseFct = baseSet.size(); int numOfBaseFct = baseSet.size();
array<RigidBodyMotion<3>,2> localSolution; std::vector<RigidBodyMotion<3> > localSolution(2);
for (int i=0; i<numOfBaseFct; i++) for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)]; localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
......
src/rodassembler.hh
+ 0
− 16
View file @ 94ecd801
...@@ -43,17 +43,6 @@ ...@@ -43,17 +43,6 @@
/** \brief The stress-free configuration */ /** \brief The stress-free configuration */
std::vector<RigidBodyMotion<3> > referenceConfiguration_; 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: public:
//! ??? //! ???
...@@ -121,11 +110,6 @@ ...@@ -121,11 +110,6 @@
void assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol, void assembleMatrix(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix) const; 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, void assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const; Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
......
src/rodlocalstiffness.hh
+ 195
− 20
View file @ 94ecd801
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
#include <dune/istl/matrixindexset.hh> #include <dune/istl/matrixindexset.hh>
#include <dune/istl/matrix.hh> #include <dune/istl/matrix.hh>
#include <dune/disc/operators/localstiffness.hh> #include <dune/disc/operators/localstiffness.hh>
#include<dune/disc/operators/boundaryconditions.hh>
#include "rigidbodymotion.hh" #include "rigidbodymotion.hh"
...@@ -13,6 +14,8 @@ template<class GridView, class RT> ...@@ -13,6 +14,8 @@ template<class GridView, class RT>
class RodLocalStiffness class RodLocalStiffness
: public Dune::LocalStiffness<GridView,RT,6> : public Dune::LocalStiffness<GridView,RT,6>
{ {
typedef RigidBodyMotion<3> TargetSpace;
// grid types // grid types
typedef typename GridView::Grid::ctype DT; typedef typename GridView::Grid::ctype DT;
typedef typename GridView::template Codim<0>::Entity Entity; typedef typename GridView::template Codim<0>::Entity Entity;
...@@ -27,6 +30,22 @@ class RodLocalStiffness ...@@ -27,6 +30,22 @@ class RodLocalStiffness
// Quadrature order used for the bending and torsion energy // Quadrature order used for the bending and torsion energy
enum {bendingQuadOrder = 2}; 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: public:
//! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d //! Each block is x, y, theta in 2d, T (R^3 \times SO(3)) in 3d
...@@ -61,17 +80,11 @@ public: ...@@ -61,17 +80,11 @@ public:
A_[i] = A[i]; A_[i] = A[i];
} }
} }
void assemble(const Entity& e,
const std::vector<TargetSpace>& localSolution);
//! assemble local stiffness matrix for given element and order /** \brief 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
*/ */
void assemble (const Entity& e, void assemble (const Entity& e,
const Dune::BlockVector<Dune::FieldVector<double, 6> >& localSolution, const Dune::BlockVector<Dune::FieldVector<double, 6> >& localSolution,
...@@ -93,8 +106,8 @@ public: ...@@ -93,8 +106,8 @@ public:
RT energy (const Entity& e, RT energy (const Entity& e,
const Dune::array<RigidBodyMotion<3>,2>& localSolution, const std::vector<RigidBodyMotion<3> >& localSolution,
const Dune::array<RigidBodyMotion<3>,2>& localReferenceConfiguration, const std::vector<RigidBodyMotion<3> >& localReferenceConfiguration,
int k=1); int k=1);
static void interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s, static void interpolationDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
...@@ -103,14 +116,14 @@ public: ...@@ -103,14 +116,14 @@ public:
static void interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s, static void interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& q1, double s,
double intervalLength, Dune::array<Quaternion<double>,6>& grad); 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 Entity& element,
const Dune::FieldVector<double,1>& pos) const; const Dune::FieldVector<double,1>& pos) const;
/** \brief Assemble the element gradient of the energy functional */ /** \brief Assemble the element gradient of the energy functional */
void assembleGradient(const Entity& element, void assembleGradient(const Entity& element,
const Dune::array<RigidBodyMotion<3>,2>& solution, const std::vector<RigidBodyMotion<3> >& solution,
const Dune::array<RigidBodyMotion<3>,2>& referenceConfiguration, const std::vector<RigidBodyMotion<3> >& referenceConfiguration,
Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const; Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const;
template <class T> template <class T>
...@@ -130,8 +143,8 @@ public: ...@@ -130,8 +143,8 @@ public:
template <class GridType, class RT> template <class GridType, class RT>
RT RodLocalStiffness<GridType, RT>:: RT RodLocalStiffness<GridType, RT>::
energy(const Entity& element, energy(const Entity& element,
const Dune::array<RigidBodyMotion<3>,2>& localSolution, const std::vector<RigidBodyMotion<3> >& localSolution,
const Dune::array<RigidBodyMotion<3>,2>& localReferenceConfiguration, const std::vector<RigidBodyMotion<3> >& localReferenceConfiguration,
int k) int k)
{ {
RT energy = 0; RT energy = 0;
...@@ -403,7 +416,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>& ...@@ -403,7 +416,7 @@ interpolationVelocityDerivative(const Rotation<3,RT>& q0, const Rotation<3,RT>&
template <class GridType, class RT> template <class GridType, class RT>
Dune::FieldVector<double, 6> RodLocalStiffness<GridType, 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 Entity& element,
const Dune::FieldVector<double,1>& pos) const const Dune::FieldVector<double,1>& pos) const
{ {
...@@ -479,8 +492,8 @@ getStrain(const Dune::array<RigidBodyMotion<3>,2>& localSolution, ...@@ -479,8 +492,8 @@ getStrain(const Dune::array<RigidBodyMotion<3>,2>& localSolution,
template <class GridType, class RT> template <class GridType, class RT>
void RodLocalStiffness<GridType, RT>:: void RodLocalStiffness<GridType, RT>::
assembleGradient(const Entity& element, assembleGradient(const Entity& element,
const Dune::array<RigidBodyMotion<3>,2>& solution, const std::vector<RigidBodyMotion<3> >& solution,
const Dune::array<RigidBodyMotion<3>,2>& referenceConfiguration, const std::vector<RigidBodyMotion<3> >& referenceConfiguration,
Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const
{ {
using namespace Dune; using namespace Dune;
...@@ -666,5 +679,167 @@ assembleGradient(const Entity& element, ...@@ -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 #endif
src/rodsolver.cc
+ 2
− 2
View file @ 94ecd801
...@@ -189,8 +189,8 @@ void RodSolver<GridType>::solve() ...@@ -189,8 +189,8 @@ void RodSolver<GridType>::solve()
corr = 0; corr = 0;
rodAssembler_->assembleGradient(x_, rhs); rodAssembler_->assembleGradient(x_, rhs);
//rodAssembler_->assembleMatrix(x_, *hessianMatrix_); rodAssembler_->assembleMatrix(x_, *hessianMatrix_);
rodAssembler_->assembleMatrixFD(x_, *hessianMatrix_); //rodAssembler_->assembleMatrixFD(x_, *hessianMatrix_);
//gradientFDCheck(x_, rhs, *rodAssembler_); //gradientFDCheck(x_, rhs, *rodAssembler_);
//hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_); //hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment