Commit 2163484a authored by Andreas Nüßing's avatar Andreas Nüßing

Merge branch 'feature/udg-make-compartments-dynamic' into 'master'

[UDG] make dipole and electrode compartment dynamic

See merge request duneuro/duneuro!43
parents d06dd4ea be41926b
......@@ -52,9 +52,8 @@ namespace duneuro
void setSourceModel(const Dune::ParameterTree& config, DataTree dataTree = DataTree())
{
denseSourceModel_ =
UDGSourceModelFactory::template createDense<compartments - 1,
typename Traits::RangeDOFVector>(
*solver_, subTriangulation_, search_, config);
UDGSourceModelFactory::template createDense<typename Traits::RangeDOFVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
}
void bind(const typename Traits::DipoleType& dipole)
......
......@@ -16,7 +16,7 @@
namespace duneuro
{
template <class GFS, int child, class ST, class V>
template <class GFS, class ST, class V>
class UDGPartialIntegrationSourceModel
: public SourceModelBase<typename GFS::Traits::GridViewType, V>
{
......@@ -34,26 +34,31 @@ namespace duneuro
using UCache = Dune::PDELab::LFSIndexCache<ULFS>;
UDGPartialIntegrationSourceModel(const GFS& gfs, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<typename BaseT::SearchType> search)
: BaseT(search), subTriangulation_(subTriangulation), ulfs_(gfs), ucache_(ulfs_)
std::shared_ptr<typename BaseT::SearchType> search,
std::size_t child)
: BaseT(search)
, subTriangulation_(subTriangulation)
, child_(child)
, ulfs_(gfs)
, ucache_(ulfs_)
{
}
virtual void assembleRightHandSide(VectorType& vector) const
{
using ChildLFS = typename ULFS::template Child<child>::Type;
using ChildLFS = typename ULFS::template Child<0>::Type;
using FESwitch =
Dune::FiniteElementInterfaceSwitch<typename ChildLFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
ChildLFS& childLfs = ulfs_.child(child);
ChildLFS& childLfs = ulfs_.child(child_);
UST ust(subTriangulation_->gridView(), *subTriangulation_);
ust.create(this->dipoleElement());
bool foundCompartment = false;
for (const auto& ep : ust) {
if (ep.domainIndex() != child)
if (ep.domainIndex() != child_)
continue;
foundCompartment = true;
ulfs_.bind(ep.subEntity(), true);
......@@ -98,12 +103,14 @@ namespace duneuro
if (!foundCompartment) {
DUNE_THROW(Dune::Exception,
"dipole should be in compartment "
<< child << " but no such compartment was found in the fundamental element");
<< child_
<< " but no such compartment was found in the fundamental element");
}
}
private:
std::shared_ptr<ST> subTriangulation_;
std::size_t child_;
mutable ULFS ulfs_;
mutable UCache ucache_;
};
......
......@@ -21,7 +21,7 @@
namespace duneuro
{
template <class GFS, int child, class ST, class V>
template <class GFS, class ST, class V>
class UDGPatchBasedVenantSourceModel
: public SourceModelBase<typename GFS::Traits::GridViewType, V>
{
......@@ -39,9 +39,10 @@ namespace duneuro
UDGPatchBasedVenantSourceModel(const GFS& gfs, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<typename BaseT::SearchType> search,
const Dune::ParameterTree& params)
std::size_t child, const Dune::ParameterTree& params)
: BaseT(search)
, subTriangulation_(subTriangulation)
, child_(child)
, elementNeighborhoodMap_(std::make_shared<ElementNeighborhoodMap<GV>>(gfs.gridView()))
, gfs_(gfs)
, monopolarVenant_(params)
......@@ -56,7 +57,7 @@ namespace duneuro
// select monopoles within each element based on a quadrature rule
std::vector<Dune::FieldVector<Real, dim>> positions;
for (const auto& element : elements) {
const auto& bbox = subTriangulation_->BBox(child, element);
const auto& bbox = subTriangulation_->BBox(child_, element);
const auto& rule =
Dune::QuadratureRules<Real, dim>::rule(bbox.type(), quadratureRuleOrder_);
for (const auto& qp : rule) {
......@@ -75,14 +76,14 @@ namespace duneuro
std::size_t offset = 0;
using ULFS = Dune::PDELab::UnfittedLocalFunctionSpace<GFS>;
using UCache = Dune::PDELab::LFSIndexCache<ULFS>;
using ChildLFS = typename ULFS::template Child<child>::Type;
using ChildLFS = typename ULFS::template Child<0>::Type;
using FESwitch =
Dune::FiniteElementInterfaceSwitch<typename ChildLFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
ULFS lfs(gfs_);
UCache cache(lfs);
ChildLFS& childLfs = lfs.child(child);
ChildLFS& childLfs = lfs.child(child_);
using RangeType = typename BasisSwitch::Range;
......@@ -91,7 +92,7 @@ namespace duneuro
const auto& element = elements[i];
ust.create(element);
for (const auto& ep : ust) {
if (ep.domainIndex() != child)
if (ep.domainIndex() != child_)
continue;
lfs.bind(ep.subEntity(), true);
cache.update();
......@@ -116,13 +117,14 @@ namespace duneuro
{
auto elementPatch =
make_element_patch(subTriangulation_, elementNeighborhoodMap_, this->elementSearch(),
this->dipole().position(), child, config_);
this->dipole().position(), child_, config_);
interpolate(elementPatch->elements(), this->dipole(), vector);
}
private:
std::shared_ptr<ST> subTriangulation_;
std::size_t child_;
std::shared_ptr<ElementNeighborhoodMap<GV>> elementNeighborhoodMap_;
const GFS& gfs_;
MonopolarVenant<Real, dim> monopolarVenant_;
......
......@@ -12,46 +12,46 @@
namespace duneuro
{
struct UDGSourceModelFactory {
template <int dipoleCompartment, class Vector, class Solver, class ST>
template <class Vector, class Solver, class ST>
static std::unique_ptr<SourceModelInterface<typename Solver::Traits::RangeField,
Solver::Traits::dimension, Vector>>
createDense(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
return Dune::Std::make_unique<UDGPartialIntegrationSourceModel<
typename Solver::Traits::FunctionSpace::GFS, dipoleCompartment, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search);
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment);
} else if (type == "patch_based_venant") {
return Dune::Std::make_unique<UDGPatchBasedVenantSourceModel<
typename Solver::Traits::FunctionSpace::GFS, dipoleCompartment, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, config);
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, config);
} else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\"");
}
}
template <int dipoleCompartment, class Vector, class Solver, class ST>
template <class Vector, class Solver, class ST>
static std::unique_ptr<SourceModelInterface<typename Solver::Traits::RangeField,
Solver::Traits::dimension, Vector>>
createSparse(
const Solver& solver, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<KDTreeElementSearch<typename Solver::Traits::FundamentalGridView>> search,
const Dune::ParameterTree& config)
std::size_t dipoleCompartment, const Dune::ParameterTree& config)
{
const auto type = config.get<std::string>("type");
if (type == "partial_integration") {
return Dune::Std::make_unique<UDGPartialIntegrationSourceModel<
typename Solver::Traits::FunctionSpace::GFS, dipoleCompartment, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search);
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment);
} else if (type == "patch_based_venant") {
return Dune::Std::make_unique<UDGPatchBasedVenantSourceModel<
typename Solver::Traits::FunctionSpace::GFS, dipoleCompartment, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, config);
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, config);
} else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\"");
......
......@@ -13,7 +13,7 @@
namespace duneuro
{
template <class GFS, int child, class ST>
template <class GFS, class ST>
class UDGTransferMatrixRHS
{
public:
......@@ -27,7 +27,8 @@ namespace duneuro
using ULFS = Dune::PDELab::UnfittedLocalFunctionSpace<GFS>;
using UCache = Dune::PDELab::LFSIndexCache<ULFS>;
UDGTransferMatrixRHS(const GFS& gfs, std::shared_ptr<ST> st) : st_(st), gfs_(gfs)
UDGTransferMatrixRHS(const GFS& gfs, std::shared_ptr<ST> st, std::size_t child)
: st_(st), gfs_(gfs), child_(child)
{
}
......@@ -35,7 +36,7 @@ namespace duneuro
const Element& electrodeElement, const Coordinate& electrodeLocal,
DOFVector& output)
{
using ChildLFS = typename ULFS::template Child<child>::Type;
using ChildLFS = typename ULFS::template Child<0>::Type;
using FESwitch =
Dune::FiniteElementInterfaceSwitch<typename ChildLFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
......@@ -46,7 +47,7 @@ namespace duneuro
ULFS ulfs(gfs_);
UCache ucache(ulfs);
ChildLFS& childLfs = ulfs.child(child);
ChildLFS& childLfs = ulfs.child(child_);
ust.create(electrodeElement);
if (ust.begin() == ust.end()) {
......@@ -54,7 +55,7 @@ namespace duneuro
}
for (const auto& ep : ust) {
if (ep.domainIndex() != child)
if (ep.domainIndex() != child_)
continue;
ulfs.bind(ep.subEntity(), true);
ucache.update();
......@@ -78,7 +79,7 @@ namespace duneuro
ust.create(referenceElement);
for (const auto& ep : ust) {
if (ep.domainIndex() != child)
if (ep.domainIndex() != child_)
continue;
ulfs.bind(ep.subEntity(), true);
ucache.update();
......@@ -101,6 +102,7 @@ namespace duneuro
private:
std::shared_ptr<ST> st_;
const GFS& gfs_;
std::size_t child_;
};
}
......
......@@ -120,9 +120,9 @@ namespace duneuro
Dune::Timer timer;
rightHandSideVector = 0.0;
// assemble right hand side
UDGTransferMatrixRHS<typename Traits::FunctionSpace::GFS, 0,
typename Traits::SubTriangulation>
rhsAssembler(solver_->functionSpace().getGFS(), subTriangulation_);
UDGTransferMatrixRHS<typename Traits::FunctionSpace::GFS, typename Traits::SubTriangulation>
rhsAssembler(solver_->functionSpace().getGFS(), subTriangulation_,
config.get<std::size_t>("compartment"));
rhsAssembler.assembleRightHandSide(reference.element, reference.localPosition,
electrode.element, electrode.localPosition,
rightHandSideVector);
......
......@@ -55,14 +55,14 @@ namespace duneuro
density_ = source_model_default_density(config);
if (density_ == VectorDensity::sparse) {
sparseSourceModel_ =
UDGSourceModelFactory::template createSparse<compartments - 1,
typename Traits::SparseRHSVector>(
*solver_, subTriangulation_, search_, config);
UDGSourceModelFactory::template createSparse<typename Traits::SparseRHSVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"),
config);
} else {
denseSourceModel_ =
UDGSourceModelFactory::template createDense<compartments - 1,
typename Traits::DenseRHSVector>(
*solver_, subTriangulation_, search_, config);
UDGSourceModelFactory::template createDense<typename Traits::DenseRHSVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"),
config);
}
}
......
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