From 2415d3434ecaf108099d8f1eb2250acc23fb61e6 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 11 Jul 2011 13:09:58 +0000
Subject: [PATCH] Use Dune::array instead of std::vector to store coefficients
and barycentric coordinates. This is easy because we know the number of
element corners at compile time (I currently don't see support for
non-simplex elements at all). This saves time, but not much: about 1%.
[[Imported from SVN: r7550]]
---
dune/gfe/averagedistanceassembler.hh | 11 ++++++-----
dune/gfe/cosseratenergystiffness.hh | 4 ++--
dune/gfe/geodesicfeassembler.hh | 8 +++-----
dune/gfe/localgeodesicfefunction.hh | 26 ++++++++++++--------------
dune/gfe/localgeodesicfestiffness.hh | 22 +++++++++++-----------
dune/gfe/targetspacertrsolver.cc | 10 +++++-----
dune/gfe/targetspacertrsolver.hh | 11 +++++++----
7 files changed, 46 insertions(+), 46 deletions(-)
diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index a940fe65..e92ac0ce 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -3,7 +3,8 @@
#include <vector>
-template <class TargetSpace>
+/** \tparam N Number of coefficients (i.e., simplex corners) */
+template <class TargetSpace, int N>
class AverageDistanceAssembler
{
static const int size = TargetSpace::TangentVector::size;
@@ -11,8 +12,8 @@ class AverageDistanceAssembler
public:
- AverageDistanceAssembler(const std::vector<TargetSpace>& coefficients,
- const std::vector<double>& weights)
+ AverageDistanceAssembler(const Dune::array<TargetSpace,N>& coefficients,
+ const Dune::array<double,N>& weights)
: coefficients_(coefficients),
weights_(weights)
{}
@@ -74,9 +75,9 @@ public:
}
- const std::vector<TargetSpace> coefficients_;
+ const Dune::array<TargetSpace,N> coefficients_;
- const std::vector<double> weights_;
+ const Dune::array<double,N> weights_;
};
diff --git a/dune/gfe/cosseratenergystiffness.hh b/dune/gfe/cosseratenergystiffness.hh
index b8f53ac6..4dcc4d5f 100644
--- a/dune/gfe/cosseratenergystiffness.hh
+++ b/dune/gfe/cosseratenergystiffness.hh
@@ -95,7 +95,7 @@ public:
/** \brief Assemble the energy for a single element */
RT energy (const Entity& e,
- const std::vector<TargetSpace>& localSolution) const;
+ const Dune::array<TargetSpace, gridDim+1>& localSolution) const;
RT quadraticMembraneEnergy(const Dune::FieldMatrix<double,3,3>& U) const
{
@@ -162,7 +162,7 @@ public:
template <class GridView, int dim>
typename CosseratEnergyLocalStiffness<GridView, dim>::RT CosseratEnergyLocalStiffness<GridView, dim>::
energy(const Entity& element,
- const std::vector<RigidBodyMotion<dim> >& localSolution) const
+ const Dune::array<RigidBodyMotion<dim>, gridDim+1>& localSolution) const
{
RT energy = 0;
diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index 2a13170a..288d07e0 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -120,7 +120,7 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
const int numOfBaseFct = it->template count<gridDim>();
// Extract local solution
- std::vector<TargetSpace> localSolution(numOfBaseFct);
+ Dune::array<TargetSpace,gridDim+1> localSolution;
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -168,7 +168,7 @@ assembleGradient(const std::vector<TargetSpace>& sol,
const int nDofs = it->template count<gridDim>();
// Extract local solution
- std::vector<TargetSpace> localSolution(nDofs);
+ Dune::array<TargetSpace,gridDim+1> localSolution;
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
@@ -198,7 +198,7 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
if (sol.size()!=indexSet.size(gridDim))
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
- std::vector<TargetSpace> localSolution;
+ Dune::array<TargetSpace,gridDim+1> localSolution;
ElementIterator it = gridView_.template begin<0>();
ElementIterator endIt = gridView_.template end<0>();
@@ -206,8 +206,6 @@ computeEnergy(const std::vector<TargetSpace>& sol) const
// Loop over all elements
for (; it!=endIt; ++it) {
- localSolution.resize(it->template count<gridDim>());
-
for (int i=0; i<it->template count<gridDim>(); i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index bf7f0f67..2a6760e2 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -31,11 +31,9 @@ class LocalGeodesicFEFunction
public:
/** \brief Constructor */
- LocalGeodesicFEFunction(const std::vector<TargetSpace>& coefficients)
+ LocalGeodesicFEFunction(const Dune::array<TargetSpace,dim+1>& coefficients)
: coefficients_(coefficients)
- {
- assert(coefficients_.size() == dim+1);
- }
+ {}
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const;
@@ -80,8 +78,8 @@ private:
return B;
}
- static std::vector<ctype> barycentricCoordinates(const Dune::FieldVector<ctype,dim>& local) {
- std::vector<ctype> result(dim+1);
+ static Dune::array<ctype,dim+1> barycentricCoordinates(const Dune::FieldVector<ctype,dim>& local) {
+ Dune::array<ctype,dim+1> result;
result[0] = 1;
for (int i=0; i<dim; i++) {
result[0] -= local[i];
@@ -144,7 +142,7 @@ private:
}
/** \brief The coefficient vector */
- std::vector<TargetSpace> coefficients_;
+ Dune::array<TargetSpace,dim+1> coefficients_;
};
@@ -153,11 +151,11 @@ TargetSpace LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
// First compute the coordinates on the standard simplex (in R^{n+1})
- std::vector<ctype> w = barycentricCoordinates(local);
+ Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
- AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+ AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
- TargetSpaceRiemannianTRSolver<TargetSpace> solver;
+ TargetSpaceRiemannianTRSolver<TargetSpace,dim+1> solver;
solver.setup(&assembler,
coefficients_[0], // initial iterate
@@ -209,8 +207,8 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
Dune::FieldMatrix<ctype,embeddedDim,dim> RHS = dFdw * B;
// the actual system matrix
- std::vector<ctype> w = barycentricCoordinates(local);
- AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+ Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
+ AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleEmbeddedHessian(q,dFdq);
@@ -302,7 +300,7 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
// dFdq
std::vector<ctype> w = barycentricCoordinates(local);
- AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+ AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleEmbeddedHessian(q,dFdq);
@@ -405,7 +403,7 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
// the actual system matrix
std::vector<ctype> w = barycentricCoordinates(local);
- AverageDistanceAssembler<TargetSpace> assembler(coefficients_, w);
+ AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleEmbeddedHessian(q,dFdq);
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 566630e6..27ec09ed 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -40,17 +40,17 @@ public:
*/
virtual void assembleHessian(const Entity& e,
- const std::vector<TargetSpace>& localSolution);
+ const Dune::array<TargetSpace,gridDim+1>& localSolution);
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e,
- const std::vector<TargetSpace>& localSolution) const = 0;
+ const Dune::array<TargetSpace,gridDim+1>& localSolution) const = 0;
/** \brief Assemble the element gradient of the energy functional
The default implementation in this class uses a finite difference approximation */
virtual void assembleGradient(const Entity& element,
- const std::vector<TargetSpace>& solution,
+ const Dune::array<TargetSpace,gridDim+1>& solution,
std::vector<typename TargetSpace::TangentVector>& gradient) const;
// assembled data
@@ -62,7 +62,7 @@ public:
template <class GridView, class TargetSpace>
void LocalGeodesicFEStiffness<GridView, TargetSpace>::
assembleGradient(const Entity& element,
- const std::vector<TargetSpace>& localSolution,
+ const Dune::array<TargetSpace,gridDim+1>& localSolution,
std::vector<typename TargetSpace::TangentVector>& localGradient) const
{
@@ -74,8 +74,8 @@ assembleGradient(const Entity& element,
localGradient.resize(localSolution.size());
- std::vector<TargetSpace> forwardSolution = localSolution;
- std::vector<TargetSpace> backwardSolution = localSolution;
+ Dune::array<TargetSpace,gridDim+1> forwardSolution = localSolution;
+ Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
for (size_t i=0; i<localSolution.size(); i++) {
@@ -111,7 +111,7 @@ assembleGradient(const Entity& element,
template <class GridType, class TargetSpace>
void LocalGeodesicFEStiffness<GridType, TargetSpace>::
assembleHessian(const Entity& element,
- const std::vector<TargetSpace>& localSolution)
+ const Dune::array<TargetSpace,gridDim+1>& localSolution)
{
// 1 degree of freedom per element vertex
size_t nDofs = element.template count<gridDim>();
@@ -144,8 +144,8 @@ assembleHessian(const Entity& element,
Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi;
minusEpsXi *= -1;
- std::vector<TargetSpace> forwardSolution = localSolution;
- std::vector<TargetSpace> backwardSolution = localSolution;
+ Dune::array<TargetSpace,gridDim+1> forwardSolution = localSolution;
+ Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
forwardSolution[i] = TargetSpace::exp(localSolution[i],epsXi);
backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
@@ -169,8 +169,8 @@ assembleHessian(const Entity& element,
Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi; minusEpsXi *= -1;
Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta; minusEpsEta *= -1;
- std::vector<TargetSpace> forwardSolutionXiEta = localSolution;
- std::vector<TargetSpace> backwardSolutionXiEta = localSolution;
+ Dune::array<TargetSpace,gridDim+1> forwardSolutionXiEta = localSolution;
+ Dune::array<TargetSpace,gridDim+1> backwardSolutionXiEta = localSolution;
if (i==j)
forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index e05bcbf5..4cc9caa8 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -6,9 +6,9 @@
#include "maxnormtrustregion.hh"
-template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::
-setup(const AverageDistanceAssembler<TargetSpace>* assembler,
+template <class TargetSpace, int N>
+void TargetSpaceRiemannianTRSolver<TargetSpace,N>::
+setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
const TargetSpace& x,
double tolerance,
int maxTrustRegionSteps,
@@ -45,8 +45,8 @@ setup(const AverageDistanceAssembler<TargetSpace>* assembler,
}
-template <class TargetSpace>
-void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
+template <class TargetSpace, int N>
+void TargetSpaceRiemannianTRSolver<TargetSpace,N>::solve()
{
MaxNormTrustRegion<blocksize> trustRegion(1, // we have only one block
initialTrustRegionRadius_);
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 37997001..55c3e483 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -8,8 +8,11 @@
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
-/** \brief Riemannian trust-region solver for geodesic finite-element problems */
-template <class TargetSpace>
+/** \brief Riemannian trust-region solver for geodesic finite-element problems
+ \tparam TargetSpace The manifold that our functions take values in
+ \tparam N Number of simplex vertices
+ */
+template <class TargetSpace, int N>
class TargetSpaceRiemannianTRSolver
: public NumProc
// : public IterativeSolver<std::vector<TargetSpace>,
@@ -33,7 +36,7 @@ public:
{}
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
- void setup(const AverageDistanceAssembler<TargetSpace>* assembler,
+ void setup(const AverageDistanceAssembler<TargetSpace,N>* assembler,
const TargetSpace& x,
double tolerance,
int maxTrustRegionSteps,
@@ -70,7 +73,7 @@ protected:
double innerTolerance_;
/** \brief The assembler for the average-distance functional */
- const AverageDistanceAssembler<TargetSpace>* assembler_;
+ const AverageDistanceAssembler<TargetSpace,N>* assembler_;
/** \brief The solver for the quadratic inner problems */
std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;
--
GitLab