diff --git a/dune.module b/dune.module
index e1179a02c978bbc88772b4ba37263371a8ea2736..f42cc5b4f7cefa8b219d6bf4e9b07bb7c06f0383 100644
--- a/dune.module
+++ b/dune.module
@@ -7,5 +7,5 @@ Module: dune-solvers
Version: 2.6-git
Maintainer: oliver.sander@tu-dresden.de
#depending on
-Depends: dune-common dune-grid dune-istl dune-localfunctions dune-matrix-vector
+Depends: dune-common dune-grid dune-istl dune-localfunctions dune-matrix-vector dune-functions
Suggests: dune-alugrid
diff --git a/dune/solvers/transferoperators/compressedmultigridtransfer.hh b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
index 1cdd9ebd8891da21db9062fb4ec3fc30310dde3c..44088600daab085b022ca51bb067714a8301f520 100644
--- a/dune/solvers/transferoperators/compressedmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/compressedmultigridtransfer.hh
@@ -64,11 +64,12 @@ public:
* \param grid The grid
* \param cL The coarse grid level
* \param fL The fine grid level
+ * \param order The order of the lagrange basis
*/
- template <class GridType>
- void setup(const GridType& grid, int cL, int fL)
+ template <class GridType, int k = 1>
+ void setup(const GridType& grid, int cL, int fL, std::integral_constant<int,k> order = {})
{
- GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL);
+ GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL, order);
}
@@ -121,7 +122,7 @@ public:
}
- /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
+ /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
*
* Set occupation of Galerkin restricted coarse matrix. Call this one before
* galerkinRestrict to ensure all non-zeroes are present
@@ -195,11 +196,12 @@ public:
* \param grid The grid
* \param cL The coarse grid level
* \param fL The fine grid level
+ * \param order The order of the lagrange basis
*/
- template <class GridType>
- void setup(const GridType& grid, int cL, int fL)
+ template <class GridType, int k = 1>
+ void setup(const GridType& grid, int cL, int fL, std::integral_constant<int,k> order = {})
{
- GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL);
+ GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL, order);
}
@@ -271,7 +273,7 @@ public:
}
- /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
+ /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
*
* Set occupation of Galerkin restricted coarse matrix. Call this one before
* galerkinRestrict to ensure all non-zeroes are present
diff --git a/dune/solvers/transferoperators/densemultigridtransfer.hh b/dune/solvers/transferoperators/densemultigridtransfer.hh
index 994a0036eee4facee3fc5ed473439e60d22f2e43..cc85a6c015b1cd40673308e03140acb0741276d3 100644
--- a/dune/solvers/transferoperators/densemultigridtransfer.hh
+++ b/dune/solvers/transferoperators/densemultigridtransfer.hh
@@ -46,11 +46,12 @@ public:
* \param grid The grid
* \param cL The coarse grid level
* \param fL The fine grid level
+ * \param order The order of the lagrange basis
*/
- template <class GridType>
- void setup(const GridType& grid, int cL, int fL)
+ template <class GridType, int k = 1>
+ void setup(const GridType& grid, int cL, int fL, std::integral_constant<int,k> order = {})
{
- GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(matrix_, grid, cL, fL);
+ GenericMultigridTransfer::setup<TransferOperatorType, GridType, field_type>(*matrix_, grid, cL, fL, order);
}
@@ -103,7 +104,7 @@ public:
}
- /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
+ /** \brief Set Occupation of Galerkin restricted coarse stiffness matrix
*
* Set occupation of Galerkin restricted coarse matrix. Call this one before
* galerkinRestrict to ensure all non-zeroes are present
diff --git a/dune/solvers/transferoperators/genericmultigridtransfer.hh b/dune/solvers/transferoperators/genericmultigridtransfer.hh
index 64352907d7f8edab35a3507727e3b5ee78be1a8e..07031d42f20c5a8ff5f962c2f189a9a613b8d8aa 100644
--- a/dune/solvers/transferoperators/genericmultigridtransfer.hh
+++ b/dune/solvers/transferoperators/genericmultigridtransfer.hh
@@ -8,6 +8,8 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/bitsetvector.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
@@ -23,7 +25,7 @@
* as a BCRSMatrix. Therefore, the template parameter DiscFuncType
* has to comply with the ISTL requirements.
- * \todo Currently only works for first-order Lagrangian elements!
+ * \todo Currently only works for Lagrangian elements!
*/
//template<class DiscFuncType, class MatrixBlock, class TransferOperatorType>
class GenericMultigridTransfer {
@@ -63,27 +65,46 @@ public:
}
- template<class TransferOperatorType, class GridType, class field_type>
- static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL)
+ /** \brief Set up transfer operator between two directly connected grid levels of the same grid
+
+ The 'fine' grid level (fL) is expected to be the first refinement of the 'coarse' one (cL).
+
+ \param matrix The matrix object to assemble into
+ \param grid The hierarchial grid
+ \param cL Coarse level
+ \param fL Fine level ( = cL+1 )
+ \tparam k compile time parameter: order of the larange basis
+ */
+ template<class TransferOperatorType, class GridType, class field_type, int k = 1>
+ static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL, std::integral_constant<int,k> = {})
{
- typedef typename TransferOperatorType::block_type TransferMatrixBlock;
- typedef typename GridType::ctype ctype;
+ using TransferMatrixBlock = typename TransferOperatorType::block_type;
+ using ctype = typename GridType::ctype;
if (fL != cL+1)
DUNE_THROW(Dune::Exception, "The two function spaces don't belong to consecutive levels!");
const int dim = GridType::dimension;
- const typename GridType::Traits::LevelIndexSet& coarseIndexSet = grid.levelIndexSet(cL);
- const typename GridType::Traits::LevelIndexSet& fineIndexSet = grid.levelIndexSet(fL);
+ const auto& coarseIndexSet = grid.levelIndexSet(cL);
+ const auto& fineIndexSet = grid.levelIndexSet(fL);
- int rows = grid.size(fL, dim);
- int cols = grid.size(cL, dim);
+ const auto& coarseLevelGridView = grid.levelGridView(cL);
+ const auto& fineLevelGridView = grid.levelGridView(fL);
+
+ // construct a global basis to determine the size of the transfer matrix
+ Dune::Functions::PQkNodalBasis< typename GridType::LevelGridView, k> coarseBasis(coarseLevelGridView);
+ Dune::Functions::PQkNodalBasis< typename GridType::LevelGridView, k> fineBasis(fineLevelGridView);
+
+ const int rows = fineBasis.size();
+ const int cols = coarseBasis.size();
// A factory for the shape functions
- typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
- typedef typename P1FECache::FiniteElementType FEType;
- P1FECache p1FECache;
+ using PkFECache = typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, k>;
+ PkFECache pkFECache;
+
+ // minimal value to be determined in the support
+ double tol = 0.001;
matrix.setSize(rows,cols);
matrix.setBuildMode(TransferOperatorType::random);
@@ -95,11 +116,8 @@ public:
for (const auto& cIt : elements(grid.levelGridView(cL))) {
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
// Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(cIt.type());
-
+ const auto& coarseBaseSet = pkFECache.get(cIt.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
@@ -111,18 +129,15 @@ public:
continue;
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
+ const auto& fGeometryInFather = fIt.geometryInFather();
// Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());
-
+ const auto& fineBaseSet = pkFECache.get(fIt.type());
const size_t numFineBaseFct = fineBaseSet.localBasis().size();
-
for (size_t j=0; j<numFineBaseFct; j++)
{
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
+ const auto& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
@@ -134,17 +149,15 @@ public:
for (size_t i=0; i<numCoarseBaseFct; i++)
{
- if (values[i] > 0.001)
+ if (values[i] > tol)
{
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
+ const auto& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(cIt, iLocalKey.subEntity(), iLocalKey.codim());
indices.add(globalFine, globalCoarse);
}
}
}
-
}
-
}
indices.exportIdx(matrix);
@@ -156,12 +169,9 @@ public:
for (const auto& cIt : elements(grid.levelGridView(cL))) {
// Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(cIt.type());
-
+ const auto& coarseBaseSet = pkFECache.get(cIt.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
- typedef typename GridType::template Codim<0>::Entity EntityType;
-
// preallocate vector for function evaluations
std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
@@ -171,18 +181,15 @@ public:
continue;
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
-
- const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
+ const auto& fGeometryInFather = fIt.geometryInFather();
// Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());
-
+ const auto& fineBaseSet = pkFECache.get(fIt.type());
const size_t numFineBaseFct = fineBaseSet.localBasis().size();
for (size_t j=0; j<numFineBaseFct; j++)
{
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
-
+ const auto& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim());
@@ -193,11 +200,10 @@ public:
for (size_t i=0; i<numCoarseBaseFct; i++)
{
- if (values[i] > 0.001)
+ if (values[i] > tol)
{
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
+ const auto& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(cIt, iLocalKey.subEntity(), iLocalKey.codim());
-
matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
}
}
@@ -212,25 +218,34 @@ public:
The 'fine' grid view gvFine is expected to be a refinement of the 'coarse' one gvCoarse.
\param matrix The matrix object to assemble into
+ \param gvCoarse The coarse grid
+ \param gvFine The fine grid level
+ \tparam k compile time parameter: order of the larange basis
*/
- template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type>
- static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine)
+ template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type, int k = 1>
+ static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine, std::integral_constant<int,k> = {})
{
- typedef typename TransferOperatorType::block_type TransferMatrixBlock;
- typedef typename GridViewCoarse::ctype ctype;
+ using TransferMatrixBlock = typename TransferOperatorType::block_type;
+ using ctype = typename GridViewCoarse::ctype;
const int dim = GridViewCoarse::dimension;
- const typename GridViewCoarse::IndexSet& coarseIndexSet = gvCoarse.indexSet();
- const typename GridViewFine::IndexSet& fineIndexSet = gvFine.indexSet();
+ const auto& coarseIndexSet = gvCoarse.indexSet();
+ const auto& fineIndexSet = gvFine.indexSet();
- int rows = gvFine.size(dim);
- int cols = gvCoarse.size(dim);
+ // construct a global basis to determine the size of the transfer matrix
+ Dune::Functions::PQkNodalBasis< GridViewCoarse, k> coarseBasis(gvCoarse);
+ Dune::Functions::PQkNodalBasis< GridViewFine , k> fineBasis(gvFine);
+
+ const int rows = fineBasis.size();
+ const int cols = coarseBasis.size();
// A factory for the shape functions
- typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
- typedef typename P1FECache::FiniteElementType FEType;
- P1FECache p1FECache;
+ using PkFECache = typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, k>;
+ PkFECache pkFECache;
+
+ // minimal value to be determined in the support
+ double tol = 0.001;
matrix.setSize(rows,cols);
matrix.setBuildMode(TransferOperatorType::random);
@@ -246,8 +261,7 @@ public:
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
// Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());;
-
+ const auto& fineBaseSet = pkFECache.get(fIt.type());;
const size_t numFineBaseFct = fineBaseSet.localBasis().size();
std::vector<Dune::FieldVector<ctype,dim> > local(numFineBaseFct);
@@ -264,7 +278,6 @@ public:
while (not gvCoarse.contains(ancestor) && ancestor.level() != 0 ) {
auto geometryInFather = ancestor.geometryInFather();
-
for (size_t i=0; i<numFineBaseFct; i++)
local[i] = geometryInFather.global(local[i]);
ancestor = ancestor.father();
@@ -276,8 +289,7 @@ public:
for (size_t j=0; j<numFineBaseFct; j++)
{
// Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
-
+ const auto& coarseBaseSet = pkFECache.get(ancestor.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
@@ -291,15 +303,14 @@ public:
for (size_t i=0; i<numCoarseBaseFct; i++)
{
- if (values[i] > 0.001)
+ if (values[i] > tol)
{
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
+ const auto& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
indices.add(globalFine, globalCoarse);
}
}
}
-
}
indices.exportIdx(matrix);
@@ -309,19 +320,17 @@ public:
// /////////////////////////////////////////////
for (const auto& fIt : elements(gvFine)) {
-
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
// Get local finite element
- const FEType& fineBaseSet = p1FECache.get(fIt.type());;
-
+ const auto& fineBaseSet = pkFECache.get(fIt.type());;
const size_t numFineBaseFct = fineBaseSet.localBasis().size();
std::vector<Dune::FieldVector<ctype,dim> > local(numFineBaseFct);
for (size_t i=0; i<numFineBaseFct; i++)
{
- const Dune::LocalKey& iLocalKey = fineBaseSet.localCoefficients().localKey(i);
+ const auto& iLocalKey = fineBaseSet.localCoefficients().localKey(i);
local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
}
@@ -343,8 +352,7 @@ public:
for (size_t j=0; j<numFineBaseFct; j++)
{
// Get local finite element
- const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
-
+ const auto& coarseBaseSet = pkFECache.get(ancestor.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
@@ -353,16 +361,15 @@ public:
// Evaluate coarse grid base functions
coarseBaseSet.localBasis().evaluateFunction(local[j], values);
- const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
+ const auto& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
for (size_t i=0; i<numCoarseBaseFct; i++)
{
- if (values[i] > 0.001)
+ if (values[i] > tol)
{
- const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
+ const auto& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
-
matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
}
}