Skip to content
Snippets Groups Projects
Commit 547b5acd authored by Patrick Jaap's avatar Patrick Jaap
Browse files

WIP: use generic method in transfer operators

parent b3a40c8d
Branches
Tags
No related merge requests found
Pipeline #28205 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
dune/solvers/transferoperators/genericmultigridtransfer.hh
+ 15
− 292
View file @ 547b5acd
...@@ -8,6 +8,9 @@ ...@@ -8,6 +8,9 @@
#include <dune/common/fmatrix.hh> #include <dune/common/fmatrix.hh>
#include <dune/common/bitsetvector.hh> #include <dune/common/bitsetvector.hh>
#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/geometry/referenceelements.hh> #include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh> #include <dune/localfunctions/lagrange/pqkfactory.hh>
...@@ -64,146 +67,12 @@ public: ...@@ -64,146 +67,12 @@ public:
template<class TransferOperatorType, class GridType, class field_type> template<class TransferOperatorType, class GridType, class field_type>
static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL) static void setup(TransferOperatorType& matrix, const GridType& grid, int cL, int fL, int order = 1)
{ {
typedef typename TransferOperatorType::block_type TransferMatrixBlock; // call function overload with two gridViews
typedef typename GridType::ctype ctype; auto fineGridView = grid.levelGridView(fL);
auto coarseGridView = grid.levelGridView(cL);
if (fL != cL+1) setup(matrix, coarseGridView, fineGridView, order);
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);
int rows = grid.size(fL, dim);
int cols = grid.size(cL, dim);
// A factory for the shape functions
typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
typedef typename P1FECache::FiniteElementType FEType;
P1FECache p1FECache;
matrix.setSize(rows,cols);
matrix.setBuildMode(TransferOperatorType::random);
// ///////////////////////////////////////////
// Determine the indices present in the matrix
// /////////////////////////////////////////////////
Dune::MatrixIndexSet indices(rows, cols);
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 size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
for (const auto& fIt : descendantElements(cIt,fL)) {
if (fIt.level()==cIt.level())
continue;
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
// Get local finite element
const FEType& fineBaseSet = p1FECache.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);
int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim());
Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
// Evaluate coarse grid base functions
coarseBaseSet.localBasis().evaluateFunction(local, values);
for (size_t i=0; i<numCoarseBaseFct; i++)
{
if (values[i] > 0.001)
{
const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(cIt, iLocalKey.subEntity(), iLocalKey.codim());
indices.add(globalFine, globalCoarse);
}
}
}
}
}
indices.exportIdx(matrix);
// /////////////////////////////////////////////
// Compute the matrix
// /////////////////////////////////////////////
for (const auto& cIt : elements(grid.levelGridView(cL))) {
// Get local finite element
const FEType& coarseBaseSet = p1FECache.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);
for (const auto& fIt : descendantElements(cIt,fL)) {
if (fIt.level()==cIt.level())
continue;
const auto& fineRefElement = Dune::ReferenceElements<ctype, dim>::general(fIt.type());
const typename EntityType::LocalGeometry& fGeometryInFather = fIt.geometryInFather();
// Get local finite element
const FEType& fineBaseSet = p1FECache.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);
int globalFine = fineIndexSet.subIndex(fIt, jLocalKey.subEntity(), jLocalKey.codim());
Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), jLocalKey.codim());
Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
// Evaluate coarse grid base functions
coarseBaseSet.localBasis().evaluateFunction(local, values);
for (size_t i=0; i<numCoarseBaseFct; i++)
{
if (values[i] > 0.001)
{
const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(cIt, iLocalKey.subEntity(), iLocalKey.codim());
matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
}
}
}
}
}
} }
...@@ -213,161 +82,15 @@ public: ...@@ -213,161 +82,15 @@ public:
\param matrix The matrix object to assemble into \param matrix The matrix object to assemble into
*/ */
template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type> template<class TransferOperatorType, class GridViewCoarse, class GridViewFine, class field_type = double>
static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine) static void setup(TransferOperatorType& matrix, const GridViewCoarse& gvCoarse, const GridViewFine& gvFine, int order = 1)
{ {
typedef typename TransferOperatorType::block_type TransferMatrixBlock; // create Lagrange bases of given order for fine and coarse grid level
typedef typename GridViewCoarse::ctype ctype; auto fineBasis = Dune::Functions::BasisFactory::makeBasis(gvFine, Dune::Functions::BasisFactory::lagrange(order));
auto coarseBasis = Dune::Functions::BasisFactory::makeBasis(gvCoarse,Dune::Functions::BasisFactory::lagrange(order));
const int dim = GridViewCoarse::dimension;
const typename GridViewCoarse::IndexSet& coarseIndexSet = gvCoarse.indexSet();
const typename GridViewFine::IndexSet& fineIndexSet = gvFine.indexSet();
int rows = gvFine.size(dim);
int cols = gvCoarse.size(dim);
// A factory for the shape functions // call generic method to assemble the transfer matrix
typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache; assembleGlobalBasisTransferMatrix(matrix,coarseBasis,fineBasis);
typedef typename P1FECache::FiniteElementType FEType;
P1FECache p1FECache;
matrix.setSize(rows,cols);
matrix.setBuildMode(TransferOperatorType::random);
// ///////////////////////////////////////////
// Determine the indices present in the matrix
// /////////////////////////////////////////////////
Dune::MatrixIndexSet indices(rows, cols);
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 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);
local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
}
// Get ancestor element in the coarse grid view, and the local position there.
auto ancestor = fIt;
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();
}
assert(gvCoarse.contains(ancestor));
for (size_t j=0; j<numFineBaseFct; j++)
{
// Get local finite element
const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
// Evaluate coarse grid base functions
coarseBaseSet.localBasis().evaluateFunction(local[j], values);
const Dune::LocalKey& 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)
{
const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
indices.add(globalFine, globalCoarse);
}
}
}
}
indices.exportIdx(matrix);
// /////////////////////////////////////////////
// Compute the matrix
// /////////////////////////////////////////////
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 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);
local[i] = fineRefElement.position(iLocalKey.subEntity(), iLocalKey.codim());
}
// Get ancestor element in the coarse grid view, and the local position there.
auto ancestor = fIt;
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();
}
assert(gvCoarse.contains(ancestor));
for (size_t j=0; j<numFineBaseFct; j++)
{
// Get local finite element
const FEType& coarseBaseSet = p1FECache.get(ancestor.type());
const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
// preallocate vector for function evaluations
std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
// Evaluate coarse grid base functions
coarseBaseSet.localBasis().evaluateFunction(local[j], values);
const Dune::LocalKey& 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)
{
const Dune::LocalKey& iLocalKey = coarseBaseSet.localCoefficients().localKey(i);
int globalCoarse = coarseIndexSet.subIndex(ancestor, iLocalKey.subEntity(), iLocalKey.codim());
matrix[globalFine][globalCoarse] = Dune::ScaledIdentityMatrix<ctype,TransferMatrixBlock::rows>(values[i]);
}
}
}
}
} }
//! Multiply the vector f from the right to the prolongation matrix //! Multiply the vector f from the right to the prolongation matrix
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment