Commit 41a4d6ec authored by Andreas Nüßing's avatar Andreas Nüßing

Merge branch 'feature/udg-subtraction' into 'master'

[UDG] Subtraction

See merge request duneuro/duneuro!51
parents 91d41192 4deb2a0a
......@@ -12,7 +12,7 @@ namespace duneuro
{
public:
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
typedef typename VC::GridType::LeafGridView GV;
typedef typename VC::GridView GV;
typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV, typename GV::ctype> Traits;
explicit ConvectionDiffusion_DG_DefaultParameter(std::shared_ptr<VC> volumeConductor)
......
......@@ -16,6 +16,7 @@ namespace duneuro
typedef typename ES::template Codim<0>::Entity EntityType;
typedef Dune::FieldMatrix<ctype, dim, dim> TensorType;
typedef ES EntitySet;
typedef ES GridView;
EntitySetVolumeConductor(const ES& entitySet, const std::vector<TensorType>& tensors)
: entitySet_(entitySet), tensors_(tensors), elementMapper_(entitySet_)
......
......@@ -4,16 +4,15 @@
#include <dune/common/parametertree.hh>
#include <duneuro/common/exceptions.hh>
#include <duneuro/eeg/fitted_subtraction_source_model.hh>
#include <duneuro/eeg/partial_integration_source_model.hh>
#include <duneuro/eeg/patch_based_venant_source_model.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/spatial_venant_source_model.hh>
#include <duneuro/eeg/subtraction_source_model.hh>
#include <duneuro/eeg/truncated_spatial_venant_source_model.hh>
#include <duneuro/eeg/vertex_based_venant_source_model.hh>
#include <duneuro/eeg/whitney_source_model.hh>
namespace duneuro
{
struct CGSourceModelFactory {
......@@ -50,7 +49,7 @@ namespace duneuro
Vector>>(
volumeConductor, solver.functionSpace().getGFS(), search, config);
} else if (type == "subtraction") {
return std::make_shared<SubtractionSourceModel<
return std::make_shared<FittedSubtractionSourceModel<
typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace,
Vector, SubtractionContinuityType::continuous>>(volumeConductor, solver.functionSpace(),
search, config, solverConfig);
......
......@@ -17,7 +17,8 @@ namespace duneuro
createDense(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
......@@ -36,7 +37,8 @@ namespace duneuro
createSparse(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
......
......@@ -4,10 +4,10 @@
#include <dune/common/parametertree.hh>
#include <duneuro/common/exceptions.hh>
#include <duneuro/eeg/fitted_subtraction_source_model.hh>
#include <duneuro/eeg/localized_subtraction_source_model.hh>
#include <duneuro/eeg/partial_integration_source_model.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/subtraction_source_model.hh>
namespace duneuro
{
......@@ -29,10 +29,10 @@ namespace duneuro
V>>(
volumeConductor, solver.functionSpace().getGFS(), search, config);
} else if (type == "subtraction") {
return std::make_shared<SubtractionSourceModel<typename Solver::Traits::VolumeConductor,
typename Solver::Traits::FunctionSpace, V,
SubtractionContinuityType::discontinuous>>(
volumeConductor, solver.functionSpace(), search, config, solverConfig);
return std::make_shared<FittedSubtractionSourceModel<
typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace, V,
SubtractionContinuityType::discontinuous>>(volumeConductor, solver.functionSpace(),
search, config, solverConfig);
} else if (type == "localized_subtraction") {
return std::make_shared<LocalizedSubtractionSourceModel<
typename Solver::Traits::VolumeConductor, typename Solver::Traits::FunctionSpace, V>>(
......
#ifndef DUNEURO_SUBTRACTIONDGRESIDUAL_HH
#define DUNEURO_SUBTRACTIONDGRESIDUAL_HH
#ifndef DUNEURO_FITTED_SUBTRACTION_SOURCE_MODEL_HH
#define DUNEURO_FITTED_SUBTRACTION_SOURCE_MODEL_HH
#include <dune/common/parametertree.hh>
......@@ -15,7 +15,8 @@
namespace duneuro
{
template <class VC, class FS, class V, SubtractionContinuityType continuityType>
class SubtractionSourceModel : public SourceModelBase<typename FS::GFS::Traits::GridViewType, V>
class FittedSubtractionSourceModel
: public SourceModelBase<typename FS::GFS::Traits::GridViewType, V>
{
public:
using BaseT = SourceModelBase<typename FS::GFS::Traits::GridViewType, V>;
......@@ -32,9 +33,10 @@ namespace duneuro
using VectorType = typename BaseT::VectorType;
using SearchType = typename BaseT::SearchType;
SubtractionSourceModel(std::shared_ptr<VC> volumeConductor, const FS& fs,
std::shared_ptr<SearchType> search, const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
FittedSubtractionSourceModel(std::shared_ptr<VC> volumeConductor, const FS& fs,
std::shared_ptr<SearchType> search,
const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
: BaseT(search)
, problem_(volumeConductor->gridView(), volumeConductor)
, edgeNormProvider_(solverConfig.get<std::string>("edge_norm_type", "houston"), 1.0)
......@@ -94,4 +96,4 @@ namespace duneuro
};
}
#endif // DUNEURO_SUBTRACTIONDGRESIDUAL_HH
#endif // DUNEURO_FITTED_SUBTRACTION_SOURCE_MODEL_HH
......@@ -16,10 +16,10 @@
/**** includes ****/
#include <dune/common/timer.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
/**** our includes ****/
#include <duneuro/common/convection_diffusion_dg_default_parameter.hh>
#include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/dipole.hh>
#include <duneuro/eeg/subtraction_dg_uinfty.hh>
......@@ -28,12 +28,13 @@ namespace duneuro
{
/**** class definition ****/
template <typename GV, typename RF, typename VC>
class SubtractionDGDefaultParameter
class SubtractionDGDefaultParameter : public ConvectionDiffusion_DG_DefaultParameter<VC>
{
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
public:
/*** Typedefs ***/
using BaseT = ConvectionDiffusion_DG_DefaultParameter<VC>;
typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV, RF> Traits;
typedef typename Traits::GridViewType::Traits::IndexSet IndexSet;
typedef typename Traits::GridViewType::Traits::Grid GridType;
......@@ -41,114 +42,10 @@ namespace duneuro
/*** Constructor ***/
SubtractionDGDefaultParameter(const typename Traits::GridViewType& gv_,
std::shared_ptr<VC> volumeConductor)
: gv(gv_), volumeConductor_(volumeConductor), mapper(gv), u_infty(gv), grad_u_infty(gv)
: BaseT(volumeConductor), gv(gv_), u_infty(gv), grad_u_infty(gv)
{
}
/*** functions required by the ConvectionDiffusionParameterInterface ***/
/** source/reaction term **/
typename Traits::RangeType f(const typename Traits::ElementType& e,
const typename Traits::DomainType& x) const
{
return typename Traits::RangeType(0.0);
}
/** nonlinearity under the gradient **/
typename Traits::RangeType b(const typename Traits::ElementType& e,
const typename Traits::DomainType& x) const
{
typename Traits::RangeType ret;
ret[0] = 0.0;
ret[1] = 0.0;
ret[2] = 0.0;
return ret;
}
template <class IG>
typename Traits::RangeType b(const IG& ig, const typename Traits::IntersectionDomainType& x,
ConvectionDiffusion_DG_Side::Type side) const
{
switch (side) {
case ConvectionDiffusion_DG_Side::inside: return typename Traits::RangeType(0.0);
default: return typename Traits::RangeType(0.0);
}
}
typename Traits::RangeFieldType c(const typename Traits::ElementType& e,
const typename Traits::DomainType& x) const
{
return 0.0;
}
/** tensor diffusion coefficient **/
template <class EG>
typename Traits::PermTensorType A(const EG& e, const typename Traits::DomainType& x) const
{
typename Traits::PermTensorType ret = sigma(e.entity());
/* see the definition of the problem in convectiondiffusiondg.hh */
// ret *= -1;
return ret;
}
typename Traits::PermTensorType A(const typename Traits::ElementType& e,
const typename Traits::DomainType& x) const
{
typename Traits::PermTensorType ret = sigma(e);
/* see the definition of the problem in convectiondiffusiondg.hh */
// ret *= -1;
return ret;
}
template <class IG>
typename Traits::PermTensorType A(const IG& ig,
const typename Traits::IntersectionDomainType& x,
ConvectionDiffusion_DG_Side::Type side) const
{
switch (side) {
case ConvectionDiffusion_DG_Side::inside:
return A(ig.inside(), ig.geometryInInside().global(x));
default: return A(ig.outside(), ig.geometryInOutside().global(x));
}
}
/** dirichlet boundary or not **/
template <typename I>
bool isDirichlet(const I& intersection,
const Dune::FieldVector<typename I::ctype, I::dimension - 1>& coord) const
{
return false;
}
/** dirichlet boundary or not **/
template <typename I>
BCType bctype(const I& intersection,
const Dune::FieldVector<typename I::ctype, I::dimension - 1>& coord) const
{
return BCType::Neumann;
}
/** dirichlet boundary condition value **/
/** irrelevant for us(see above) but might be needed by the implementation?! **/
typename Traits::RangeFieldType g(const typename Traits::ElementType& e,
const typename Traits::DomainType& x) const
{
return 0.0;
}
template <typename I>
typename Traits::RangeFieldType
g(const I& intersection,
const Dune::FieldVector<typename I::ctype, I::dimension - 1>& coord) const
{
return 0.0;
}
template <typename I>
typename Traits::RangeFieldType
o(const I& intersection,
const Dune::FieldVector<typename I::ctype, I::dimension - 1>& coord) const
{
return 0.0;
}
/** Neumann boundary condition **/
typename Traits::RangeFieldType j(const typename Traits::IntersectionType& e,
const typename Traits::IntersectionDomainType& x)
{
......@@ -174,22 +71,6 @@ namespace duneuro
return temp * normal;
}
/*** functions ***/
/** returns the conductivity label for a given grid element **/
typename Traits::PermTensorType sigma(const typename Traits::ElementType& e) const
{
return volumeConductor_->tensor(e);
}
/** get sigma_corr as a PermTensorType **/
typename Traits::PermTensorType sigma_corr(const typename Traits::ElementType& e) const
{
typename Traits::PermTensorType ret(sigma(e));
ret -= sigma_infty;
return ret;
}
/** evaluate graduinfty and uinfty at global coordinates **/
typename Traits::RangeType get_grad_u_infty(const typename Traits::DomainType& x) const
{
typename Traits::RangeType ret;
......@@ -214,11 +95,6 @@ namespace duneuro
return sigma_infty_inv;
}
typename Traits::PermTensorType get_sigma(const typename Traits::ElementType& e)
{
return sigma(e);
}
typename Traits::GridViewType& get_gridview()
{
return gv;
......@@ -230,7 +106,7 @@ namespace duneuro
const typename Traits::DomainType& dipoleMoment)
{
/* find the new source's surrounding and the conductivity in it */
sigma_infty = sigma(element);
sigma_infty = this->A(element, localDipolePosition);
sigma_infty_inv = sigma_infty;
sigma_infty_inv.invert();
......@@ -252,15 +128,9 @@ namespace duneuro
typename Traits::GridViewType gv;
/*** the parameters for the forward simulation(dipole moment, position, etc.) ***/
std::shared_ptr<VC> volumeConductor_;
typename Traits::PermTensorType sigma_infty;
typename Traits::PermTensorType sigma_infty_inv;
/*** mapper to obtain the index of an element ***/
const Dune::MultipleCodimMultipleGeomTypeMapper<typename Traits::GridViewType,
Dune::MCMGElementLayout>
mapper;
/*** u_infty and its gradient as analytic grid functions ***/
InfinityPotential<typename Traits::GridViewType, typename Traits::RangeFieldType> u_infty;
InfinityPotentialGradient<typename Traits::GridViewType, typename Traits::RangeFieldType>
......
This diff is collapsed.
#ifndef DUNEURO_SUBTRACTION_UDG_DEFAULT_PARAMETER_HH
#define DUNEURO_SUBTRACTION_UDG_DEFAULT_PARAMETER_HH
#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
#include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/convection_diffusion_udg_default_parameter.hh>
#include <duneuro/common/dipole.hh>
#include <duneuro/eeg/subtraction_dg_uinfty.hh>
namespace duneuro
{
template <typename GV, typename RF>
class SubtractionUDGDefaultParameter : public ConvectionDiffusion_UDG_DefaultParameter<GV>
{
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
public:
typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV, RF> Traits;
using BaseT = ConvectionDiffusion_UDG_DefaultParameter<GV>;
explicit SubtractionUDGDefaultParameter(const GV& gv_,
const std::vector<double>& conductivities,
unsigned int dipolePhase)
: BaseT(conductivities), gv(gv_), dipolePhase_(dipolePhase), u_infty(gv), grad_u_infty(gv)
{
}
/** Neumann boundary condition **/
template <class IG>
typename Traits::RangeFieldType j(const IG& ig,
const typename Traits::IntersectionDomainType& local) const
{
/* map position on the intersection x(2D coordinates) to global position global_x in
* the grid(3D coordinates) for evaluation of graduinfty.
*/
typename Traits::DomainType global_x = ig.geometry().global(local);
/* evaluate graduinfty*/
Dune::FieldVector<typename Traits::RangeFieldType, GV::dimension> graduinfty;
grad_u_infty.evaluateGlobal(global_x, graduinfty);
/* temporary variable for the calculations */
Dune::FieldVector<typename Traits::RangeFieldType, GV::dimension> temp;
/* sigma^{\infty}\cdot\nabla u^{\infty} */
sigma_infty.mv(graduinfty, temp);
/* normal vector at the integration point */
Dune::FieldVector<typename Traits::RangeFieldType, GV::dimension> normal =
ig.unitOuterNormal(local);
return temp * normal;
}
/** evaluate graduinfty and uinfty at global coordinates **/
typename Traits::RangeType get_grad_u_infty(const typename Traits::DomainType& x) const
{
typename Traits::RangeType ret;
grad_u_infty.evaluateGlobal(x, ret);
return ret;
}
typename Traits::RangeFieldType get_u_infty(const typename Traits::DomainType& x) const
{
typename Dune::FieldVector<double, 1> ret;
u_infty.evaluateGlobal(x, ret);
return ret;
}
/** multiple helper functions that return private variables **/
typename Traits::PermTensorType get_sigma_infty()
{
return sigma_infty;
}
typename Traits::PermTensorType get_sigma_infty_inv()
{
return sigma_infty_inv;
}
/** set the the dipole position and moment **/
void bind(const typename Traits::ElementType& element,
const typename Traits::DomainType& localDipolePosition,
const typename Traits::DomainType& dipoleMoment)
{
/* find the new source's surrounding and the conductivity in it */
sigma_infty = this->A(dipolePhase_);
sigma_infty_inv = sigma_infty;
sigma_infty_inv.invert();
auto global = element.geometry().global(localDipolePosition);
/** set the values for the analytic grid function u_infty and its gradient **/
u_infty.set_parameters(dipoleMoment, global, sigma_infty, sigma_infty_inv);
grad_u_infty.set_parameters(dipoleMoment, global, sigma_infty, sigma_infty_inv);
}
const InfinityPotential<typename Traits::GridViewType, typename Traits::RangeFieldType>&
get_u_infty() const
{
return u_infty;
}
private:
GV gv;
unsigned int dipolePhase_;
/*** the parameters for the forward simulation(dipole moment, position, etc.) ***/
typename Traits::PermTensorType sigma_infty;
typename Traits::PermTensorType sigma_infty_inv;
/*** u_infty and its gradient as analytic grid functions ***/
InfinityPotential<typename Traits::GridViewType, typename Traits::RangeFieldType> u_infty;
InfinityPotentialGradient<typename Traits::GridViewType, typename Traits::RangeFieldType>
grad_u_infty;
};
}
#endif // DUNEURO_SUBTRACTION_UDG_DEFAULT_PARAMETER_HH
......@@ -6,6 +6,7 @@
#include <duneuro/common/exceptions.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/udg_subtraction_source_model.hh>
#include <duneuro/eeg/unfitted_partial_integration_source_model.hh>
#include <duneuro/eeg/unfitted_patch_based_venant_source_model.hh>
......@@ -18,7 +19,8 @@ namespace duneuro
createDense(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
......@@ -29,6 +31,11 @@ namespace duneuro
return Dune::Std::make_unique<UnfittedPatchBasedVenantSourceModel<
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, config);
} else if (type == "subtraction") {
return Dune::Std::make_unique<UDGSubtractionSourceModel<
typename Solver::Traits::FunctionSpace, ST, Vector>>(
solver.functionSpace(), subTriangulation, search, dipoleCompartment, config,
solverConfig);
} else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\"");
......@@ -41,7 +48,8 @@ namespace duneuro
createSparse(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
......
#ifndef DUNEURO_UDG_SUBTRACTION_SOURCE_MODEL_HH
#define DUNEURO_UDG_SUBTRACTION_SOURCE_MODEL_HH
#include <dune/common/parametertree.hh>
#include <dune/pdelab/backend/interface.hh>
#include <dune/pdelab/boilerplate/pdelab.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/subtraction_dg_operator.hh>
#include <duneuro/eeg/subtraction_udg_default_parameter.hh>
namespace duneuro
{
template <class FS, class ST, class V>
class UDGSubtractionSourceModel
: public SourceModelBase<typename FS::GFS::Traits::GridViewType, V>
{
public:
using GV = typename FS::GFS::Traits::GridViewType;
using DF = typename GV::ctype;
using RF = typename V::field_type;
using JF = RF;
using BaseT = SourceModelBase<GV, V>;
enum { dim = GV::dimension };
using Problem = SubtractionUDGDefaultParameter<GV, RF>;
using EdgeNormProvider = MultiEdgeNormProvider;
using PenaltyFluxWeighting = UnfittedDynamicPenaltyFluxWeights;
using LOP = SubtractionDG<Problem, EdgeNormProvider, PenaltyFluxWeighting,
SubtractionContinuityType::discontinuous>;
using WLOP = Dune::UDG::MultiPhaseLocalOperatorWrapper<LOP>;
using DOF = typename FS::DOF;
using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<GV>;
using MatrixBackend = Dune::PDELab::istl::BCRSMatrixBackend<>;
using GridOperator =
Dune::UDG::UDGGridOperator<typename FS::GFS, typename FS::GFS, WLOP, MatrixBackend, DF, RF,
JF, UnfittedSubTriangulation>;
using ElementType = typename BaseT::ElementType;
using CoordinateType = typename BaseT::CoordinateType;
using VectorType = typename BaseT::VectorType;
using SearchType = typename BaseT::SearchType;
UDGSubtractionSourceModel(const FS& fs, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<SearchType> search, unsigned int dipolePhase,
const Dune::ParameterTree& config,
const Dune::ParameterTree& solverConfig)
: BaseT(search)
, subTriangulation_(subTriangulation)
, problem_(subTriangulation->gridView(),
solverConfig.get<std::vector<double>>("conductivities"), dipolePhase)
, edgeNormProvider_(solverConfig.get<std::string>("edge_norm_type", "houston"), 1.0)
, weighting_(solverConfig.get<std::string>("weights", "tensorOnly"))
, lop_(problem_, weighting_, config.get<unsigned int>("intorderadd"),
config.get<unsigned int>("intorderadd_lb"))
, wlop_(lop_)
, ust_(subTriangulation->gridView(), *subTriangulation_)
, gridOperator_(fs.getGFS(), fs.getGFS(), ust_, wlop_, MatrixBackend(2 * dim + 1))
, x_(fs.getGFS(), 0.0)
, res_(fs.getGFS(), 0.0)
, interp_(fs.getGFS(), 0.0)
{
}
virtual void bind(const typename BaseT::DipoleType& dipole,
DataTree dataTree = DataTree()) override
{
BaseT::bind(dipole);
problem_.bind(this->dipoleElement(), this->localDipolePosition(), this->dipole().moment());
}
virtual void assembleRightHandSide(VectorType& vector) const override
{
x_ = 0.0;
gridOperator_.residual(x_, vector);
vector *= -1.0;
}
virtual void postProcessSolution(VectorType& vector) const override
{
DUNE_THROW(Dune::NotImplemented, "post processing of function not yet implemented");
}
virtual void
postProcessSolution(const std::vector<CoordinateType>& electrodes,
std::vector<typename VectorType::field_type>& vector) const override
{
assert(electrodes.size() == vector.size());
Dune::FieldVector<typename Problem::Traits::RangeFieldType, 1> result;
for (unsigned int i = 0; i < electrodes.size(); ++i) {
problem_.get_u_infty().evaluateGlobal(electrodes[i], result);
vector[i] += result;
}
}
private:
std::shared_ptr<ST> subTriangulation_;
Problem problem_;
EdgeNormProvider edgeNormProvider_;
PenaltyFluxWeighting weighting_;
LOP lop_;
WLOP wlop_;
UnfittedSubTriangulation ust_;
GridOperator gridOperator_;
mutable DOF x_;
mutable DOF res_;
mutable DOF interp_;
};
}
#endif // DUNEURO_UDG_SUBTRACTION_SOURCE_MODEL_HH
......@@ -43,10 +43,12 @@ namespace duneuro
{
}
void setSourceModel(const Dune::ParameterTree& config, DataTree dataTree = DataTree())
void setSourceModel(const Dune::ParameterTree& config, const Dune::ParameterTree& solverConfig,
DataTree dataTree = DataTree())
{
denseSourceModel_ = SMF::template createDense<typename Traits::RangeDOFVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config,
solverConfig);
}
void bind(const typename Traits::DipoleType& dipole, DataTree dataTree)
......
......@@ -43,17 +43,20 @@ namespace duneuro
{
}
void setSourceModel(const Dune::ParameterTree& config, DataTree dataTree = DataTree())
void setSourceModel(const Dune::ParameterTree& config, const Dune::ParameterTree& solverConfig,
DataTree dataTree = DataTree())
{
sparseSourceModel_.reset();
denseSourceModel_.reset();
density_ = source_model_default_density(config);
if (density_ == VectorDensity::sparse) {
sparseSourceModel_ = SMF::template createSparse<typename Traits::SparseRHSVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config,
solverConfig);
} else {
denseSourceModel_ = SMF::template createDense<typename Traits::DenseRHSVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config,
solverConfig);
}
}
......
......@@ -100,6 +100,7 @@ namespace duneuro
explicit UnfittedMEEGDriver(UnfittedMEEGDriverData<dim> data, const Dune::ParameterTree& config)
: data_(data)
, config_(config)
, grid_(make_structured_grid<dim>(config.sub("volume_conductor.grid")))
, fundamentalGridView_(grid_->levelGridView(0))
, levelSetGridView_(grid_->levelGridView(grid_->maxLevel()))
......@@ -123,7 +124,7 @@ namespace duneuro
Function& solution, const Dune::ParameterTree& config,
DataTree dataTree = DataTree()) override
{
eegForwardSolver_.setSourceModel(config.sub("source_model"));
eegForwardSolver_.setSourceModel(config.sub("source_model"), config_.sub("solver"));
eegForwardSolver_.bind(dipole, dataTree);
#if HAVE_TBB
eegForwardSolver_.solve(solverBackend_.local().get(),
......@@ -272,7 +273,7 @@ namespace duneuro
[&](const tbb::blocked_range<std::size_t>& range) {
User myUser(subTriangulation_, solver_, elementSearch_,
config.sub("solver"));
myUser.setSourceModel(config.sub("source_model"));
myUser.setSourceModel(config.sub("source_model"), config_.sub("solver"));
for (std::size_t index = range.begin(); index != range.end(); ++index) {
auto dt = dataTree.sub("dipole_" + std::to_string(index));
myUser.bind(dipoles[index], dt);
......@@ -288,7 +289,7 @@ namespace duneuro
});
#else
User myUser(subTriangulation_, solver_, elementSearch_, config.sub("solver"));
myUser.setSourceModel(config.sub("source_model"));
myUser.setSourceModel(config.sub("source_model"), config_.sub("solver"));
for (std::size_t index = 0; index < dipoles.size(); ++index) {
auto dt = dataTree.sub("dipole_" + std::to_string(index));
myUser.bind(dipoles[index], dt);
......@@ -323,7 +324,7 @@ namespace duneuro
[&](const tbb::blocked_range<std::size_t>& range) {
User myUser(subTriangulation_, solver_, elementSearch_,
config.sub("solver"));
myUser.setSourceModel(config.sub("source_model"));
myUser.setSourceModel(config.sub("source_model"), config_.sub("solver"));
for (std::size_t index = range.begin(); index != range.end(); ++index) {
auto dt = dataTree.sub("dipole_" + std::to_string(index));
myUser.bind(dipoles[index], dt);
......@@ -332,7 +333,7 @@ namespace duneuro