Commit 75558daf authored by Andreas Nüßing's avatar Andreas Nüßing

[Transfer] extract rhs

We introduce an interface for the transfer matrix right hand side and
modify the existing transfer matrix rhs (fitted and unfitted) to conform
to this inteface
parent 5b2e1178
......@@ -129,9 +129,9 @@ namespace duneuro
rightHandSideVector = 0.0;
TransferMatrixRHS<typename Traits::FunctionSpace::GFS> rhsAssembler(
solver_->functionSpace().getGFS());
rhsAssembler.assembleRightHandSide(reference.element, reference.localPosition,
electrode.element, electrode.localPosition,
rightHandSideVector);
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler.assembleRightHandSide(rightHandSideVector);
timer.stop();
dataTree.set("time_rhs_assembly", timer.lastElapsed());
timer.start();
......
......@@ -7,60 +7,74 @@
#include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
#include <duneuro/eeg/transfer_matrix_rhs_interface.hh>
namespace duneuro
{
template <class GFS>
class TransferMatrixRHS
: public TransferMatrixRHSInterface<
typename GFS::Traits::GridViewType,
Dune::PDELab::Backend::Vector<GFS, typename GFS::Traits::GridViewType::ctype>>
{
public:
using GV = typename GFS::Traits::GridViewType;
enum { dim = GV::dimension };
using Real = typename GV::ctype;
using Coordinate = Dune::FieldVector<Real, dim>;
using Element = typename GV::template Codim<0>::Entity;
using DOFVector = Dune::PDELab::Backend::Vector<GFS, Real>;
using BaseT = TransferMatrixRHSInterface<GV, Dune::PDELab::Backend::Vector<GFS, Real>>;
using LocalCoordinate = typename BaseT::LocalCoordinate;
using Element = typename BaseT::Element;
using VectorType = typename BaseT::VectorType;
using LFS = Dune::PDELab::LocalFunctionSpace<GFS>;
using Cache = Dune::PDELab::LFSIndexCache<LFS>;
using FESwitch = Dune::FiniteElementInterfaceSwitch<typename LFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
using RangeType = typename BasisSwitch::Range;
TransferMatrixRHS(const GFS& gfs) : gfs_(gfs)
TransferMatrixRHS(const GFS& gfs)
: gfs_(gfs)
, lfsReference_(gfs_)
, cacheReference_(lfsReference_)
, lfsElectrode_(gfs_)
, cacheElectrode_(lfsElectrode_)
{
}
void assembleRightHandSide(const Element& referenceElement, const Coordinate& referenceLocal,
const Element& electrodeElement, const Coordinate& electrodeLocal,
DOFVector& output)
virtual void bind(const Element& referenceElement, const LocalCoordinate& referenceLocal,
const Element& electrodeElement,
const LocalCoordinate& electrodeLocal) override
{
using FESwitch = Dune::FiniteElementInterfaceSwitch<typename LFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
using RangeType = typename BasisSwitch::Range;
// evaluate basis functions at reference position
LFS lfs(gfs_);
Cache cache(lfs);
lfs.bind(electrodeElement);
cache.update();
std::vector<RangeType> phi(lfs.size());
FESwitch::basis(lfs.finiteElement()).evaluateFunction(electrodeLocal, phi);
bind(lfsReference_, cacheReference_, phiReference_, referenceElement, referenceLocal);
bind(lfsElectrode_, cacheElectrode_, phiElectrode_, electrodeElement, electrodeLocal);
}
for (unsigned int i = 0; i < lfs.size(); ++i) {
output[cache.containerIndex(i)] = phi[i];
virtual void assembleRightHandSide(VectorType& output) const override
{
for (unsigned int i = 0; i < cacheElectrode_.size(); ++i) {
output[cacheElectrode_.containerIndex(i)] = phiElectrode_[i];
}
// evaluate basis functions at electrode
lfs.bind(referenceElement);
cache.update();
phi.resize(lfs.size());
FESwitch::basis(lfs.finiteElement()).evaluateFunction(referenceLocal, phi);
for (unsigned int i = 0; i < lfs.size(); ++i) {
output[cache.containerIndex(i)] -= phi[i];
for (unsigned int i = 0; i < cacheReference_.size(); ++i) {
output[cacheReference_.containerIndex(i)] -= phiReference_[i];
}
}
private:
const GFS& gfs_;
LFS lfsReference_;
Cache cacheReference_;
std::vector<RangeType> phiReference_;
LFS lfsElectrode_;
Cache cacheElectrode_;
std::vector<RangeType> phiElectrode_;
void bind(LFS& lfs, Cache& cache, std::vector<RangeType>& phi, const Element& element,
const LocalCoordinate& local)
{
lfs.bind(element);
cache.update();
phi.resize(lfs.size());
FESwitch::basis(lfs.finiteElement()).evaluateFunction(local, phi);
}
};
}
......
#ifndef DUNEURO_TRANSFER_MATRIX_RHS_INTERFACE_HH
#define DUNEURO_TRANSFER_MATRIX_RHS_INTERFACE_HH
#include <dune/common/fvector.hh>
namespace duneuro
{
template <class GV, class V>
struct TransferMatrixRHSInterface {
public:
using Element = typename GV::template Codim<0>::Entity;
using LocalCoordinate = Dune::FieldVector<typename GV::ctype, GV::dimension>;
using VectorType = V;
virtual void bind(const Element& referenceElement, const LocalCoordinate& referenceLocal,
const Element& electrodeElement, const LocalCoordinate& electrodeLocal) = 0;
virtual void assembleRightHandSide(VectorType& vector) const = 0;
virtual ~TransferMatrixRHSInterface(){};
};
}
#endif // DUNEURO_TRANSFER_MATRIX_RHS_INTERFACE_HH
......@@ -11,104 +11,104 @@
#include <dune/udg/pdelab/assembler/ulocalfunctionspace.hh>
#include <dune/udg/pdelab/subtriangulation.hh>
#include <duneuro/eeg/transfer_matrix_rhs_interface.hh>
namespace duneuro
{
template <class GFS, class ST>
class UnfittedTransferMatrixRHS
: public TransferMatrixRHSInterface<
typename GFS::Traits::GridViewType,
Dune::PDELab::Backend::Vector<GFS, typename GFS::Traits::GridViewType::ctype>>
{
public:
using GV = typename GFS::Traits::GridViewType;
enum { dim = GV::dimension };
using Real = typename GV::ctype;
using BaseT = TransferMatrixRHSInterface<GV, Dune::PDELab::Backend::Vector<GFS, Real>>;
using Coordinate = Dune::FieldVector<Real, dim>;
using Element = typename GV::template Codim<0>::Entity;
using DOFVector = Dune::PDELab::Backend::Vector<GFS, Real>;
using LocalCoordinate = typename BaseT::LocalCoordinate;
using Element = typename BaseT::Element;
using VectorType = typename BaseT::VectorType;
using UST = Dune::PDELab::UnfittedSubTriangulation<GV>;
using ULFS = Dune::PDELab::UnfittedLocalFunctionSpace<GFS>;
using UCache = Dune::PDELab::LFSIndexCache<ULFS>;
using ChildLFS = typename ULFS::template Child<0>::Type;
using FESwitch =
Dune::FiniteElementInterfaceSwitch<typename ChildLFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
using RangeType = typename BasisSwitch::Range;
UnfittedTransferMatrixRHS(const GFS& gfs, std::shared_ptr<ST> st, std::size_t child,
bool scaleToBBox)
: st_(st), gfs_(gfs), child_(child), scaleToBBox_(scaleToBBox)
: st_(st)
, gfs_(gfs)
, child_(child)
, scaleToBBox_(scaleToBBox)
, ulfsReference_(gfs_)
, cacheReference_(ulfsReference_)
, ulfsElectrode_(gfs_)
, cacheElectrode_(ulfsElectrode_)
{
}
void assembleRightHandSide(const Element& referenceElement, const Coordinate& referenceLocal,
const Element& electrodeElement, const Coordinate& electrodeLocal,
DOFVector& output)
virtual void bind(const Element& referenceElement, const LocalCoordinate& referenceLocal,
const Element& electrodeElement,
const LocalCoordinate& electrodeLocal) override
{
using ChildLFS = typename ULFS::template Child<0>::Type;
using FESwitch =
Dune::FiniteElementInterfaceSwitch<typename ChildLFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
using RangeType = typename BasisSwitch::Range;
UST ust(st_->gridView(), *st_);
bind(ulfsReference_, cacheReference_, phiReference_, referenceElement, referenceLocal);
bind(ulfsElectrode_, cacheElectrode_, phiElectrode_, electrodeElement, electrodeLocal);
}
ULFS ulfs(gfs_);
UCache ucache(ulfs);
virtual void assembleRightHandSide(VectorType& output) const override
{
const ChildLFS& childElectrode = ulfsElectrode_.child(child_);
for (unsigned int i = 0; i < childElectrode.size(); ++i) {
output[cacheElectrode_.containerIndex(childElectrode.localIndex(i))] = phiElectrode_[i];
}
const ChildLFS& childReference = ulfsReference_.child(child_);
for (unsigned int i = 0; i < childReference.size(); ++i) {
output[cacheReference_.containerIndex(childReference.localIndex(i))] -= phiReference_[i];
}
}
private:
std::shared_ptr<ST> st_;
const GFS& gfs_;
std::size_t child_;
bool scaleToBBox_;
ULFS ulfsReference_;
UCache cacheReference_;
std::vector<RangeType> phiReference_;
ULFS ulfsElectrode_;
UCache cacheElectrode_;
std::vector<RangeType> phiElectrode_;
void bind(ULFS& ulfs, UCache& ucache, std::vector<RangeType>& phi, const Element& element,
const LocalCoordinate& local)
{
UST ust(st_->gridView(), *st_);
ChildLFS& childLfs = ulfs.child(child_);
ust.create(electrodeElement);
ust.create(element);
if (ust.begin() == ust.end()) {
DUNE_THROW(Dune::Exception, "subtriangulation has no parts");
}
for (const auto& ep : ust) {
if (ep.domainIndex() != child_)
continue;
ulfs.bind(ep.subEntity(), true);
ucache.update();
FESwitch::basis(childLfs.finiteElement()).reset();
if (childLfs.size() == 0) {
DUNE_THROW(Dune::Exception, "child lfs has zero size");
}
std::vector<RangeType> phi(childLfs.size());
auto local =
scaleToBBox_ ?
ep.boundingBox().local(electrodeElement.geometry().global(electrodeLocal)) :
electrodeLocal;
FESwitch::basis(childLfs.finiteElement()).evaluateFunction(local, phi);
for (unsigned int i = 0; i < childLfs.size(); ++i) {
output[ucache.containerIndex(childLfs.localIndex(i))] = phi[i];
}
break;
}
ust.create(referenceElement);
for (const auto& ep : ust) {
if (ep.domainIndex() != child_)
continue;
ulfs.bind(ep.subEntity(), true);
ucache.update();
FESwitch::basis(childLfs.finiteElement()).reset();
assert(childLfs.size() > 0);
std::vector<RangeType> phi(childLfs.size());
auto local =
scaleToBBox_ ?
ep.boundingBox().local(referenceElement.geometry().global(referenceLocal)) :
referenceLocal;
FESwitch::basis(childLfs.finiteElement()).evaluateFunction(local, phi);
for (unsigned int i = 0; i < childLfs.size(); ++i) {
output[ucache.containerIndex(childLfs.localIndex(i))] -= phi[i];
}
phi.resize(childLfs.size());
auto felocal =
scaleToBBox_ ? ep.boundingBox().local(element.geometry().global(local)) : local;
FESwitch::basis(childLfs.finiteElement()).evaluateFunction(felocal, phi);
break;
}
}
private:
std::shared_ptr<ST> st_;
const GFS& gfs_;
std::size_t child_;
bool scaleToBBox_;
};
}
......
......@@ -124,16 +124,16 @@ namespace duneuro
rhsAssembler(solver_->functionSpace().getGFS(), subTriangulation_,
config.get<std::size_t>("compartment"), scaleToBBox_);
#if HAVE_TBB
{
{ // mutex, as the finite element map is not thread safe
tbb::mutex::scoped_lock lock(solver_->functionSpaceMutex());
rhsAssembler.assembleRightHandSide(reference.element, reference.localPosition,
electrode.element, electrode.localPosition,
rightHandSideVector);
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
}
rhsAssembler.assembleRightHandSide(rightHandSideVector);
#else
rhsAssembler.assembleRightHandSide(reference.element, reference.localPosition,
electrode.element, electrode.localPosition,
rightHandSideVector);
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler.assembleRightHandSide(rightHandSideVector);
#endif
timer.stop();
dataTree.set("time_rhs_assembly", timer.lastElapsed());
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment