Commit cccede23 authored by Andreas Nüßing's avatar Andreas Nüßing

[Transfer] create rhs factories

Instead of directly creating the transfer matrix rhs, we use a factory
for its constrution. This is basicaly a preparation for the unified
transfer matrix solver, als well as for different kinds of transfer
matrix rhs (such as cem)
parent 75558daf
......@@ -138,6 +138,11 @@ namespace duneuro
return search_;
}
bool scaleToBBox() const
{
return false;
}
#if HAVE_TBB
tbb::mutex& functionSpaceMutex()
{
......
......@@ -134,6 +134,11 @@ namespace duneuro
return search_;
}
bool scaleToBBox() const
{
return true;
}
#if HAVE_TBB
tbb::mutex& functionSpaceMutex()
{
......
#ifndef DUNEURO_FITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
#define DUNEURO_FITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
#include <dune/common/parametertree.hh>
#include <duneuro/eeg/transfer_matrix_rhs.hh>
#include <duneuro/eeg/transfer_matrix_rhs_interface.hh>
namespace duneuro
{
struct FittedTransferMatrixRHSFactory {
template <class Vector, class Solver>
static std::unique_ptr<TransferMatrixRHSInterface<
typename Solver::Traits::VolumeConductor::GridView, Vector>>
create(const Solver& solver, const Dune::ParameterTree& config)
{
auto type = config.get<std::string>("type", "point");
if (type == "point") {
return std::make_unique<TransferMatrixRHS<typename Solver::Traits::FunctionSpace::GFS,
Vector>>(solver.functionSpace().getGFS());
} else {
DUNE_THROW(Dune::Exception, "unknown transfer matrix type \"" << type << "\"");
}
}
};
}
#endif // DUNEURO_FITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
......@@ -26,7 +26,7 @@ namespace duneuro
using ProjectedPosition = ProjectedElectrode<typename VolumeConductor::GridView>;
};
template <class S>
template <class S, class RHSFactory>
class FittedTransferMatrixSolver
{
public:
......@@ -127,11 +127,11 @@ namespace duneuro
Dune::Timer timer;
// assemble right hand side
rightHandSideVector = 0.0;
TransferMatrixRHS<typename Traits::FunctionSpace::GFS> rhsAssembler(
solver_->functionSpace().getGFS());
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler.assembleRightHandSide(rightHandSideVector);
auto rhsAssembler =
RHSFactory::template create<typename Traits::RangeDOFVector>(*solver_, config);
rhsAssembler->bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler->assembleRightHandSide(rightHandSideVector);
timer.stop();
dataTree.set("time_rhs_assembly", timer.lastElapsed());
timer.start();
......
......@@ -11,17 +11,14 @@
namespace duneuro
{
template <class GFS>
class TransferMatrixRHS
: public TransferMatrixRHSInterface<
typename GFS::Traits::GridViewType,
Dune::PDELab::Backend::Vector<GFS, typename GFS::Traits::GridViewType::ctype>>
template <class GFS, class V>
class TransferMatrixRHS : public TransferMatrixRHSInterface<typename GFS::Traits::GridViewType, V>
{
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 BaseT = TransferMatrixRHSInterface<GV, V>;
using LocalCoordinate = typename BaseT::LocalCoordinate;
using Element = typename BaseT::Element;
using VectorType = typename BaseT::VectorType;
......
......@@ -15,17 +15,15 @@
namespace duneuro
{
template <class GFS, class ST>
template <class GFS, class ST, class V>
class UnfittedTransferMatrixRHS
: public TransferMatrixRHSInterface<
typename GFS::Traits::GridViewType,
Dune::PDELab::Backend::Vector<GFS, typename GFS::Traits::GridViewType::ctype>>
: public TransferMatrixRHSInterface<typename GFS::Traits::GridViewType, V>
{
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 BaseT = TransferMatrixRHSInterface<GV, V>;
using Coordinate = Dune::FieldVector<Real, dim>;
using LocalCoordinate = typename BaseT::LocalCoordinate;
using Element = typename BaseT::Element;
......@@ -39,7 +37,7 @@ namespace duneuro
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,
UnfittedTransferMatrixRHS(const GFS& gfs, std::shared_ptr<const ST> st, std::size_t child,
bool scaleToBBox)
: st_(st)
, gfs_(gfs)
......@@ -73,7 +71,7 @@ namespace duneuro
}
private:
std::shared_ptr<ST> st_;
std::shared_ptr<const ST> st_;
const GFS& gfs_;
std::size_t child_;
bool scaleToBBox_;
......
#ifndef DUNEURO_UNFITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
#define DUNEURO_UNFITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
#include <dune/common/parametertree.hh>
#include <duneuro/eeg/transfer_matrix_rhs_interface.hh>
#include <duneuro/eeg/unfitted_transfer_matrix_rhs.hh>
namespace duneuro
{
struct UnfittedTransferMatrixRHSFactory {
template <class Vector, class Solver>
static std::unique_ptr<TransferMatrixRHSInterface<typename Solver::Traits::FundamentalGridView,
Vector>>
create(const Solver& solver, const Dune::ParameterTree& config)
{
auto type = config.get<std::string>("type", "point");
if (type == "point") {
return std::make_unique<UnfittedTransferMatrixRHS<
typename Solver::Traits::FunctionSpace::GFS, typename Solver::Traits::SubTriangulation,
Vector>>(solver.functionSpace().getGFS(), solver.subTriangulation(),
config.get<std::size_t>("compartment"), solver.scaleToBBox());
} else {
DUNE_THROW(Dune::Exception, "unknown transfer matrix type \"" << type << "\"");
}
}
};
}
#endif // DUNEURO_UNFITTED_TRANSFER_MATRIX_RHS_FACTORY_HH
......@@ -24,7 +24,7 @@ namespace duneuro
using ProjectedPosition = duneuro::ProjectedPosition<Element, Coordinate>;
};
template <class S>
template <class S, class RHSFactory>
class UnfittedTransferMatrixSolver
{
public:
......@@ -32,11 +32,9 @@ namespace duneuro
UnfittedTransferMatrixSolver(
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<typename Traits::Solver> solver, bool scaleToBBox,
const Dune::ParameterTree& config)
std::shared_ptr<typename Traits::Solver> solver, const Dune::ParameterTree& config)
: subTriangulation_(subTriangulation)
, solver_(solver)
, scaleToBBox_(scaleToBBox)
, rightHandSideVector_(solver_->functionSpace().getGFS(), 0.0)
, config_(config)
{
......@@ -101,7 +99,6 @@ namespace duneuro
private:
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_;
std::shared_ptr<typename Traits::Solver> solver_;
bool scaleToBBox_;
#if HAVE_TBB
tbb::enumerable_thread_specific<typename Traits::RangeDOFVector> rightHandSideVector_;
#else
......@@ -119,21 +116,19 @@ namespace duneuro
Dune::Timer timer;
rightHandSideVector = 0.0;
// assemble right hand side
UnfittedTransferMatrixRHS<typename Traits::FunctionSpace::GFS,
typename Traits::SubTriangulation>
rhsAssembler(solver_->functionSpace().getGFS(), subTriangulation_,
config.get<std::size_t>("compartment"), scaleToBBox_);
auto rhsAssembler =
RHSFactory::template create<typename Traits::RangeDOFVector>(*solver_, config);
#if HAVE_TBB
{ // mutex, as the finite element map is not thread safe
tbb::mutex::scoped_lock lock(solver_->functionSpaceMutex());
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler->bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
}
rhsAssembler.assembleRightHandSide(rightHandSideVector);
rhsAssembler->assembleRightHandSide(rightHandSideVector);
#else
rhsAssembler.bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler.assembleRightHandSide(rightHandSideVector);
rhsAssembler->bind(reference.element, reference.localPosition, electrode.element,
electrode.localPosition);
rhsAssembler->assembleRightHandSide(rightHandSideVector);
#endif
timer.stop();
dataTree.set("time_rhs_assembly", timer.lastElapsed());
......
......@@ -27,6 +27,7 @@
#include <duneuro/eeg/dg_source_model_factory.hh>
#include <duneuro/eeg/eeg_forward_solver.hh>
#include <duneuro/eeg/electrode_projection_factory.hh>
#include <duneuro/eeg/fitted_transfer_matrix_rhs_factory.hh>
#include <duneuro/eeg/fitted_transfer_matrix_solver.hh>
#include <duneuro/eeg/fitted_transfer_matrix_user.hh>
#include <duneuro/io/fitted_tensor_vtk_functor.hh>
......@@ -66,6 +67,7 @@ namespace duneuro
typename SelectFittedSolver<solverType, VC, elementType, degree>::SolverBackendType;
using SourceModelFactory =
typename SelectFittedSolver<solverType, VC, elementType, degree>::SourceModelFactoryType;
using TransferMatrixRHSFactory = FittedTransferMatrixRHSFactory;
using DomainDOFVector = typename Solver::Traits::DomainDOFVector;
using ElementSearch = KDTreeElementSearch<typename VC::GridView>;
};
......@@ -430,7 +432,8 @@ namespace duneuro
#else
typename Traits::SolverBackend solverBackend_;
#endif
FittedTransferMatrixSolver<typename Traits::Solver> eegTransferMatrixSolver_;
FittedTransferMatrixSolver<typename Traits::Solver, typename Traits::TransferMatrixRHSFactory>
eegTransferMatrixSolver_;
FittedMEGTransferMatrixSolver<typename Traits::Solver> megTransferMatrixSolver_;
EEGForwardSolver<typename Traits::Solver, typename Traits::SourceModelFactory>
eegForwardSolver_;
......
......@@ -21,6 +21,7 @@
#include <duneuro/eeg/eeg_forward_solver.hh>
#include <duneuro/eeg/projected_electrodes.hh>
#include <duneuro/eeg/udg_source_model_factory.hh>
#include <duneuro/eeg/unfitted_transfer_matrix_rhs_factory.hh>
#include <duneuro/eeg/unfitted_transfer_matrix_solver.hh>
#include <duneuro/eeg/unfitted_transfer_matrix_user.hh>
#include <duneuro/io/refined_vtk_writer.hh>
......@@ -74,7 +75,8 @@ namespace duneuro
using Solver = typename SelectUnfittedSolver<solverType, dim, degree, compartments>::SolverType;
using SourceModelFactory = typename SelectUnfittedSolver<solverType, dim, degree,
compartments>::SourceModelFactoryType;
using EEGTransferMatrixSolver = UnfittedTransferMatrixSolver<Solver>;
using TransferMatrixRHSFactory = UnfittedTransferMatrixRHSFactory;
using EEGTransferMatrixSolver = UnfittedTransferMatrixSolver<Solver, TransferMatrixRHSFactory>;
using TransferMatrixUser = UnfittedTransferMatrixUser<Solver, SourceModelFactory>;
using SolverBackend =
typename SelectUnfittedSolver<solverType, dim, degree, compartments>::SolverBackendType;
......@@ -112,8 +114,7 @@ namespace duneuro
config.sub("solver")))
, solverBackend_(solver_,
config.hasSub("solver") ? config.sub("solver") : Dune::ParameterTree())
, eegTransferMatrixSolver_(subTriangulation_, solver_, Traits::scaleToBBox(),
config.sub("solver"))
, eegTransferMatrixSolver_(subTriangulation_, solver_, config.sub("solver"))
, eegForwardSolver_(solver_)
, conductivities_(config.get<std::vector<double>>("solver.conductivities"))
{
......
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