diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index 2d9e4802931cd7b96a6b7c3f063793c4b810e318..0efff862ffd04972b9fd32de277df72e676de414 100644
--- a/dune/gfe/Makefile.am
+++ b/dune/gfe/Makefile.am
@@ -26,7 +26,6 @@ srcinclude_HEADERS = adolcnamespaceinjections.hh \
localprojectedfefunction.hh \
maxnormtrustregion.hh \
orthogonalmatrix.hh \
- pktop1mgtransfer.hh \
quaternion.hh \
realtuple.hh \
riemanniantrsolver.hh \
diff --git a/dune/gfe/pktop1mgtransfer.hh b/dune/gfe/pktop1mgtransfer.hh
deleted file mode 100644
index 0ea4f11865f32078d48a39e96f80906db26feee6..0000000000000000000000000000000000000000
--- a/dune/gfe/pktop1mgtransfer.hh
+++ /dev/null
@@ -1,150 +0,0 @@
-#ifndef PK_TO_P1_MG_TRANSFER_HH
-#define PK_TO_P1_MG_TRANSFER_HH
-
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/bitsetvector.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
-
-#include <dune/solvers/transferoperators/truncatedmgtransfer.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-
-/** \brief Restriction and prolongation operator for truncated multigrid
- * from a Pk nodal basis to a p1 one.
- *
- * \note This class is a temporary hack. Certainly there are better ways to do
- * this, and of course they should be in dune-solvers. But I need this quick
- * for some rapid prototyping.
- */
-
-template<
- class VectorType,
- class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension>,
- class MatrixType = Dune::BCRSMatrix< typename Dune::FieldMatrix<
- typename VectorType::field_type, VectorType::block_type::dimension, VectorType::block_type::dimension> > >
-class PKtoP1MGTransfer :
- public TruncatedCompressedMGTransfer<VectorType, BitVectorType, MatrixType>
-{
-
- enum {blocksize = VectorType::block_type::dimension};
-
- typedef typename VectorType::field_type field_type;
-
-public:
-
- typedef typename CompressedMultigridTransfer<VectorType, BitVectorType, MatrixType>::TransferOperatorType TransferOperatorType;
-
-
- /** \brief Default constructor */
- PKtoP1MGTransfer()
- {}
-
- /** \brief Constructor with a given transfer matrix
- */
- PKtoP1MGTransfer(std::shared_ptr<TransferOperatorType> matrix)
- : TruncatedCompressedMGTransfer<VectorType,BitVectorType,MatrixType>(matrix)
- {}
-
-
- template <class Basis, class GridView>
- void setup(const Basis& fineBasis,
- const P1NodalBasis<GridView>& p1Basis)
- {
- typedef typename TransferOperatorType::block_type TransferMatrixBlock;
-
- int rows = fineBasis.size();
- int cols = p1Basis.size();
-
- this->matrix_->setSize(rows,cols);
- this->matrix_->setBuildMode(TransferOperatorType::random);
-
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
-
- const GridView& gridView = p1Basis.getGridView();
-
- ElementIterator it = gridView.template begin<0>();
- ElementIterator endIt = gridView.template end<0>();
-
-
- // ///////////////////////////////////////////
- // Determine the indices present in the matrix
- // /////////////////////////////////////////////////
- Dune::MatrixIndexSet indices(rows, cols);
-
- for (; it != endIt; ++it) {
-
- typedef typename P1NodalBasis<GridView,double>::LocalFiniteElement LP1FE;
- const LP1FE& coarseBaseSet = p1Basis.getLocalFiniteElement(*it);
-
- typedef typename Dune::LocalFiniteElementFunctionBase<LP1FE>::type FunctionBaseClass;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- AlienElementLocalBasisFunction<Element, LP1FE, FunctionBaseClass> f(*it, *it, coarseBaseSet.localBasis());
-
- for (int i=0; i<coarseBaseSet.localBasis().size(); i++) {
-
- f.setIndex(i);
- std::vector<double> values;
- fineBasis.getLocalFiniteElement(*it).localInterpolation().interpolate(f,values);
-
- for (size_t j=0; j<values.size(); j++) {
-
- if (values[j] > 0.001)
- {
- size_t globalFine = fineBasis.index(*it, j);
- size_t globalCoarse = p1Basis.index(*it, i);
-
- indices.add(globalFine, globalCoarse);
- }
-
- }
-
- }
-
- }
-
- indices.exportIdx(*this->matrix_);
-
- // /////////////////////////////////////////////
- // Compute the matrix
- // /////////////////////////////////////////////
-
- for (it = gridView.template begin<0>(); it != endIt; ++it) {
-
- typedef typename P1NodalBasis<GridView,double>::LocalFiniteElement LP1FE;
- const LP1FE& coarseBaseSet = p1Basis.getLocalFiniteElement(*it);
-
- typedef typename Dune::LocalFiniteElementFunctionBase<LP1FE>::type FunctionBaseClass;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- AlienElementLocalBasisFunction<Element, LP1FE, FunctionBaseClass> f(*it, *it, coarseBaseSet.localBasis());
-
- for (int i=0; i<coarseBaseSet.localBasis().size(); i++) {
-
- f.setIndex(i);
- std::vector<double> values;
- fineBasis.getLocalFiniteElement(*it).localInterpolation().interpolate(f,values);
-
- for (size_t j=0; j<values.size(); j++) {
-
- if (values[j] > 0.001)
- {
- size_t globalFine = fineBasis.index(*it, j);
- size_t globalCoarse = p1Basis.index(*it, i);
-
- (*this->matrix_)[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<double,TransferMatrixBlock::rows>(values[j]);
- }
-
- }
-
- }
-
- }
-
- }
-
-};
-
-#endif
diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index bdd7f45c33d1e05c9d375e9ccd133ec3efaece5e..b1780ea9dd2d963ffbeb49be53d241800259d38d 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -12,14 +12,12 @@
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
// Using a monotone multigrid as the inner solver
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#include <dune/solvers/iterationsteps/mmgstep.hh>
#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
-#if defined THIRD_ORDER || defined SECOND_ORDER
-#include <dune/gfe/pktop1mgtransfer.hh>
-#endif
#include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include "maxnormtrustregion.hh"
@@ -189,10 +187,13 @@ setup(const GridType& grid,
#if defined THIRD_ORDER || defined SECOND_ORDER
{
+ typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
P1NodalBasis<typename GridType::LeafGridView,double> p1Basis(grid_->leafGridView());
- PKtoP1MGTransfer<CorrectionType>* topTransferOp = new PKtoP1MGTransfer<CorrectionType>;
- topTransferOp->setup(basis,p1Basis);
+ TransferOperatorType pkToP1TransferMatrix;
+ assembleBasisInterpolationMatrix<TransferOperatorType,
+ P1NodalBasis<typename GridType::LeafGridView,double>,
+ BasisType>(pkToP1TransferMatrix,p1Basis,basis);
// If we are on more than 1 processors, join all local transfer matrices on rank 0,
// and construct a single global transfer operator there.
@@ -202,7 +203,6 @@ setup(const GridType& grid,
typedef Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView, Dune::MCMGVertexLayout> LeafP1LocalMapper;
LeafP1LocalMapper leafP1LocalMapper(grid_->leafGridView());
- typedef typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType TransferOperatorType;
MatrixCommunicator<GlobalMapper,
typename GridType::LeafGridView,
typename GridType::LeafGridView,
@@ -211,9 +211,9 @@ setup(const GridType& grid,
LeafP1LocalMapper,
GlobalLeafP1Mapper> matrixComm(*globalMapper_, p1Index, grid_->leafGridView(), grid_->leafGridView(), localMapper, leafP1LocalMapper, 0);
- mmgStep->mgTransfer_.back() = new PKtoP1MGTransfer<CorrectionType>;
- Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(topTransferOp->getMatrix()));
- dynamic_cast<PKtoP1MGTransfer<CorrectionType>*>(mmgStep->mgTransfer_.back())->setMatrix(topTransferOperator);
+ mmgStep->mgTransfer_.back() = new TruncatedCompressedMGTransfer<CorrectionType>;
+ Dune::shared_ptr<TransferOperatorType> topTransferOperator = Dune::make_shared<TransferOperatorType>(matrixComm.reduceCopy(pkToP1TransferMatrix));
+ dynamic_cast<TruncatedCompressedMGTransfer<CorrectionType>*>(mmgStep->mgTransfer_.back())->setMatrix(topTransferOperator);
}
#endif