Commit 68a86a38 authored by Andreas Nüßing's avatar Andreas Nüßing

[SubtractionUDG] add first implementation

parent a9287a4b
......@@ -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_)
......
......@@ -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>
......
......@@ -54,7 +54,7 @@ namespace duneuro
const int dim = IG::dimension;
/** intorder **/
const int intorder = intorderadd_lb + 2 * lfsv.finiteElement().localBasis().order();
const int intorder = intorderadd_lb + 2 * FESwitch::basis(lfsv.finiteElement()).order();
const auto& geometryInInside = ig.geometryInInside();
const auto& geometry = ig.geometry();
......@@ -86,16 +86,11 @@ namespace duneuro
void lambda_volume(const EG& eg, const LFSV& lfsv, R& r) const
{
/** domain and range field type **/
typedef
typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType
DF;
typedef
typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType
RF;
typedef typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType
RangeType;
typedef typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType
JacobianType;
typedef Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType> FESwitch;
typedef Dune::BasisInterfaceSwitch<typename FESwitch::Basis> BasisSwitch;
typedef typename BasisSwitch::DomainField DF;
typedef typename BasisSwitch::RangeField RF;
typedef typename BasisSwitch::Range RangeType;
typedef typename LFSV::Traits::SizeType size_type;
......@@ -105,33 +100,27 @@ namespace duneuro
const auto& geometry = eg.geometry();
const int intorder = intorderadd
+ 2 * std::max(lfsv.finiteElement().localBasis().order(),
lfsv.finiteElement().localBasis().order());
const int intorder = intorderadd + 2 * FESwitch::basis(lfsv.finiteElement()).order();
/** select quadrature rule**/
Dune::GeometryType gt = geometry.type();
const auto& rule = Dune::QuadratureRules<DF, dim>::rule(gt, intorder);
typename PROBLEMDATA::Traits::PermTensorType sigma_corr(param.sigma_corr(eg.entity()));
const auto& localcenter = Dune::ReferenceElements<DF, dim>::general(gt).position(0, 0);
typename PROBLEMDATA::Traits::PermTensorType sigma_corr = param.A(eg, localcenter);
sigma_corr -= param.get_sigma_infty();
std::vector<RangeType> phi(lfsv.size());
std::vector<JacobianType> js(lfsv.size());
Dune::FieldMatrix<DF, dimw, dim> jac;
std::vector<Dune::FieldVector<RF, dim>> gradphi(lfsv.size());
std::vector<Dune::FieldMatrix<RF, 1, dim>> gradphi(lfsv.size());
/** loop over quadrature points **/
for (const auto& qp : rule) {
/* evaluate shape functions */
lfsv.finiteElement().localBasis().evaluateFunction(qp.position(), phi);
FESwitch::basis(lfsv.finiteElement()).evaluateFunction(qp.position(), phi);
/* evaluate gradient of basis functions on reference element */
lfsv.finiteElement().localBasis().evaluateJacobian(qp.position(), js);
/* transform gradients from reference element to real element */
jac = geometry.jacobianInverseTransposed(qp.position());
for (size_type i = 0; i < lfsv.size(); i++)
jac.mv(js[i][0], gradphi[i]);
BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()), geometry, qp.position(),
gradphi);
/* get global coordinates of the quadrature point */
typename PROBLEMDATA::Traits::DomainType x = geometry.global(qp.position());
......@@ -152,7 +141,7 @@ namespace duneuro
DUNE_THROW(Dune::Exception, "lambda_volume f nan");
}
for (size_type i = 0; i < lfsv.size(); i++)
r.accumulate(lfsv, i, (f * gradphi[i] * factor));
r.accumulate(lfsv, i, (f * gradphi[i][0] * factor));
}
}
......@@ -161,14 +150,11 @@ namespace duneuro
void lambda_skeleton(const IG& ig, const LFSV& lfsv_s, const LFSV& lfsv_n, R& r_s, R& r_n) const
{
/** domain and range field type **/
typedef
typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::DomainFieldType
DF;
typedef
typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType
RF;
typedef typename LFSV::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType
RangeType;
typedef Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType> FESwitch;
typedef Dune::BasisInterfaceSwitch<typename FESwitch::Basis> BasisSwitch;
typedef typename BasisSwitch::DomainField DF;
typedef typename BasisSwitch::RangeField RF;
typedef typename BasisSwitch::Range RangeType;
typedef typename LFSV::Traits::SizeType size_type;
/** dimension **/
......@@ -176,8 +162,8 @@ namespace duneuro
/** intorder **/
const int intorder = intorderadd
+ 2 * std::max(lfsv_s.finiteElement().localBasis().order(),
lfsv_n.finiteElement().localBasis().order());
+ 2 * std::max(FESwitch::basis(lfsv_s.finiteElement()).order(),
FESwitch::basis(lfsv_n.finiteElement()).order());
const auto& geometryInInside = ig.geometryInInside();
const auto& geometryInOutside = ig.geometryInOutside();
......@@ -188,18 +174,13 @@ namespace duneuro
const Dune::QuadratureRule<DF, dim - 1>& rule =
Dune::QuadratureRules<DF, dim - 1>::rule(gtface, intorder);
const auto& insideEntity = ig.inside();
const auto& outsideEntity = ig.outside();
/** compute weights **/
/* evaluate permability tensor */
const Dune::FieldVector<DF, dim>& inside_local =
Dune::ReferenceElements<DF, dim>::general(insideEntity.type()).position(0, 0);
const Dune::FieldVector<DF, dim>& outside_local =
Dune::ReferenceElements<DF, dim>::general(outsideEntity.type()).position(0, 0);
typename PROBLEMDATA::Traits::PermTensorType A_s, A_n;
A_s = param.A(insideEntity, inside_local);
A_n = param.A(outsideEntity, outside_local);
const Dune::FieldVector<DF, dim - 1>& localcenter =
Dune::ReferenceElements<DF, dim - 1>::general(ig.geometry().type()).position(0, 0);
A_s = param.A(ig, localcenter, ConvectionDiffusion_DG_Side::inside);
A_n = param.A(ig, localcenter, ConvectionDiffusion_DG_Side::outside);
/* tensor times normal */
auto weights = weighting(ig, A_s, A_n);
......@@ -216,13 +197,15 @@ namespace duneuro
Dune::FieldVector<DF, dim> iplocal_n = geometryInOutside.global(qp.position());
/* evaluate basis functions */
lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s, phi_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(iplocal_n, phi_n);
FESwitch::basis(lfsv_s.finiteElement()).evaluateFunction(iplocal_s, phi_s);
FESwitch::basis(lfsv_n.finiteElement()).evaluateFunction(iplocal_n, phi_n);
/* get sigma^corr for both sides of the interface */
typename PROBLEMDATA::Traits::PermTensorType sigma_corr_s, sigma_corr_n;
sigma_corr_s = param.sigma_corr(insideEntity);
sigma_corr_n = param.sigma_corr(outsideEntity);
sigma_corr_s = param.A(ig, localcenter, ConvectionDiffusion_DG_Side::inside);
sigma_corr_s -= param.get_sigma_infty();
sigma_corr_n = param.A(ig, localcenter, ConvectionDiffusion_DG_Side::outside);
sigma_corr_n -= param.get_sigma_infty();
const auto& global = geometry.global(qp.position());
......
#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/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 LOP = SubtractionDG<Problem, EdgeNormProvider, 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)
, lop_(problem_,
solverConfig.get<bool>("weights", true) ? ConvectionDiffusion_DG_Weights::weightsOn :
ConvectionDiffusion_DG_Weights::weightsOff,
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_;
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);
}
}
......
......@@ -123,7 +123,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 +272,7 @@ namespace duneuro
[&](const tbb::blocked_range<std::size_t>& range) {
User myUser(subTriangulation_, solver_, elementSearch_,
config.sub("solver"));