From ffea6cdf71e3270f892ae61e96afdb8b2c799cf6 Mon Sep 17 00:00:00 2001
From: Patrick Jaap <patrick.jaap@tu-dresden.de>
Date: Tue, 9 Jan 2018 15:36:27 +0100
Subject: [PATCH] Implement 'setup' for all lagrange element transfer operators
Indead of limiting this to P1 elements, we can now build transfer matrices for
lagrange elements of any order. This is achived by a compile time parameter "k".
For back compatibility, it has a default value k = 1.
To determine the size of the matrix, we need a global basis from dune-functions.
---
dune.module | 2 +-
.../compressedmultigridtransfer.hh | 18 ++-
.../densemultigridtransfer.hh | 9 +-
.../genericmultigridtransfer.hh | 143 +++++++++---------
4 files changed, 91 insertions(+), 81 deletions(-)
diff --git a/dune.module b/dune.module
index e1179a0..f42cc5b 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 1cdd9eb..4408860 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 994a003..cc85a6c 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 6435290..07031d4 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]);
}
}
--
GitLab