diff --git a/dune/gfe/geodesicfeassembler.hh b/dune/gfe/geodesicfeassembler.hh
index cf217226d773b3f123a2e7981120be48a685daf7..327176c5762917c1ac13f86a5de7a5d2c11c3c60 100644
--- a/dune/gfe/geodesicfeassembler.hh
+++ b/dune/gfe/geodesicfeassembler.hh
@@ -11,13 +11,12 @@
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
*/
-template <class GridView, class TargetSpace>
+template <class Basis, class TargetSpace>
class GeodesicFEAssembler {
-
- typedef typename GridView::template Codim<0>::Entity EntityType;
- typedef typename GridView::template Codim<0>::EntityPointer EntityPointer;
+
+ typedef typename Basis::GridView GridView;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
+
//! Dimension of the grid.
enum { gridDim = GridView::dimension };
@@ -26,19 +25,21 @@ class GeodesicFEAssembler {
//!
typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
+
protected:
- const GridView gridView_;
+ const Basis basis_;
- LocalGeodesicFEStiffness<GridView,TargetSpace>* localStiffness_;
+ LocalGeodesicFEStiffness<GridView,
+ typename Basis::LocalFiniteElement,
+ TargetSpace>* localStiffness_;
public:
/** \brief Constructor for a given grid */
- GeodesicFEAssembler(const GridView& gridView,
- LocalGeodesicFEStiffness<GridView,TargetSpace>* localStiffness)
- : gridView_(gridView),
+ GeodesicFEAssembler(const Basis& basis,
+ LocalGeodesicFEStiffness<GridView,typename Basis::LocalFiniteElement, TargetSpace>* localStiffness)
+ : basis_(basis),
localStiffness_(localStiffness)
{}
@@ -62,27 +63,27 @@ public:
-template <class GridView, class TargetSpace>
-void GeodesicFEAssembler<GridView,TargetSpace>::
+template <class Basis, class TargetSpace>
+void GeodesicFEAssembler<Basis,TargetSpace>::
getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
{
- const typename GridView::IndexSet& indexSet = gridView_.indexSet();
-
- int n = gridView_.size(gridDim);
+ int n = basis_.size();
nb.resize(n, n);
- ElementIterator it = gridView_.template begin<0>();
- ElementIterator endit = gridView_.template end<0> ();
+ ElementIterator it = basis_.getGridView().template begin<0>();
+ ElementIterator endit = basis_.getGridView().template end<0> ();
for (; it!=endit; ++it) {
+
+ const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(*it);
- for (int i=0; i<it->template count<gridDim>(); i++) {
+ for (int i=0; i<lfe.localBasis().size(); i++) {
- for (int j=0; j<it->template count<gridDim>(); j++) {
+ for (int j=0; j<lfe.localBasis().size(); j++) {
- int iIdx = indexSet.subIndex(*it,i,gridDim);
- int jIdx = indexSet.subIndex(*it,j,gridDim);
+ int iIdx = basis_.index(*it,i);
+ int jIdx = basis_.index(*it,j);
nb.add(iIdx, jIdx);
@@ -94,14 +95,12 @@ getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
}
-template <class GridView, class TargetSpace>
-void GeodesicFEAssembler<GridView,TargetSpace>::
+template <class Basis, class TargetSpace>
+void GeodesicFEAssembler<Basis,TargetSpace>::
assembleMatrix(const std::vector<TargetSpace>& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix,
bool computeOccupationPattern) const
{
- const typename GridView::IndexSet& indexSet = gridView_.indexSet();
-
if (computeOccupationPattern) {
Dune::MatrixIndexSet neighborsPerVertex;
@@ -112,30 +111,30 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
matrix = 0;
- ElementIterator it = gridView_.template begin<0>();
- ElementIterator endit = gridView_.template end<0> ();
+ ElementIterator it = basis_.getGridView().template begin<0>();
+ ElementIterator endit = basis_.getGridView().template end<0> ();
for( ; it != endit; ++it ) {
- const int numOfBaseFct = it->template count<gridDim>();
+ const int numOfBaseFct = basis_.getLocalFiniteElement(*it).localBasis().size();
// Extract local solution
- Dune::array<TargetSpace,gridDim+1> localSolution;
+ std::vector<TargetSpace> localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+ localSolution[i] = sol[basis_.index(*it,i)];
// setup matrix
- localStiffness_->assembleHessian(*it, localSolution);
+ localStiffness_->assembleHessian(*it, basis_.getLocalFiniteElement(*it), localSolution);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
- int row = indexSet.subIndex(*it,i,gridDim);
+ int row = basis_.index(*it,i);
for (int j=0; j<numOfBaseFct; j++ ) {
- int col = indexSet.subIndex(*it,j,gridDim);
+ int col = basis_.index(*it,j);
matrix[row][col] += localStiffness_->A_[i][j];
}
@@ -145,71 +144,70 @@ assembleMatrix(const std::vector<TargetSpace>& sol,
}
-template <class GridView, class TargetSpace>
-void GeodesicFEAssembler<GridView,TargetSpace>::
+template <class Basis, class TargetSpace>
+void GeodesicFEAssembler<Basis,TargetSpace>::
assembleGradient(const std::vector<TargetSpace>& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
- const typename GridView::IndexSet& indexSet = gridView_.indexSet();
-
- if (sol.size()!=gridView_.size(gridDim))
+ if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
- ElementIterator it = gridView_.template begin<0>();
- ElementIterator endIt = gridView_.template end<0>();
+ ElementIterator it = basis_.getGridView().template begin<0>();
+ ElementIterator endIt = basis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
// A 1d grid has two vertices
- const int nDofs = it->template count<gridDim>();
+ const int nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
// Extract local solution
- Dune::array<TargetSpace,gridDim+1> localSolution;
+ std::vector<TargetSpace> localSolution(nDofs);
for (int i=0; i<nDofs; i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+ localSolution[i] = sol[basis_.index(*it,i)];
// Assemble local gradient
std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
- localStiffness_->assembleGradient(*it, localSolution, localGradient);
+ localStiffness_->assembleGradient(*it, basis_.getLocalFiniteElement(*it), localSolution, localGradient);
// Add to global gradient
for (int i=0; i<nDofs; i++)
- grad[indexSet.subIndex(*it,i,gridDim)] += localGradient[i];
+ grad[basis_.index(*it,i)] += localGradient[i];
}
}
-template <class GridView, class TargetSpace>
-double GeodesicFEAssembler<GridView, TargetSpace>::
+template <class Basis, class TargetSpace>
+double GeodesicFEAssembler<Basis, TargetSpace>::
computeEnergy(const std::vector<TargetSpace>& sol) const
{
double energy = 0;
- const typename GridView::IndexSet& indexSet = gridView_.indexSet();
-
- if (sol.size()!=indexSet.size(gridDim))
+ if (sol.size()!=basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
- Dune::array<TargetSpace,gridDim+1> localSolution;
-
- ElementIterator it = gridView_.template begin<0>();
- ElementIterator endIt = gridView_.template end<0>();
+ ElementIterator it = basis_.getGridView().template begin<0>();
+ ElementIterator endIt = basis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
+
+ // Number of degrees of freedom on this element
+ size_t nDofs = basis_.getLocalFiniteElement(*it).localBasis().size();
+
+ std::vector<TargetSpace> localSolution(nDofs);
- for (int i=0; i<it->template count<gridDim>(); i++)
- localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
+ for (int i=0; i<nDofs; i++)
+ localSolution[i] = sol[basis_.index(*it,i)];
- energy += localStiffness_->energy(*it, localSolution);
+ energy += localStiffness_->energy(*it, basis_.getLocalFiniteElement(*it), localSolution);
}
diff --git a/dune/gfe/geodesicfefunctionadaptor.hh b/dune/gfe/geodesicfefunctionadaptor.hh
index 6d152d161b2171737371378424608491cc9d2554..48a26f3f04da9daca1630946868a9648582602b4 100644
--- a/dune/gfe/geodesicfefunctionadaptor.hh
+++ b/dune/gfe/geodesicfefunctionadaptor.hh
@@ -13,6 +13,8 @@ void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
{
const int dim = GridType::dimension;
+ assert(x.size() == grid.size(dim));
+
typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
typedef typename GridType::template Codim<dim>::LeafIterator VertexIterator;
@@ -44,6 +46,7 @@ void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
// Restore and interpolate the data
// /////////////////////////////////////////////////////
+ P1NodalBasis<typename GridType::LeafGridView> p1Basis(grid.leafView());
x.resize(grid.size(dim));
ElementIterator eIt = grid.template leafbegin<0>();
@@ -57,7 +60,9 @@ void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
for (int i=0; i<eIt->father()->template count<dim>(); i++)
coefficients[i] = dofMap.find(idSet.subId(*eIt->father(),i,dim))->second;
- LocalGeodesicFEFunction<dim,double,TargetSpace> fatherFunction(coefficients);
+ typedef typename P1NodalBasis<typename GridType::LeafGridView>::LocalFiniteElement LocalFiniteElement;
+ LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fatherFunction(p1Basis.getLocalFiniteElement(*eIt),
+ coefficients);
// The embedding of this element into the father geometry
const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
@@ -84,4 +89,127 @@ void geodesicFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
}
+
+
+/** \brief Refine a grid globally and prolong a given geodesic finite element function
+ */
+template <class GridType, class TargetSpace>
+void higherOrderGFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+{
+ const int dim = GridType::dimension;
+
+ typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
+
+ // /////////////////////////////////////////////////////
+ // Save leaf p1 data in a map
+ // /////////////////////////////////////////////////////
+
+ const typename GridType::Traits::LocalIdSet& idSet = grid.localIdSet();
+
+ // DUNE ids are not unique across all codimensions, hence the following hack. Sigh...
+ typedef std::pair<typename GridType::Traits::LocalIdSet::IdType, unsigned int> IdType;
+ std::map<IdType, TargetSpace> dofMap;
+
+ typedef P2NodalBasis<typename GridType::LeafGridView,double> P2Basis;
+ P2Basis p2Basis(grid.leafView());
+
+ assert(x.size() == p2Basis.size());
+
+ ElementIterator eIt = grid.template leafbegin<0>();
+ ElementIterator eEndIt = grid.template leafend<0>();
+
+ for (; eIt!=eEndIt; ++eIt) {
+
+ const typename P2Basis::LocalFiniteElement& lfe = p2Basis.getLocalFiniteElement(*eIt);
+ //localCoefficients = p2Basis.getLocalFiniteElement(*eIt).localCoefficients();
+
+ for (size_t i=0; i<lfe.localCoefficients().size(); i++) {
+
+ IdType id = std::make_pair(idSet.subId(*eIt,
+ lfe.localCoefficients().localKey(i).subEntity(),
+ lfe.localCoefficients().localKey(i).codim()),
+ lfe.localCoefficients().localKey(i).codim());
+
+ unsigned int idx = p2Basis.index(*eIt, i);
+
+ //std::cout << "id: (" << id.first << ", " << id.second << ")" << std::endl;
+ dofMap.insert(std::make_pair(id, x[idx]));
+
+ }
+
+ }
+
+
+ // /////////////////////////////////////////////////////
+ // Globally refine the grid
+ // /////////////////////////////////////////////////////
+
+ grid.globalRefine(1);
+
+
+ // /////////////////////////////////////////////////////
+ // Restore and interpolate the data
+ // /////////////////////////////////////////////////////
+
+ p2Basis.update(grid.leafView());
+
+ x.resize(p2Basis.size());
+
+ for (eIt=grid.template leafbegin<0>(); eIt!=eEndIt; ++eIt) {
+
+ const typename P2Basis::LocalFiniteElement& lfe = p2Basis.getLocalFiniteElement(*eIt);
+
+ std::auto_ptr<typename Dune::PQkLocalFiniteElementFactory<double,double,dim,2>::FiniteElementType> fatherLFE
+ = std::auto_ptr<typename Dune::PQkLocalFiniteElementFactory<double,double,dim,2>::FiniteElementType>(Dune::PQkLocalFiniteElementFactory<double,double,dim,2>::create(eIt->type()));
+
+ // Set up a local gfe function on the father element
+ std::vector<TargetSpace> coefficients(fatherLFE->localCoefficients().size());
+
+ for (int i=0; i<fatherLFE->localCoefficients().size(); i++) {
+
+ IdType id = std::make_pair(idSet.subId(*eIt->father(),
+ fatherLFE->localCoefficients().localKey(i).subEntity(),
+ fatherLFE->localCoefficients().localKey(i).codim()),
+ fatherLFE->localCoefficients().localKey(i).codim());
+
+ coefficients[i] = dofMap.find(id)->second;
+
+ }
+
+ LocalGeodesicFEFunction<dim,double,typename P2Basis::LocalFiniteElement,TargetSpace> fatherFunction(*fatherLFE, coefficients);
+
+ // The embedding of this element into the father geometry
+ const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
+
+ for (int i=0; i<lfe.localCoefficients().size(); i++) {
+
+ IdType id = std::make_pair(idSet.subId(*eIt,
+ lfe.localCoefficients().localKey(i).subEntity(),
+ lfe.localCoefficients().localKey(i).codim()),
+ lfe.localCoefficients().localKey(i).codim());
+
+ unsigned int idx = p2Basis.index(*eIt, i);
+
+ if (dofMap.find(id) != dofMap.end()) {
+
+ // If the vertex exists on the coarser level we take the value from there.
+ // That should be faster and more accurate than interpolating
+ x[idx] = dofMap[id];
+
+ } else {
+
+ // Interpolate
+ const Dune::GenericReferenceElement<double,dim>& refElement = Dune::GenericReferenceElements<double,dim>::general(eIt->type());
+ const Dune::FieldVector<double,dim>& pos = refElement.position(lfe.localCoefficients().localKey(i).subEntity(),
+ lfe.localCoefficients().localKey(i).codim());
+ x[idx] = fatherFunction.evaluate(geometryInFather.global(pos));
+
+ }
+
+ }
+
+ }
+
+}
+
#endif
diff --git a/dune/gfe/harmonicenergystiffness.hh b/dune/gfe/harmonicenergystiffness.hh
index 2ca8e2e575a18ae2953ae36612223923c1072a86..266ec1e2897fda8797726ca73eb3d3cf4b6b23c6 100644
--- a/dune/gfe/harmonicenergystiffness.hh
+++ b/dune/gfe/harmonicenergystiffness.hh
@@ -11,9 +11,9 @@
#warning Finite-difference approximation of the energy gradient
#endif
-template<class GridView, class TargetSpace>
+template<class GridView, class LocalFiniteElement, class TargetSpace>
class HarmonicEnergyLocalStiffness
- : public LocalGeodesicFEStiffness<GridView,TargetSpace>
+ : public LocalGeodesicFEStiffness<GridView,LocalFiniteElement,TargetSpace>
{
// grid types
typedef typename GridView::Grid::ctype DT;
@@ -30,7 +30,8 @@ public:
/** \brief Assemble the energy for a single element */
RT energy (const Entity& e,
- const Dune::array<TargetSpace,gridDim+1>& localSolution) const;
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution) const;
#ifndef HARMONIC_ENERGY_FD_GRADIENT
// The finite difference gradient method is in the base class.
@@ -46,16 +47,18 @@ public:
#endif
};
-template <class GridView, class TargetSpace>
-typename HarmonicEnergyLocalStiffness<GridView, TargetSpace>::RT HarmonicEnergyLocalStiffness<GridView, TargetSpace>::
+template <class GridView, class LocalFiniteElement, class TargetSpace>
+typename HarmonicEnergyLocalStiffness<GridView, LocalFiniteElement, TargetSpace>::RT
+HarmonicEnergyLocalStiffness<GridView, LocalFiniteElement, TargetSpace>::
energy(const Entity& element,
- const Dune::array<TargetSpace,gridDim+1>& localSolution) const
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution) const
{
- RT energy = 0;
-
- assert(element.type().isSimplex());
+ assert(element.type() == localFiniteElement.type());
- LocalGeodesicFEFunction<gridDim, double, TargetSpace> localGeodesicFEFunction(localSolution);
+ RT energy = 0;
+ LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
+ localSolution);
int quadOrder = 1;//gridDim;
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 663cbc0367cbb81646b52a46d8f1a0436aba7641..43ee577b77a302046b63322f0c7e4beaa16025e8 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -20,7 +20,7 @@
\tparam ctype Type used for coordinates on the reference element
\tparam TargetSpace The manifold that the function takes its values in
*/
-template <int dim, class ctype, class TargetSpace>
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
class LocalGeodesicFEFunction
{
@@ -32,8 +32,22 @@ class LocalGeodesicFEFunction
public:
/** \brief Constructor */
- LocalGeodesicFEFunction(const Dune::array<TargetSpace,dim+1>& coefficients)
- : coefficients_(coefficients)
+ LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
+ const Dune::array<TargetSpace,dim+1>& coefficients) DUNE_DEPRECATED
+ : localFiniteElement_(localFiniteElement),
+ coefficients_(coefficients.size())
+ {
+ for (size_t i=0; i<coefficients.size(); i++)
+ coefficients_[i] = coefficients[i];
+ }
+
+ /** \brief Constructor
+ \param type Type of the reference element
+ */
+ LocalGeodesicFEFunction(const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& coefficients)
+ : localFiniteElement_(localFiniteElement),
+ coefficients_(coefficients)
{}
/** \brief Evaluate the function */
@@ -67,28 +81,6 @@ public:
private:
- /** \brief The linear part of the map that turns coordinates on the reference simplex into coordinates on the standard simplex
- \todo A special-purpose implementation of this matrix may lead to some speed-up */
- static Dune::FieldMatrix<ctype,dim+1,dim> referenceToBarycentricLinearPart()
- {
- Dune::FieldMatrix<ctype,dim+1,dim> B;
- B[0] = -1;
- for (int i=0; i<dim; i++)
- for (int j=0; j<dim; j++)
- B[i+1][j] = (i==j);
- return B;
- }
-
- 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];
- result[i+1] = local[i];
- }
- return result;
- }
-
static Dune::FieldMatrix<double,embeddedDim,embeddedDim> pseudoInverse(const Dune::FieldMatrix<double,embeddedDim,embeddedDim>& dFdq,
const TargetSpace& q)
{
@@ -122,10 +114,10 @@ private:
}
/** \brief Compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w */
- Dune::FieldMatrix<ctype,embeddedDim,dim+1> computeDFdw(const TargetSpace& q) const
+ Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > computeDFdw(const TargetSpace& q) const
{
- Dune::FieldMatrix<ctype,embeddedDim,dim+1> dFdw;
- for (int i=0; i<dim+1; i++) {
+ Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > dFdw(embeddedDim,localFiniteElement_.localBasis().size());
+ for (int i=0; i<localFiniteElement_.localBasis().size(); i++) {
Dune::FieldVector<ctype,embeddedDim> tmp = TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], q);
for (int j=0; j<embeddedDim; j++)
dFdw[j][i] = tmp[j];
@@ -133,27 +125,38 @@ private:
return dFdw;
}
- Tensor3<ctype,embeddedDim,embeddedDim,embeddedDim> computeDqDqF(const Dune::array<ctype,dim+1>& w, const TargetSpace& q) const
+ Tensor3<ctype,embeddedDim,embeddedDim,embeddedDim> computeDqDqF(const std::vector<ctype>& w, const TargetSpace& q) const
{
Tensor3<ctype,embeddedDim,embeddedDim,embeddedDim> result;
result = 0;
- for (int i=0; i<dim+1; i++)
+ for (int i=0; i<w.size(); i++)
result.axpy(w[i], TargetSpace::thirdDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i],q));
return result;
}
+ /** \brief The scalar local finite element, which provides the weighting factors
+ \todo We really only need the local basis
+ */
+ const LocalFiniteElement& localFiniteElement_;
+
/** \brief The coefficient vector */
- Dune::array<TargetSpace,dim+1> coefficients_;
+ std::vector<TargetSpace> coefficients_;
};
-template <int dim, class ctype, class TargetSpace>
-TargetSpace LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+TargetSpace LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
- // First compute the coordinates on the standard simplex (in R^{n+1})
- Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
+ // Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
+ std::vector<Dune::FieldVector<ctype,1> > wNested;
+ localFiniteElement_.localBasis().evaluateFunction(local,wNested);
+ std::vector<ctype> w(wNested.size());
+ for (size_t i=0; i<w.size(); i++)
+ w[i] = wNested[i][0];
+
+ /** \todo The 'dim+1' parameter is a dummy here */
AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
TargetSpaceRiemannianTRSolver<TargetSpace,dim+1> solver;
@@ -172,8 +175,8 @@ evaluate(const Dune::FieldVector<ctype, dim>& local) const
return solver.getSol();
}
-template <int dim, class ctype, class TargetSpace>
-Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::dimension, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::dimension, dim> LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
{
Dune::FieldMatrix<ctype, embeddedDim, dim> result;
@@ -198,18 +201,29 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
TargetSpace q = evaluate(local);
// the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
- Dune::FieldMatrix<ctype,dim+1,dim> B = referenceToBarycentricLinearPart();
-
+ std::vector<Dune::FieldMatrix<ctype,1,dim> > BNested(coefficients_.size());
+ localFiniteElement_.localBasis().evaluateJacobian(local, BNested);
+ Dune::Matrix<Dune::FieldMatrix<double,1,1> > B(coefficients_.size(), dim);
+ for (size_t i=0; i<coefficients_.size(); i++)
+ for (size_t j=0; j<dim; j++)
+ B[i][j] = BNested[i][0][j];
+
// compute negative derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
- Dune::FieldMatrix<ctype,embeddedDim,dim+1> dFdw = computeDFdw(q);
+ Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > dFdw = computeDFdw(q);
dFdw *= -1;
// multiply the two previous matrices: the result is the right hand side
- Dune::FieldMatrix<ctype,embeddedDim,dim> RHS = dFdw * B;
+ Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > RHS = dFdw * B;
// the actual system matrix
- Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
- AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
+ std::vector<Dune::FieldVector<ctype,1> > wNested;
+ localFiniteElement_.localBasis().evaluateFunction(local, wNested);
+
+ std::vector<ctype> w(wNested.size());
+ for (size_t i=0; i<w.size(); i++)
+ w[i] = wNested[i][0];
+
+ AverageDistanceAssembler<TargetSpace,1> assembler(coefficients_, w);
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleEmbeddedHessian(q,dFdq);
@@ -263,8 +277,8 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local) const
return result;
}
-template <int dim, class ctype, class TargetSpace>
-Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::dimension, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::dimension, dim> LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
{
double eps = 1e-6;
@@ -290,8 +304,8 @@ evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
return result;
}
-template <int dim, class ctype, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
Dune::FieldMatrix<double,embeddedDim,embeddedDim>& result) const
@@ -339,8 +353,8 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
}
-template <int dim, class ctype, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
Dune::FieldMatrix<double,embeddedDim,embeddedDim>& result) const
@@ -365,8 +379,8 @@ evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& l
cornersPlus [coefficient] = TargetSpace::exp(coefficients_[coefficient], forwardVariation);
cornersMinus[coefficient] = TargetSpace::exp(coefficients_[coefficient], backwardVariation);
- LocalGeodesicFEFunction<dim,double,TargetSpace> fPlus(cornersPlus);
- LocalGeodesicFEFunction<dim,double,TargetSpace> fMinus(cornersMinus);
+ LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fPlus(cornersPlus);
+ LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fMinus(cornersMinus);
TargetSpace hPlus = fPlus.evaluate(local);
TargetSpace hMinus = fMinus.evaluate(local);
@@ -387,8 +401,8 @@ evaluateFDDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& l
}
-template <int dim, class ctype, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
Tensor3<double,embeddedDim,embeddedDim,dim>& result) const
@@ -397,13 +411,24 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
TargetSpace q = evaluate(local);
// the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
- Dune::FieldMatrix<ctype,dim+1,dim> B = referenceToBarycentricLinearPart();
+ std::vector<Dune::FieldMatrix<ctype,1,dim> > BNested(coefficients_.size());
+ localFiniteElement_.localBasis().evaluateJacobian(local, BNested);
+ Dune::Matrix<Dune::FieldMatrix<double,1,1> > B(coefficients_.size(), dim);
+ for (size_t i=0; i<coefficients_.size(); i++)
+ for (size_t j=0; j<dim; j++)
+ B[i][j] = BNested[i][0][j];
- // compute derivate of F(w,q) (the derivative of the weighted distance fctl) wrt to w
- Dune::FieldMatrix<ctype,embeddedDim,dim+1> dFdw = computeDFdw(q);
+ // compute derivative of F(w,q) (the derivative of the weighted distance fctl) wrt to w
+ Dune::Matrix<Dune::FieldMatrix<ctype,1,1> > dFdw = computeDFdw(q);
// the actual system matrix
- Dune::array<ctype,dim+1> w = barycentricCoordinates(local);
+ std::vector<Dune::FieldVector<ctype,1> > wNested;
+ localFiniteElement_.localBasis().evaluateFunction(local,wNested);
+
+ std::vector<ctype> w(wNested.size());
+ for (size_t i=0; i<w.size(); i++)
+ w[i] = wNested[i][0];
+
AverageDistanceAssembler<TargetSpace,dim+1> assembler(coefficients_, w);
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
@@ -468,8 +493,8 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
}
-template <int dim, class ctype, class TargetSpace>
-void LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
+template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
int coefficient,
Tensor3<double,embeddedDim,embeddedDim,dim>& result) const
@@ -485,8 +510,8 @@ evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>
aMinus[j] -= eps;
cornersPlus[coefficient] = TargetSpace(aPlus);
cornersMinus[coefficient] = TargetSpace(aMinus);
- LocalGeodesicFEFunction<dim,double,TargetSpace> fPlus(cornersPlus);
- LocalGeodesicFEFunction<dim,double,TargetSpace> fMinus(cornersMinus);
+ LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fPlus(cornersPlus);
+ LocalGeodesicFEFunction<dim,double,LocalFiniteElement,TargetSpace> fMinus(cornersMinus);
Dune::FieldMatrix<double,TargetSpace::EmbeddedTangentVector::size,dim> hPlus = fPlus.evaluateDerivative(local);
Dune::FieldMatrix<double,TargetSpace::EmbeddedTangentVector::size,dim> hMinus = fMinus.evaluateDerivative(local);
@@ -523,8 +548,8 @@ evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>
\tparam dim Dimension of the reference element
\tparam ctype Type used for coordinates on the reference element
*/
-template <int dim, class ctype>
-class LocalGeodesicFEFunction<dim,ctype,RigidBodyMotion<3> >
+template <int dim, class ctype, class LocalFiniteElement>
+class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,RigidBodyMotion<3> >
{
typedef RigidBodyMotion<3> TargetSpace;
@@ -545,7 +570,7 @@ public:
for (int i=0; i<dim+1; i++)
orientationCoefficients[i] = coefficients[i].q;
- orientationFEFunction_ = std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,Rotation<3,double> > > (new LocalGeodesicFEFunction<dim,ctype,Rotation<3,double> >(orientationCoefficients));
+ orientationFEFunction_ = std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> > > (new LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> >(orientationCoefficients));
}
@@ -665,7 +690,7 @@ private:
Dune::array<Dune::FieldVector<double,3>, dim+1> translationCoefficients_;
- std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,Rotation<3,double> > > orientationFEFunction_;
+ std::auto_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<3,double> > > orientationFEFunction_;
};
diff --git a/dune/gfe/localgeodesicfestiffness.hh b/dune/gfe/localgeodesicfestiffness.hh
index 36e1a52dd0073004d1e3107868d6e8dc50ba1aa5..a307574dc871c2744bb83512a0ac3fc8acd2ab14 100644
--- a/dune/gfe/localgeodesicfestiffness.hh
+++ b/dune/gfe/localgeodesicfestiffness.hh
@@ -7,7 +7,7 @@
#include <dune/istl/matrix.hh>
-template<class GridView, class TargetSpace>
+template<class GridView, class LocalFiniteElement, class TargetSpace>
class LocalGeodesicFEStiffness
{
// grid types
@@ -40,17 +40,20 @@ public:
*/
virtual void assembleHessian(const Entity& e,
- const Dune::array<TargetSpace,gridDim+1>& localSolution);
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution);
/** \brief Compute the energy at the current configuration */
virtual RT energy (const Entity& e,
- const Dune::array<TargetSpace,gridDim+1>& localSolution) const = 0;
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& 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 Dune::array<TargetSpace,gridDim+1>& solution,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& solution,
std::vector<typename TargetSpace::TangentVector>& gradient) const;
// assembled data
@@ -59,10 +62,11 @@ public:
};
-template <class GridView, class TargetSpace>
-void LocalGeodesicFEStiffness<GridView, TargetSpace>::
+template <class GridView, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEStiffness<GridView, LocalFiniteElement, TargetSpace>::
assembleGradient(const Entity& element,
- const Dune::array<TargetSpace,gridDim+1>& localSolution,
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution,
std::vector<typename TargetSpace::TangentVector>& localGradient) const
{
@@ -74,8 +78,8 @@ assembleGradient(const Entity& element,
localGradient.resize(localSolution.size());
- Dune::array<TargetSpace,gridDim+1> forwardSolution = localSolution;
- Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
for (size_t i=0; i<localSolution.size(); i++) {
@@ -93,7 +97,7 @@ assembleGradient(const Entity& element,
forwardSolution[i] = TargetSpace::exp(localSolution[i], forwardCorrection);
backwardSolution[i] = TargetSpace::exp(localSolution[i], backwardCorrection);
- localGradient[i][j] = (energy(element,forwardSolution) - energy(element,backwardSolution)) / (2*eps);
+ localGradient[i][j] = (energy(element,localFiniteElement,forwardSolution) - energy(element,localFiniteElement, backwardSolution)) / (2*eps);
}
@@ -108,22 +112,20 @@ assembleGradient(const Entity& element,
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
-template <class GridType, class TargetSpace>
-void LocalGeodesicFEStiffness<GridType, TargetSpace>::
+template <class GridType, class LocalFiniteElement, class TargetSpace>
+void LocalGeodesicFEStiffness<GridType, LocalFiniteElement, TargetSpace>::
assembleHessian(const Entity& element,
- const Dune::array<TargetSpace,gridDim+1>& localSolution)
+ const LocalFiniteElement& localFiniteElement,
+ const std::vector<TargetSpace>& localSolution)
{
- // 1 degree of freedom per element vertex
- size_t nDofs = element.template count<gridDim>();
+ // Number of degrees of freedom for this element
+ size_t nDofs = localSolution.size();
// Clear assemble data
A_.setSize(nDofs, nDofs);
A_ = 0;
-
-
-
const double eps = 1e-4;
std::vector<Dune::FieldMatrix<double,blocksize,embeddedBlocksize> > B(localSolution.size());
@@ -132,7 +134,7 @@ assembleHessian(const Entity& element,
// Precompute negative energy at the current configuration
// (negative because that is how we need it as part of the 2nd-order fd formula)
- double centerValue = -energy(element, localSolution);
+ double centerValue = -energy(element, localFiniteElement, localSolution);
// Precompute energy infinitesimal corrections in the directions of the local basis vectors
std::vector<Dune::array<double,blocksize> > forwardEnergy(nDofs);
@@ -144,14 +146,14 @@ assembleHessian(const Entity& element,
Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi;
minusEpsXi *= -1;
- Dune::array<TargetSpace,gridDim+1> forwardSolution = localSolution;
- Dune::array<TargetSpace,gridDim+1> backwardSolution = localSolution;
+ std::vector<TargetSpace> forwardSolution = localSolution;
+ std::vector<TargetSpace> backwardSolution = localSolution;
forwardSolution[i] = TargetSpace::exp(localSolution[i],epsXi);
backwardSolution[i] = TargetSpace::exp(localSolution[i],minusEpsXi);
- forwardEnergy[i][i2] = energy(element, forwardSolution);
- backwardEnergy[i][i2] = energy(element, backwardSolution);
+ forwardEnergy[i][i2] = energy(element, localFiniteElement, forwardSolution);
+ backwardEnergy[i][i2] = energy(element, localFiniteElement, backwardSolution);
}
@@ -169,8 +171,8 @@ assembleHessian(const Entity& element,
Dune::FieldVector<double,embeddedBlocksize> minusEpsXi = epsXi; minusEpsXi *= -1;
Dune::FieldVector<double,embeddedBlocksize> minusEpsEta = epsEta; minusEpsEta *= -1;
- Dune::array<TargetSpace,gridDim+1> forwardSolutionXiEta = localSolution;
- Dune::array<TargetSpace,gridDim+1> backwardSolutionXiEta = localSolution;
+ std::vector<TargetSpace> forwardSolutionXiEta = localSolution;
+ std::vector<TargetSpace> backwardSolutionXiEta = localSolution;
if (i==j)
forwardSolutionXiEta[i] = TargetSpace::exp(localSolution[i],epsXi+epsEta);
@@ -186,8 +188,8 @@ assembleHessian(const Entity& element,
backwardSolutionXiEta[j] = TargetSpace::exp(localSolution[j],minusEpsEta);
}
- double forwardValue = energy(element, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
- double backwardValue = energy(element, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
+ double forwardValue = energy(element, localFiniteElement, forwardSolutionXiEta) - forwardEnergy[i][i2] - forwardEnergy[j][j2];
+ double backwardValue = energy(element, localFiniteElement, backwardSolutionXiEta) - backwardEnergy[i][i2] - backwardEnergy[j][j2];
A_[i][j][i2][j2] = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index d9f20e12d5bec4c8169c0f282021c28bc1bbbe0c..19b327f858274513d49ad0d030eea0a874fcb3a9 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -12,6 +12,7 @@
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
+#include <dune/solvers/transferoperators/p2top1mgtransfer.hh>
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include "maxnormtrustregion.hh"
@@ -23,7 +24,11 @@
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::
setup(const GridType& grid,
- const GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace>* assembler,
+#ifdef HIGHER_ORDER
+ const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
+#else
+ const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
+#endif
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -123,9 +128,19 @@ setup(const GridType& grid,
// //////////////////////////////////////////////////////////
hasObstacle_.resize(numLevels);
+#ifdef HIGHER_ORDER
+ P2NodalBasis<typename GridType::LeafGridView,double> p2Basis(grid_->leafView());
+ P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafView());
+
+ hasObstacle_.back().resize(p2Basis.size(), true);
+
+ for (int i=0; i<hasObstacle_.size()-1; i++)
+ hasObstacle_[i].resize(grid_->size(i+1, gridDim),true);
+#else
for (int i=0; i<hasObstacle_.size(); i++)
hasObstacle_[i].resize(grid_->size(i, gridDim),true);
-
+#endif
+
// ////////////////////////////////////
// Create the transfer operators
// ////////////////////////////////////
@@ -134,12 +149,27 @@ setup(const GridType& grid,
mmgStep->mgTransfer_.resize(numLevels-1);
+#ifdef HIGHER_ORDER
+ if (numLevels>1) {
+ P2toP1MGTransfer<CorrectionType>* topTransferOp = new P2toP1MGTransfer<CorrectionType>;
+ topTransferOp->setup(p2Basis,p1Basis);
+ mmgStep->mgTransfer_.back() = topTransferOp;
+
+ for (int i=0; i<mmgStep->mgTransfer_.size()-1; i++){
+ TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
+ newTransferOp->setup(*grid_,i+1,i+2);
+ mmgStep->mgTransfer_[i] = newTransferOp;
+ }
+
+ }
+
+#else
for (int i=0; i<mmgStep->mgTransfer_.size(); i++){
TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
newTransferOp->setup(*grid_,i,i+1);
mmgStep->mgTransfer_[i] = newTransferOp;
}
-
+#endif
}
@@ -180,6 +210,10 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
for (int i=0; i<maxTrustRegionSteps_; i++) {
+/* std::cout << "current iterate:\n";
+ for (int j=0; j<x_.size(); j++)
+ std::cout << x_[j] << std::endl;*/
+
Dune::Timer totalTimer;
if (this->verbosity_ == Solver::FULL) {
std::cout << "----------------------------------------------------" << std::endl;
@@ -206,7 +240,9 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
// The right hand side is the _negative_ gradient
rhs *= -1;
-
+/* std::cout << "rhs:\n" << rhs << std::endl;
+ std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;*/
+
// //////////////////////////////////////////////////////////////////////
// Modify matrix and right-hand side to account for Dirichlet values
// //////////////////////////////////////////////////////////////////////
@@ -346,11 +382,14 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
for (int j=0; j<newIterate.size(); j++)
newIterate[j] = TargetSpace::exp(newIterate[j], corr[j]);
#else
+ //std::cout << "embedded correction:\n";
for (int j=0; j<newIterate.size(); j++) {
Dune::FieldMatrix<double,TargetSpace::TangentVector::dimension,TargetSpace::EmbeddedTangentVector::dimension> B = x_[j].orthonormalFrame();
Dune::FieldVector<double,TargetSpace::EmbeddedTangentVector::dimension> embeddedCorr(0);
+ //std::cout << "B[" << j << "]:\n" << B << std::endl;
B.mtv(corr[j], embeddedCorr);
newIterate[j] = TargetSpace::exp(newIterate[j], embeddedCorr);
+ //std::cout << embeddedCorr << " " << newIterate[j] << std::endl;
}
#endif
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 6647e9f9ad8e420daec4fd73d4bf45bdd3183fcb..9912e07dd2d7403aaf84ea3f25e14cf79780860b 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -42,7 +42,11 @@ public:
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
- const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* assembler,
+#ifdef HIGHER_ORDER
+ const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
+#else
+ const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler,
+#endif
const SolutionType& x,
const Dune::BitSetVector<blocksize>& dirichletNodes,
double tolerance,
@@ -97,7 +101,11 @@ protected:
std::auto_ptr<MatrixType> hessianMatrix_;
/** \brief The assembler for the material law */
- const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* assembler_;
+#ifdef HIGHER_ORDER
+ const GeodesicFEAssembler<P2NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler_;
+#else
+ const GeodesicFEAssembler<P1NodalBasis<typename GridType::LeafGridView,double>, TargetSpace>* assembler_;
+#endif
/** \brief The solver for the quadratic inner problems */
std::shared_ptr<Solver> innerSolver_;