From 85d3abe3c98348463ed104c27d5c95c1d38c7150 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 19 Nov 2011 20:07:49 +0000
Subject: [PATCH] geodesicfefunctionadaptor can now handle spaces of higher
than second order
[[Imported from SVN: r8240]]
---
dune/gfe/geodesicfefunctionadaptor.hh | 41 ++++++++++++---------------
harmonicmaps-eoc.cc | 15 ++++++----
2 files changed, 28 insertions(+), 28 deletions(-)
diff --git a/dune/gfe/geodesicfefunctionadaptor.hh b/dune/gfe/geodesicfefunctionadaptor.hh
index 6a504d50..72f6cc00 100644
--- a/dune/gfe/geodesicfefunctionadaptor.hh
+++ b/dune/gfe/geodesicfefunctionadaptor.hh
@@ -4,8 +4,6 @@
#include <vector>
#include <map>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
-
#include "localgeodesicfefunction.hh"
template <class Basis, class TargetSpace>
@@ -118,10 +116,15 @@ struct CoordinateFunction
/** \brief Refine a grid globally and prolong a given geodesic finite element function
+ *
+ * \tparam order Interpolation order of the function space. Kinda stupid that I
+ * have to provide this by hand. Will change...
*/
-template <class GridType>
-static void higherOrderGFEFunctionAdaptor(GridType& grid, std::vector<TargetSpace>& x)
+template <int order>
+static void higherOrderGFEFunctionAdaptor(Basis& basis,
+ typename Basis::GridView::Grid& grid, std::vector<TargetSpace>& x)
{
+ typedef typename Basis::GridView::Grid GridType;
const int dim = GridType::dimension;
typedef typename GridType::template Codim<0>::LeafIterator ElementIterator;
@@ -135,25 +138,18 @@ static void higherOrderGFEFunctionAdaptor(GridType& grid, std::vector<TargetSpac
typedef typename GridType::Traits::LocalIdSet::IdType IdType;
std::map<IdType, std::vector<TargetSpace> > dofMap;
- typedef P2NodalBasis<typename GridType::LeafGridView,double> P2Basis;
- P2Basis p2Basis(grid.leafView());
-
- assert(x.size() == p2Basis.size());
+ assert(x.size() == basis.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);
+ const typename Basis::LocalFiniteElement& lfe = basis.getLocalFiniteElement(*eIt);
std::vector<TargetSpace> coefficients(lfe.localCoefficients().size());
- for (size_t i=0; i<lfe.localCoefficients().size(); i++) {
-
- unsigned int idx = p2Basis.index(*eIt, i);
- coefficients[i] = x[idx];
-
- }
+ for (size_t i=0; i<lfe.localCoefficients().size(); i++)
+ coefficients[i] = x[basis.index(*eIt, i)];
IdType id = idSet.id(*eIt);
dofMap.insert(std::make_pair(id, coefficients));
@@ -172,22 +168,21 @@ static void higherOrderGFEFunctionAdaptor(GridType& grid, std::vector<TargetSpac
// Restore and interpolate the data
// /////////////////////////////////////////////////////
- p2Basis.update(grid.leafView());
+ basis.update(grid.leafView());
- x.resize(p2Basis.size());
+ x.resize(basis.size());
for (eIt=grid.template leafbegin<0>(); eIt!=eEndIt; ++eIt) {
- const typename P2Basis::LocalFiniteElement& lfe = p2Basis.getLocalFiniteElement(*eIt);
+ const typename Basis::LocalFiniteElement& lfe = basis.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()));
+ std::auto_ptr<typename Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::FiniteElementType> fatherLFE
+ = std::auto_ptr<typename Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::FiniteElementType>(Dune::PQkLocalFiniteElementFactory<double,double,dim,order>::create(eIt->father()->type()));
// Set up a local gfe function on the father element
- //std::vector<TargetSpace> coefficients(fatherLFE->localCoefficients().size());
std::vector<TargetSpace> coefficients = dofMap[idSet.id(*eIt->father())];
- LocalGeodesicFEFunction<dim,double,typename P2Basis::LocalFiniteElement,TargetSpace> fatherFunction(*fatherLFE, coefficients);
+ LocalGeodesicFEFunction<dim,double,typename Basis::LocalFiniteElement,TargetSpace> fatherFunction(*fatherLFE, coefficients);
// The embedding of this element into the father geometry
const typename GridType::template Codim<0>::LocalGeometry& geometryInFather = eIt->geometryInFather();
@@ -207,7 +202,7 @@ static void higherOrderGFEFunctionAdaptor(GridType& grid, std::vector<TargetSpac
for (int i=0; i<lfe.localCoefficients().size(); i++) {
- unsigned int idx = p2Basis.index(*eIt, i);
+ unsigned int idx = basis.index(*eIt, i);
x[idx] = fatherFunction.evaluate(geometryInFather.global(lagrangeNodes[i]));
diff --git a/harmonicmaps-eoc.cc b/harmonicmaps-eoc.cc
index d48dff07..a3650329 100644
--- a/harmonicmaps-eoc.cc
+++ b/harmonicmaps-eoc.cc
@@ -223,8 +223,8 @@ int main (int argc, char *argv[]) try
#else
typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
#endif
- FEBasis basis(referenceGrid->leafView());
- OperatorAssembler<FEBasis,FEBasis> operatorAssembler(basis, basis);
+ FEBasis referenceBasis(referenceGrid->leafView());
+ OperatorAssembler<FEBasis,FEBasis> operatorAssembler(referenceBasis, referenceBasis);
LaplaceAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> laplaceLocalAssembler;
MassAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> massMatrixLocalAssembler;
@@ -277,12 +277,17 @@ int main (int argc, char *argv[]) try
#endif
// Prolong solution to the very finest grid
- for (int j=i; j<numLevels; j++)
-#if defined THIRD_ORDER || defined SECOND_ORDER
- higherOrderGFEFunctionAdaptor(*grid, solution);
+ for (int j=i; j<numLevels; j++) {
+ FEBasis basis(grid->leafView());
+#if defined THIRD_ORDER
+ GeodesicFEFunctionAdaptor<FEBasis,TargetSpace>::higherOrderGFEFunctionAdaptor<3>(basis, *grid, solution);
+#elif defined SECOND_ORDER
+ GeodesicFEFunctionAdaptor<FEBasis,TargetSpace>::higherOrderGFEFunctionAdaptor<2>(basis, *grid, solution);
#else
geodesicFEFunctionAdaptor(*grid, solution);
#endif
+ }
+
// Interpret TargetSpace as isometrically embedded into an R^m, because this is
// how the corresponding Sobolev spaces are defined.
--
GitLab