From c2a078f165320c18042df8ed9dd28130dfa62d24 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sat, 27 Dec 2014 16:05:04 +0000
Subject: [PATCH] Replace PkToP1MGTransfer with
BasisInterpolationMatrixAssembler from dune-fufem
There is no need to have our own implementation for this, we had it purely for
historical reasons.
Also, this change fixes a subtle bug: Our own implementation used AlienElementLocalBasis,
which returns zero if the LocalBasis is evaluated outside of the reference element.
However, this happened sometimes due to rounding error on grids with very small quadrilateral
elements.
[[Imported from SVN: r9991]]
---
dune/gfe/Makefile.am | 1 -
dune/gfe/pktop1mgtransfer.hh | 150 ---------------------------------
dune/gfe/riemanniantrsolver.cc | 18 ++--
3 files changed, 9 insertions(+), 160 deletions(-)
delete mode 100644 dune/gfe/pktop1mgtransfer.hh
diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index 2d9e4802..0efff862 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 0ea4f118..00000000
--- 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 bdd7f45c..b1780ea9 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
--
GitLab