Commit 5b2e1178 authored by Andreas Nüßing's avatar Andreas Nüßing

Unify EEGForwardSolver

This unifies the fitted and unfitted EEGForwardSolver. Up to now, there
were two basically identical implementations, the
{Fitted,Unfitted}EEGForwardSolver. As they did basically the same thing, we
merge both implementations and parameterize them via the Solver and the
SourceModelFactory. While doing so, we move the element-search into the
Solver class and mark several variables as const that are nerver and also
should never be modified (mainly the VolumeConductor and the
ElementSearch).
parent 41a4d6ec
...@@ -69,8 +69,14 @@ namespace duneuro ...@@ -69,8 +69,14 @@ namespace duneuro
using Traits = AuxilliaryTraits<GV, bt, sc>; using Traits = AuxilliaryTraits<GV, bt, sc>;
using GFS = typename Traits::FS::GFS; using GFS = typename Traits::FS::GFS;
explicit CGFirstOrderSpace(typename Traits::Grid& grid, const typename Traits::GridView& gv) // note: the function space interface assumes a non-const reference, even though it does not
: gridview(gv), bctype(gridview, problem), fs(grid, bctype) // modify the grid as it only obtains a gridview. to avoid having to pass a non-const reference
// to this constructor, we use the const_cast below
explicit CGFirstOrderSpace(const typename Traits::Grid& grid,
const typename Traits::GridView& gv)
: gridview(gv)
, bctype(gridview, problem)
, fs(const_cast<typename Traits::Grid&>(grid), bctype)
{ {
fs.assembleConstraints(bctype); fs.assembleConstraints(bctype);
} }
......
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
#include <duneuro/common/assembler.hh> #include <duneuro/common/assembler.hh>
#include <duneuro/common/convection_diffusion_cg_default_parameter.hh> #include <duneuro/common/convection_diffusion_cg_default_parameter.hh>
#include <duneuro/common/flags.hh> #include <duneuro/common/flags.hh>
#include <duneuro/common/kdtree.hh>
#include <duneuro/common/linear_problem_solver.hh> #include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/make_dof_vector.hh> #include <duneuro/common/make_dof_vector.hh>
#include <duneuro/common/random.hh> #include <duneuro/common/random.hh>
...@@ -38,6 +39,8 @@ namespace duneuro ...@@ -38,6 +39,8 @@ namespace duneuro
struct CGSolverTraits { struct CGSolverTraits {
static const int dimension = VC::dim; static const int dimension = VC::dim;
using VolumeConductor = VC; using VolumeConductor = VC;
using CoordinateFieldType = typename VC::ctype;
using ElementSearch = KDTreeElementSearch<typename VC::GridView>;
using Problem = ConvectionDiffusionCGDefaultParameter<VC>; using Problem = ConvectionDiffusionCGDefaultParameter<VC>;
using DirichletExtension = Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem>; using DirichletExtension = Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem>;
using BoundaryCondition = Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem>; using BoundaryCondition = Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem>;
...@@ -59,12 +62,17 @@ namespace duneuro ...@@ -59,12 +62,17 @@ namespace duneuro
public: public:
using Traits = CGSolverTraits<VC, elementType, degree, DF, RF, JF>; using Traits = CGSolverTraits<VC, elementType, degree, DF, RF, JF>;
CGSolver(std::shared_ptr<VC> volumeConductor, const Dune::ParameterTree& config, CGSolver(std::shared_ptr<const VC> volumeConductor,
DataTree dataTree = DataTree()) std::shared_ptr<const typename Traits::ElementSearch> search,
: problem_(volumeConductor) const Dune::ParameterTree& config, DataTree dataTree = DataTree())
: volumeConductor_(volumeConductor)
, search_(search)
, problem_(volumeConductor)
, dirichletExtension_(volumeConductor->gridView(), problem_) , dirichletExtension_(volumeConductor->gridView(), problem_)
, boundaryCondition_(volumeConductor->gridView(), problem_) , boundaryCondition_(volumeConductor->gridView(), problem_)
, functionSpace_(volumeConductor->grid(), boundaryCondition_) , functionSpace_(const_cast<typename VC::GridType&>(volumeConductor->grid()),
boundaryCondition_)
// const_cast due to misused non-const reference in pdelab boilerplate
, localOperator_(problem_, config.get<unsigned int>("intorderadd", 0)) , localOperator_(problem_, config.get<unsigned int>("intorderadd", 0))
, assembler_(functionSpace_, localOperator_, elementType == ElementType::hexahedron ? , assembler_(functionSpace_, localOperator_, elementType == ElementType::hexahedron ?
(1 << VC::dim) + 1 : (1 << VC::dim) + 1 :
...@@ -91,7 +99,19 @@ namespace duneuro ...@@ -91,7 +99,19 @@ namespace duneuro
return functionSpace_; return functionSpace_;
} }
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor() const
{
return volumeConductor_;
}
std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
{
return search_;
}
private: private:
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
typename Traits::Problem problem_; typename Traits::Problem problem_;
typename Traits::DirichletExtension dirichletExtension_; typename Traits::DirichletExtension dirichletExtension_;
typename Traits::BoundaryCondition boundaryCondition_; typename Traits::BoundaryCondition boundaryCondition_;
......
...@@ -12,7 +12,7 @@ namespace duneuro ...@@ -12,7 +12,7 @@ namespace duneuro
using GV = typename VC::GridType::LeafGridView; using GV = typename VC::GridType::LeafGridView;
using Traits = Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype>; using Traits = Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype>;
explicit ConvectionDiffusionCGDefaultParameter(std::shared_ptr<VC> volumeConductor) explicit ConvectionDiffusionCGDefaultParameter(std::shared_ptr<const VC> volumeConductor)
: volumeConductor_(volumeConductor) : volumeConductor_(volumeConductor)
{ {
assert(volumeConductor_); assert(volumeConductor_);
...@@ -74,7 +74,7 @@ namespace duneuro ...@@ -74,7 +74,7 @@ namespace duneuro
} }
private: private:
std::shared_ptr<VC> volumeConductor_; std::shared_ptr<const VC> volumeConductor_;
}; };
} }
......
...@@ -15,7 +15,7 @@ namespace duneuro ...@@ -15,7 +15,7 @@ namespace duneuro
typedef typename VC::GridView GV; typedef typename VC::GridView GV;
typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype> Traits; typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype> Traits;
explicit ConvectionDiffusion_DG_DefaultParameter(std::shared_ptr<VC> volumeConductor) explicit ConvectionDiffusion_DG_DefaultParameter(std::shared_ptr<const VC> volumeConductor)
: volumeConductor_(volumeConductor) : volumeConductor_(volumeConductor)
{ {
assert(volumeConductor_); assert(volumeConductor_);
...@@ -94,7 +94,7 @@ namespace duneuro ...@@ -94,7 +94,7 @@ namespace duneuro
} }
private: private:
std::shared_ptr<VC> volumeConductor_; std::shared_ptr<const VC> volumeConductor_;
}; };
} }
......
...@@ -41,7 +41,7 @@ namespace duneuro ...@@ -41,7 +41,7 @@ namespace duneuro
GFS; GFS;
typedef typename Dune::PDELab::Backend::Vector<GFS, N> DOF; typedef typename Dune::PDELab::Backend::Vector<GFS, N> DOF;
CutFEMMultiPhaseSpace(const GV& gv, std::shared_ptr<SubTriangulation> subTriangulation) CutFEMMultiPhaseSpace(const GV& gv, std::shared_ptr<const SubTriangulation> subTriangulation)
: gridView_(gv) : gridView_(gv)
, entitySet_(gridView_) , entitySet_(gridView_)
, subTriangulation_(subTriangulation) , subTriangulation_(subTriangulation)
...@@ -73,7 +73,7 @@ namespace duneuro ...@@ -73,7 +73,7 @@ namespace duneuro
private: private:
GV gridView_; GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_; Dune::PDELab::AllEntitySet<GV> entitySet_;
std::shared_ptr<SubTriangulation> subTriangulation_; std::shared_ptr<const SubTriangulation> subTriangulation_;
LFE lfe_; LFE lfe_;
std::array<std::shared_ptr<FEM>, phases> fems_; std::array<std::shared_ptr<FEM>, phases> fems_;
std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_; std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_;
......
...@@ -14,6 +14,7 @@ ...@@ -14,6 +14,7 @@
#include <duneuro/common/cutfem_gridoperator.hh> #include <duneuro/common/cutfem_gridoperator.hh>
#include <duneuro/common/cutfem_multi_phase_space.hh> #include <duneuro/common/cutfem_multi_phase_space.hh>
#include <duneuro/common/edge_norm_provider.hh> #include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/kdtree.hh>
#include <duneuro/common/linear_problem_solver.hh> #include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/penalty_flux_weighting.hh> #include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/random.hh> #include <duneuro/common/random.hh>
...@@ -25,6 +26,8 @@ namespace duneuro ...@@ -25,6 +26,8 @@ namespace duneuro
struct CutFEMSolverTraits { struct CutFEMSolverTraits {
using SubTriangulation = ST; using SubTriangulation = ST;
using FundamentalGridView = typename ST::BaseT::GridView; using FundamentalGridView = typename ST::BaseT::GridView;
using CoordinateFieldType = typename FundamentalGridView::ctype;
using ElementSearch = KDTreeElementSearch<FundamentalGridView>;
static const int dimension = FundamentalGridView::dimension; static const int dimension = FundamentalGridView::dimension;
static const int compartments = comps; static const int compartments = comps;
using Problem = P; using Problem = P;
...@@ -58,18 +61,22 @@ namespace duneuro ...@@ -58,18 +61,22 @@ namespace duneuro
public: public:
using Traits = CutFEMSolverTraits<ST, compartments, degree, P, DF, RF, JF>; using Traits = CutFEMSolverTraits<ST, compartments, degree, P, DF, RF, JF>;
CutFEMSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation, CutFEMSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<const typename Traits::ElementSearch> search,
const Dune::ParameterTree& config) const Dune::ParameterTree& config)
: CutFEMSolver(subTriangulation, std::make_shared<typename Traits::Problem>( : CutFEMSolver(subTriangulation, search,
config.get<std::vector<double>>("conductivities")), std::make_shared<typename Traits::Problem>(
config.get<std::vector<double>>("conductivities")),
config) config)
{ {
} }
CutFEMSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation, CutFEMSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<const typename Traits::ElementSearch> search,
std::shared_ptr<typename Traits::Problem> problem, std::shared_ptr<typename Traits::Problem> problem,
const Dune::ParameterTree& config) const Dune::ParameterTree& config)
: subTriangulation_(subTriangulation) : subTriangulation_(subTriangulation)
, search_(search)
, problem_(problem) , problem_(problem)
, functionSpace_(subTriangulation_->gridView(), subTriangulation_) , functionSpace_(subTriangulation_->gridView(), subTriangulation_)
, edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0) , edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
...@@ -116,9 +123,9 @@ namespace duneuro ...@@ -116,9 +123,9 @@ namespace duneuro
return functionSpace_; return functionSpace_;
} }
const typename Traits::SubTriangulation& subTriangulation() const std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation() const
{ {
return *subTriangulation_; return subTriangulation_;
} }
typename Traits::Problem& problem() typename Traits::Problem& problem()
...@@ -126,6 +133,11 @@ namespace duneuro ...@@ -126,6 +133,11 @@ namespace duneuro
return *problem_; return *problem_;
} }
std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
{
return search_;
}
#if HAVE_TBB #if HAVE_TBB
tbb::mutex& functionSpaceMutex() tbb::mutex& functionSpaceMutex()
{ {
...@@ -134,7 +146,8 @@ namespace duneuro ...@@ -134,7 +146,8 @@ namespace duneuro
#endif #endif
private: private:
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_; std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
std::shared_ptr<typename Traits::Problem> problem_; std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_; typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_; typename Traits::EdgeNormProvider edgeNormProvider_;
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include <duneuro/common/convection_diffusion_dg_operator.hh> #include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/edge_norm_provider.hh> #include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/flags.hh> #include <duneuro/common/flags.hh>
#include <duneuro/common/kdtree.hh>
#include <duneuro/common/linear_problem_solver.hh> #include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/penalty_flux_weighting.hh> #include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/random.hh> #include <duneuro/common/random.hh>
...@@ -40,6 +41,8 @@ namespace duneuro ...@@ -40,6 +41,8 @@ namespace duneuro
struct DGSolverTraits { struct DGSolverTraits {
static const int dimension = VC::dim; static const int dimension = VC::dim;
using VolumeConductor = VC; using VolumeConductor = VC;
using CoordinateFieldType = typename VC::ctype;
using ElementSearch = KDTreeElementSearch<typename VC::GridView>;
using Problem = P; using Problem = P;
using FunctionSpace = typename DGFunctionSpaceTraits<VC, degree, elementType>::Type; using FunctionSpace = typename DGFunctionSpaceTraits<VC, degree, elementType>::Type;
using DomainDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, DF>; using DomainDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, DF>;
...@@ -61,9 +64,12 @@ namespace duneuro ...@@ -61,9 +64,12 @@ namespace duneuro
public: public:
using Traits = DGSolverTraits<VC, elementType, degree, P, DF, RF, JF>; using Traits = DGSolverTraits<VC, elementType, degree, P, DF, RF, JF>;
DGSolver(std::shared_ptr<VC> volumeConductor, std::shared_ptr<typename Traits::Problem> problem, DGSolver(std::shared_ptr<const VC> volumeConductor,
const Dune::ParameterTree& config, DataTree dataTree = DataTree()) std::shared_ptr<const typename Traits::ElementSearch> search,
std::shared_ptr<typename Traits::Problem> problem, const Dune::ParameterTree& config,
DataTree dataTree = DataTree())
: volumeConductor_(volumeConductor) : volumeConductor_(volumeConductor)
, search_(search)
, problem_(problem) , problem_(problem)
, functionSpace_(volumeConductor_->gridView()) , functionSpace_(volumeConductor_->gridView())
, edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0) , edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
...@@ -80,10 +86,11 @@ namespace duneuro ...@@ -80,10 +86,11 @@ namespace duneuro
dataTree.set("element_type", to_string(elementType)); dataTree.set("element_type", to_string(elementType));
} }
DGSolver(std::shared_ptr<VC> volumeConductor, const Dune::ParameterTree& config, DGSolver(std::shared_ptr<const VC> volumeConductor,
DataTree dataTree = DataTree()) std::shared_ptr<const typename Traits::ElementSearch> search,
: DGSolver(volumeConductor, std::make_shared<typename Traits::Problem>(volumeConductor), const Dune::ParameterTree& config, DataTree dataTree = DataTree())
config, dataTree) : DGSolver(volumeConductor, search,
std::make_shared<typename Traits::Problem>(volumeConductor), config, dataTree)
{ {
} }
...@@ -113,14 +120,9 @@ namespace duneuro ...@@ -113,14 +120,9 @@ namespace duneuro
return functionSpace_; return functionSpace_;
} }
const typename Traits::VolumeConductor& volumeConductor() const std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor() const
{
return *volumeConductor_;
}
typename Traits::VolumeConductor& volumeConductor()
{ {
return *volumeConductor_; return volumeConductor_;
} }
const typename Traits::Assembler& assembler() const const typename Traits::Assembler& assembler() const
...@@ -143,8 +145,14 @@ namespace duneuro ...@@ -143,8 +145,14 @@ namespace duneuro
return *problem_; return *problem_;
} }
std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
{
return search_;
}
private: private:
std::shared_ptr<VC> volumeConductor_; std::shared_ptr<const VC> volumeConductor_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
std::shared_ptr<typename Traits::Problem> problem_; std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_; typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_; typename Traits::EdgeNormProvider edgeNormProvider_;
......
...@@ -25,7 +25,7 @@ namespace duneuro ...@@ -25,7 +25,7 @@ namespace duneuro
explicit DGSolverBackend(std::shared_ptr<Solver> solver, const Dune::ParameterTree& config) explicit DGSolverBackend(std::shared_ptr<Solver> solver, const Dune::ParameterTree& config)
: solver_(solver) : solver_(solver)
, firstOrderSpace_(solver->volumeConductor().grid(), solver->volumeConductor().gridView()) , firstOrderSpace_(solver->volumeConductor()->grid(), solver->volumeConductor()->gridView())
, solverBackend_(*solver->assembler(), firstOrderSpace_.getGFS(), config) , solverBackend_(*solver->assembler(), firstOrderSpace_.getGFS(), config)
, config_(config) , config_(config)
{ {
......
...@@ -201,7 +201,7 @@ namespace duneuro ...@@ -201,7 +201,7 @@ namespace duneuro
template <class VC, class ES> template <class VC, class ES>
std::function<bool(typename VC::EntityType)> std::function<bool(typename VC::EntityType)>
make_element_filter(std::shared_ptr<VC> volumeConductor, const ES& elementSearch, make_element_filter(std::shared_ptr<const VC> volumeConductor, const ES& elementSearch,
const Dune::FieldVector<typename VC::ctype, VC::dim>& position, bool restrict) const Dune::FieldVector<typename VC::ctype, VC::dim>& position, bool restrict)
{ {
if (restrict) { if (restrict) {
...@@ -218,7 +218,7 @@ namespace duneuro ...@@ -218,7 +218,7 @@ namespace duneuro
template <class VC, class ES> template <class VC, class ES>
std::unique_ptr<ElementPatch<typename VC::GridView>> make_element_patch( std::unique_ptr<ElementPatch<typename VC::GridView>> make_element_patch(
std::shared_ptr<VC> volumeConductor, std::shared_ptr<const VC> volumeConductor,
std::shared_ptr<ElementNeighborhoodMap<typename VC::GridView>> elementNeighborhoodMap, std::shared_ptr<ElementNeighborhoodMap<typename VC::GridView>> elementNeighborhoodMap,
const ES& elementSearch, const Dune::FieldVector<typename VC::ctype, VC::dim>& position, const ES& elementSearch, const Dune::FieldVector<typename VC::ctype, VC::dim>& position,
const Dune::ParameterTree& config) const Dune::ParameterTree& config)
......
...@@ -42,7 +42,7 @@ namespace duneuro ...@@ -42,7 +42,7 @@ namespace duneuro
typedef Dune::PDELab::PowerGridFunctionSpace<DomainGFS, phases, PVBE, OrderingTag> GFS; typedef Dune::PDELab::PowerGridFunctionSpace<DomainGFS, phases, PVBE, OrderingTag> GFS;
typedef typename Dune::PDELab::Backend::Vector<GFS, N> DOF; typedef typename Dune::PDELab::Backend::Vector<GFS, N> DOF;
UDGQkMultiPhaseSpace(const GV& gv, std::shared_ptr<SubTriangulation> subTriangulation) UDGQkMultiPhaseSpace(const GV& gv, std::shared_ptr<const SubTriangulation> subTriangulation)
: gridView_(gv), entitySet_(gridView_), subTriangulation_(subTriangulation) : gridView_(gv), entitySet_(gridView_), subTriangulation_(subTriangulation)
{ {
for (unsigned int i = 0; i < phases; ++i) { for (unsigned int i = 0; i < phases; ++i) {
...@@ -71,7 +71,7 @@ namespace duneuro ...@@ -71,7 +71,7 @@ namespace duneuro
private: private:
GV gridView_; GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_; Dune::PDELab::AllEntitySet<GV> entitySet_;
std::shared_ptr<SubTriangulation> subTriangulation_; std::shared_ptr<const SubTriangulation> subTriangulation_;
LFE lfe_; LFE lfe_;
std::array<std::shared_ptr<FEM>, phases> fems_; std::array<std::shared_ptr<FEM>, phases> fems_;
std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_; std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_;
......
...@@ -15,6 +15,7 @@ ...@@ -15,6 +15,7 @@
#include <duneuro/common/convection_diffusion_dg_operator.hh> #include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/convection_diffusion_udg_default_parameter.hh> #include <duneuro/common/convection_diffusion_udg_default_parameter.hh>
#include <duneuro/common/edge_norm_provider.hh> #include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/kdtree.hh>
#include <duneuro/common/linear_problem_solver.hh> #include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/penalty_flux_weighting.hh> #include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/udg_multi_phase_space.hh> #include <duneuro/common/udg_multi_phase_space.hh>
...@@ -26,6 +27,8 @@ namespace duneuro ...@@ -26,6 +27,8 @@ namespace duneuro
struct UDGSolverTraits { struct UDGSolverTraits {
using SubTriangulation = ST; using SubTriangulation = ST;
using FundamentalGridView = typename ST::BaseT::GridView; using FundamentalGridView = typename ST::BaseT::GridView;
using CoordinateFieldType = typename FundamentalGridView::ctype;
using ElementSearch = KDTreeElementSearch<FundamentalGridView>;
static const int dimension = FundamentalGridView::dimension; static const int dimension = FundamentalGridView::dimension;
static const int compartments = comps; static const int compartments = comps;
using Problem = P; using Problem = P;
...@@ -56,17 +59,21 @@ namespace duneuro ...@@ -56,17 +59,21 @@ namespace duneuro
public: public:
using Traits = UDGSolverTraits<ST, compartments, degree, P, DF, RF, JF>; using Traits = UDGSolverTraits<ST, compartments, degree, P, DF, RF, JF>;
UDGSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation, UDGSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<const typename Traits::ElementSearch> search,
const Dune::ParameterTree& config) const Dune::ParameterTree& config)
: UDGSolver(subTriangulation, std::make_shared<typename Traits::Problem>( : UDGSolver(subTriangulation, search,
config.get<std::vector<double>>("conductivities")), std::make_shared<typename Traits::Problem>(
config.get<std::vector<double>>("conductivities")),
config) config)
{ {
} }
UDGSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation, UDGSolver(std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<const typename Traits::ElementSearch> search,
std::shared_ptr<typename Traits::Problem> problem, const Dune::ParameterTree& config) std::shared_ptr<typename Traits::Problem> problem, const Dune::ParameterTree& config)
: subTriangulation_(subTriangulation) : subTriangulation_(subTriangulation)
, search_(search)
, problem_(problem) , problem_(problem)
, functionSpace_(subTriangulation_->gridView(), subTriangulation_) , functionSpace_(subTriangulation_->gridView(), subTriangulation_)
, edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0) , edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
...@@ -112,9 +119,9 @@ namespace duneuro ...@@ -112,9 +119,9 @@ namespace duneuro
return functionSpace_; return functionSpace_;
} }
const typename Traits::SubTriangulation& subTriangulation() const std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation() const
{ {
return *subTriangulation_; return subTriangulation_;
} }
typename Traits::Problem& problem() typename Traits::Problem& problem()
...@@ -122,6 +129,11 @@ namespace duneuro ...@@ -122,6 +129,11 @@ namespace duneuro
return *problem_; return *problem_;
} }
std::shared_ptr<const typename Traits::ElementSearch> elementSearch() const
{
return search_;
}
#if HAVE_TBB #if HAVE_TBB
tbb::mutex& functionSpaceMutex() tbb::mutex& functionSpaceMutex()
{ {
...@@ -130,7 +142,8 @@ namespace duneuro ...@@ -130,7 +142,8 @@ namespace duneuro
#endif #endif
private: private:
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_; std::shared_ptr<const typename Traits::SubTriangulation> subTriangulation_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
std::shared_ptr<typename Traits::Problem> problem_; std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_; typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_; typename Traits::EdgeNormProvider edgeNormProvider_;
......
...@@ -16,88 +16,93 @@ ...@@ -16,88 +16,93 @@
namespace duneuro namespace duneuro
{ {
struct CGSourceModelFactory { struct CGSourceModelFactory {
template <class Vector, class VC, class Solver> template <class Vector, class Solver>
static std::shared_ptr<SourceModelInterface<typename VC::ctype, VC::dim, Vector>> static std::shared_ptr<SourceModelInterface<typename Solver::Traits::VolumeConductor::ctype,
createDense(std::shared_ptr<VC> volumeConductor, const Solver& solver, Solver::Traits::VolumeConductor::dim, Vector>>
std::shared_ptr<KDTreeElementSearch<typename VC::GridView>> search, createDense(const Solver& solver, const Dune::ParameterTree& config,
const Dune::ParameterTree& config, const Dune::ParameterTree& solverConfig) const Dune::ParameterTree& solverConfig)
{ {
const auto type = config.get<std::string>("type"); const auto type = config.get<std::string>("type");
if (type == "partial_integration") { if (type == "partial_integration") {
return std::make_shared<PartialIntegrationSourceModel< return std::make_shared<PartialIntegrationSourceModel<
typename Solver::Traits::FunctionSpace::GFS, Vector>>(solver.functionSpace().getGFS(), typename Solver::Traits::FunctionSpace::GFS, Vector>>(solver.functionSpace().getGFS(),
search); solver.elementSearch());
} else if (type == "venant") { } else if (type == "venant") {
return std::make_shared<VertexBasedVenantSourceModel<VC, typename Solver::Traits:: return std::make_shared<VertexBasedVenantSourceModel<
FunctionSpace::GFS, typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace::GFS,
Vector>>( Vector>>(solver.volumeConductor(), solver.functionSpace().getGFS(),
volumeConductor, solver.functionSpace().getGFS(), search, config); solver.elementSearch(), config);
} else if (type == "patch_based_venant") { } else if (type == "patch_based_venant") {
return std:: return std::make_shared<PatchBasedVenantSourceModel<
make_shared<PatchBasedVenantSourceModel<VC, typename Solver::Traits::FunctionSpace::GFS, typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace::GFS,
Vector>>( Vector>>(solver.volumeConductor(), solver.functionSpace().getGFS(),
volumeConductor, solver.functionSpace().getGFS(), search, config); solver.elementSearch(), config);
} else if (type == "spatial_venant") { } else if (type == "spatial_venant") {
return std:: return std::make_shared<SpatialVenantSourceModel<
make_shared<SpatialVenantSourceModel<VC, typename Solver::Traits::FunctionSpace::GFS, typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace::GFS,
Vector>>( Vector>>(solver.volumeConductor(), solver.functionSpace().getGFS(),
volumeConductor, solver.functionSpace().getGFS(), search, config); solver.elementSearch(), config);
} else if (type == "truncated_spatial_venant") { } else if (type == "truncated_spatial_venant") {
return std::make_shared<TruncatedSpatialVenantSourceModel<VC, typename Solver::Traits:: return std::make_shared<TruncatedSpatialVenantSourceModel<
FunctionSpace::GFS, typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace::GFS,
Vector>>( Vector>>(solver.volumeConductor(), solver.functionSpace().getGFS(),
volumeConductor, solver.functionSpace().getGFS(), search, config); solver.elementSearch(), config);
} else if (type == "subtraction") { } else if (type == "subtraction") {
return std::make_shared<FittedSubtractionSourceModel< return std::make_shared<FittedSubtractionSourceModel<
typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace, typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace,
Vector, SubtractionContinuityType::continuous>>(volumeConductor, solver.functionSpace(), Vector, SubtractionContinuityType::continuous>>(
search, config, solverConfig); solver.volumeConductor(), solver.functionSpace(), solver.elementSearch(), config,
solverConfig);
} else if (type == "whitney") { } else if (type == "whitney") {
return std::make_shared<WhitneySourceModel<VC, typename Solver::Traits::FunctionSpace::GFS, return std::make_shared<WhitneySourceModel<typename Solver::Traits::VolumeConductor,
Vector>>( typename Solver::Traits::FunctionSpace::GFS,
volumeConductor, solver.functionSpace().getGFS(), search, config); Vector>>(solver.volumeConductor(),
solver.functionSpace().getGFS(),
solver.elementSearch(), config);
} else { } else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\""); << "\"");
} }
} }
template <class Vector, class VC, class Solver> template <class Vector, class Solver>
static std::shared_ptr<SourceModelInterface<typename VC::ctype, VC::dim,