Select Git revision
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
dunefunctionsoperatorassembler.hh 9.86 KiB
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include "dune/fufem/functionspacebases/functionspacebasis.hh"
namespace Dune {
namespace Fufem {
//! Generic global assembler for operators on a gridview
template <class TrialBasis, class AnsatzBasis>
class DuneFunctionsOperatorAssembler
{
using GridView = typename TrialBasis::GridView;
public:
//! create assembler for grid
DuneFunctionsOperatorAssembler(const TrialBasis& tBasis, const AnsatzBasis& aBasis) :
trialBasis_(tBasis),
ansatzBasis_(aBasis)
{}
template <class PatternBuilder>
void assembleBulkPattern(PatternBuilder& patternBuilder) const
{
patternBuilder.resize(trialBasis_, ansatzBasis_);
auto trialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView();
for (const auto& element : elements(trialBasis_.gridView()))
{
trialLocalView.bind(element);
ansatzLocalView.bind(element);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t i=0; i<trialLocalView.size(); ++i)
{
auto row = trialLocalView.index(i);
for (size_t j=0; j<ansatzLocalView.size(); ++j)
{
auto col = ansatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
}
}
template <class PatternBuilder>
void assembleSkeletonPattern(PatternBuilder& patternBuilder) const
{
patternBuilder.resize(trialBasis_, ansatzBasis_);
auto insideTrialLocalView = trialBasis_.localView();
auto insideAnsatzLocalView = ansatzBasis_.localView();
auto outsideAnsatzLocalView = ansatzBasis_.localView();
for (const auto& element : elements(trialBasis_.gridView()))
{
insideTrialLocalView.bind(element);
insideAnsatzLocalView.bind(element);
/* Coupling on the same element */
for (size_t i = 0; i < insideTrialLocalView.size(); i++) {
auto row = insideTrialLocalView.index(i);
for (size_t j = 0; j < insideAnsatzLocalView.size(); j++) {
auto col = insideAnsatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
for (const auto& is : intersections(trialBasis_.gridView(), element))
{
/* Elements on the boundary have been treated above, iterate along the inner edges */
if (is.neighbor())
{
// Current outside element
const auto& outsideElement = is.outside();
// Bind the outer parts to the outer element
outsideAnsatzLocalView.bind(outsideElement);
// We assume that all basis functions of the inner element couple with all basis functions from the outer one
for (size_t i=0; i<insideTrialLocalView.size(); ++i)
{
auto row = insideTrialLocalView.index(i);
for (size_t j=0; j<outsideAnsatzLocalView.size(); ++j)
{
auto col = outsideAnsatzLocalView.index(j);
patternBuilder.insertEntry(row,col);
}
}
}
}
}
}
template <class MatrixBackend, class LocalAssembler>
void assembleBulkEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
{
auto trialLocalView = trialBasis_.localView();
auto ansatzLocalView = ansatzBasis_.localView();
using Field = std::decay_t<decltype(matrixBackend(trialLocalView.index(0), ansatzLocalView.index(0)))>;
using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
auto localMatrix = LocalMatrix();
for (const auto& element : elements(trialBasis_.gridView()))
{
trialLocalView.bind(element);
ansatzLocalView.bind(element);
localMatrix.setSize(trialLocalView.size(), ansatzLocalView.size());
localAssembler(element, localMatrix, trialLocalView, ansatzLocalView);
// Add element stiffness matrix onto the global stiffness matrix
for (size_t localRow=0; localRow<trialLocalView.size(); ++localRow)
{
auto row = trialLocalView.index(localRow);
for (size_t localCol=0; localCol<ansatzLocalView.size(); ++localCol)
{
auto col = ansatzLocalView.index(localCol);
matrixBackend(row,col) += localMatrix[localRow][localCol];
}
}
}
}
/* This variant of the IntersectionOperatorAssembler that works
* with dune-functions bases.
*
* Currently not supported are constrained bases and lumping.
*
* Note that this was written (and tested) having discontinuous Galerkin bases in mind.
* There may be cases where things do not work as expected for standard (continuous) Galerkin spaces.
*/
template <class MatrixBackend, class LocalAssembler, class LocalBoundaryAssembler>
void assembleSkeletonEntries(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler, LocalBoundaryAssembler&& localBoundaryAssembler) const
{
Dune::MultipleCodimMultipleGeomTypeMapper<GridView> faceMapper(trialBasis_.gridView(), mcmgElementLayout());
auto insideTrialLocalView = trialBasis_.localView();
auto insideAnsatzLocalView = ansatzBasis_.localView();
auto outsideTrialLocalView = trialBasis_.localView();
auto outsideAnsatzLocalView = ansatzBasis_.localView();
using Field = std::decay_t<decltype(matrixBackend(insideTrialLocalView.index(0), insideAnsatzLocalView.index(0)))>;
using LocalMatrix = Dune::Matrix<Dune::FieldMatrix<Field,1,1>>;
using MatrixContainer = Dune::Matrix<LocalMatrix>;
auto matrixContrainer = MatrixContainer(2,2);
for (const auto& element : elements(trialBasis_.gridView()))
{
insideTrialLocalView.bind(element);
insideAnsatzLocalView.bind(element);
// Resize
matrixContrainer[0][0].setSize(insideTrialLocalView.size(), insideAnsatzLocalView.size());
for (const auto& is : intersections(trialBasis_.gridView(), element))
{
/* Assemble (depending on whether we have a boundary edge or not) */
if (is.boundary())
{
auto& localMatrix = matrixContrainer[0][0];
localBoundaryAssembler(is, localMatrix, insideTrialLocalView, insideAnsatzLocalView);
for (size_t i=0; i< insideTrialLocalView.size(); i++) {
for (size_t j=0; j< insideAnsatzLocalView.size(); j++)
matrixBackend(insideTrialLocalView.index(i), insideAnsatzLocalView.index(j))+=localMatrix[i][j];
}
}
else if (is.neighbor())
{
/* Only handle every intersection once; we choose the case where the inner element has the lower index
*
* This should also work with grids where elements can have different geometries. TODO: Actually test this somewhere!*/
const auto& outsideElement = is.outside();
if (faceMapper.index(element) > faceMapper.index(outsideElement)) continue;
// Bind the outer parts to the outer element
outsideTrialLocalView.bind(outsideElement);
outsideAnsatzLocalView.bind(outsideElement);
matrixContrainer[0][1].setSize(insideTrialLocalView.size(), outsideAnsatzLocalView.size());
matrixContrainer[1][0].setSize(outsideTrialLocalView.size(), insideAnsatzLocalView.size());
matrixContrainer[1][1].setSize(outsideTrialLocalView.size(), outsideAnsatzLocalView.size());
localAssembler(is, matrixContrainer, insideTrialLocalView, insideAnsatzLocalView, outsideTrialLocalView, outsideAnsatzLocalView);
for (size_t i=0; i < insideTrialLocalView.size(); i++)
{
auto row = insideTrialLocalView.index(i);
// inside x inside
for (size_t j=0; j < insideTrialLocalView.size(); j++)
{
auto col = insideTrialLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[0][0][i][j];
}
// inside x outside
for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
{
auto col = outsideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[0][1][i][j];
}
}
for (size_t i=0; i < outsideTrialLocalView.size(); i++)
{
auto row = outsideTrialLocalView.index(i);
// outside x inside
for (size_t j=0; j < insideAnsatzLocalView.size(); j++)
{
auto col = insideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[1][0][i][j];
}
// outside x outside
for (size_t j=0; j < outsideAnsatzLocalView.size(); j++)
{
auto col = outsideAnsatzLocalView.index(j);
matrixBackend(row, col) += matrixContrainer[1][1][i][j];
}
}
}
}
}
}
template <class MatrixBackend, class LocalAssembler>
void assembleBulk(MatrixBackend&& matrixBackend, LocalAssembler&& localAssembler) const
{
auto patternBuilder = matrixBackend.patternBuilder();
assembleBulkPattern(patternBuilder);
patternBuilder.setupMatrix();
assembleBulkEntries(matrixBackend, std::forward<LocalAssembler>(localAssembler));
}
protected:
const TrialBasis& trialBasis_;
const AnsatzBasis& ansatzBasis_;
};
/**
* \brief Create DuneFunctionsOperatorAssembler */
template <class TrialBasis, class AnsatzBasis>
auto duneFunctionsOperatorAssembler(const TrialBasis& trialBasis, const AnsatzBasis& ansatzBasis)
{
return DuneFunctionsOperatorAssembler<TrialBasis, AnsatzBasis>(trialBasis, ansatzBasis);
}
} // namespace Fufem
} // namespace Dune
#endif // DUNE_FUFEM_ASSEMBLERS_DUNEFUNCTIONSOPERATORASSEMBLER_HH