From 5b5aeb8e93449097823662281aa1d6242d841f53 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 15 Aug 2007 09:46:13 +0000
Subject: [PATCH] a method to compute the gradient by fd approximation
[[Imported from SVN: r1524]]
---
src/rodassembler.cc | 100 ++++++++++++++++++++++++++++++++++++++++++++
src/rodassembler.hh | 10 ++---
2 files changed, 103 insertions(+), 7 deletions(-)
diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index bb09d1e1..93a525cc 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -1183,6 +1183,106 @@ assembleGradient(const std::vector<Configuration>& sol,
}
+template <class GridType>
+void RodAssembler<GridType>::
+assembleGradientFD(const std::vector<Configuration>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
+{
+ double eps = 1e-2;
+ if (grad.size() != sol.size())
+ grad.resize(sol.size());
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = sol;
+ std::vector<Configuration> backwardSolution = 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 (size_t i=0; i<sol.size(); i++) {
+
+ for (int j=0; j<6; j++) {
+
+ double forwardEnergy = 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] = 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] = backwardSolution[i-1];
+ localSolution[1] = backwardSolution[i];
+
+ backwardEnergy += localStiffness.energy(elements[i-1], localSolution, localReferenceConfiguration);
+
+ }
+
+ // Second derivative
+ grad[i][j] = (forwardEnergy - backwardEnergy) / (2*eps);
+
+ forwardSolution[i] = sol[i];
+ backwardSolution[i] = sol[i];
+
+ }
+
+ }
+
+
+}
+
template <class GridType>
double RodAssembler<GridType>::
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index d35dbc51..eaecc939 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -108,13 +108,6 @@ public:
return u;
}
-
- //! should allow to assmble boundary conditions only
-// template<typename Tag>
-// void assembleBoundaryCondition (const Entity& e, int k=1)
-// {
-// }
-
};
@@ -253,6 +246,9 @@ public:
void assembleGradient(const std::vector<Configuration>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+ void assembleGradientFD(const std::vector<Configuration>& sol,
+ Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const;
+
/** \brief Compute the energy of a deformation state */
double computeEnergy(const std::vector<Configuration>& sol) const;
--
GitLab