...
 
Commits (6)
......@@ -34,15 +34,15 @@ namespace duneuro
Dune::PDELab::MeshType::conforming, Dune::SolverCategory::sequential>;
};
template <class VC, ElementType elementType, unsigned int degree, class DF = double,
class RF = double, class JF = double>
template <class VC, ElementType elementType, unsigned int degree, class P, class DF = double,
class RF = double , class JF = double>
struct CGSolverTraits {
static const int dimension = VC::dim;
using VolumeConductor = VC;
using GridView = typename VC::GridView;
using CoordinateFieldType = typename VC::ctype;
using ElementSearch = KDTreeElementSearch<GridView>;
using Problem = ConvectionDiffusionCGDefaultParameter<VC>;
using Problem = P;
using DirichletExtension = Dune::PDELab::ConvectionDiffusionDirichletExtensionAdapter<Problem>;
using BoundaryCondition = Dune::PDELab::ConvectionDiffusionBoundaryConditionAdapter<Problem>;
using FunctionSpace =
......@@ -56,34 +56,46 @@ namespace duneuro
LinearProblemSolver<typename Assembler::GO, DomainDOFVector, RangeDOFVector>;
};
template <class VC, ElementType elementType, unsigned int degree, class DF = double,
template <class VC, ElementType elementType, unsigned int degree,
class P = ConvectionDiffusionCGDefaultParameter<VC>, class DF = double,
class RF = double, class JF = double>
class CGSolver
{
public:
using Traits = CGSolverTraits<VC, elementType, degree, DF, RF, JF>;
using Traits = CGSolverTraits<VC, elementType, degree, P, DF, RF, JF>;
CGSolver(std::shared_ptr<const VC> volumeConductor,
std::shared_ptr<const typename Traits::ElementSearch> search,
std::shared_ptr<typename Traits::Problem> problem,
const Dune::ParameterTree& config, DataTree dataTree = DataTree())
: volumeConductor_(volumeConductor)
, search_(search)
, problem_(volumeConductor)
, dirichletExtension_(volumeConductor->gridView(), problem_)
, boundaryCondition_(volumeConductor->gridView(), problem_)
, problem_(problem)
, dirichletExtension_(volumeConductor->gridView(),*problem_)
, boundaryCondition_(volumeConductor->gridView(), *problem_)
, 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 ?
(1 << VC::dim) + 1 :
Dune::StaticPower<3, VC::dim>::power)
, linearSolver_(assembler_.getGO(), config)
{
assert(volumeConductor_);
dataTree.set("degree", degree);
dataTree.set("element_type", to_string(elementType));
functionSpace_.assembleConstraints(boundaryCondition_);
}
CGSolver(std::shared_ptr<const VC> volumeConductor,
std::shared_ptr<const typename Traits::ElementSearch> search,
const Dune::ParameterTree& config, DataTree dataTree = DataTree())
: CGSolver(volumeConductor, search,
std::make_shared<typename Traits::Problem>(volumeConductor), config, dataTree)
{
}
template <typename SolverBackend>
void solve(SolverBackend& solverBackend, const typename Traits::RangeDOFVector& rightHandSide,
typename Traits::DomainDOFVector& solution, const Dune::ParameterTree& config,
......@@ -94,6 +106,19 @@ namespace duneuro
linearSolver_.apply(solverBackend, solution, rightHandSide, config, dataTree);
dataTree.set("time", timer.elapsed());
}
template <typename SolverBackend>
void solve(SolverBackend& solverBackend, typename Traits::DomainDOFVector& solution,
const Dune::ParameterTree& config, DataTree dataTree = DataTree())
{
Dune::Timer timer;
//randomize_uniform(Dune::PDELab::Backend::native(solution), DF(-1.0), DF(1.0));
Dune::PDELab::interpolate(dirichletExtension_, functionSpace_.getGFS(), solution);
functionSpace_.setNonConstrainedDOFS(solution, 0.0);
linearSolver_.apply(solverBackend, solution, config, dataTree);
dataTree.set("time", timer.elapsed());
}
const typename Traits::FunctionSpace& functionSpace() const
{
......@@ -111,9 +136,9 @@ namespace duneuro
}
private:
std::shared_ptr<const typename Traits::VolumeConductor> volumeConductor_;
std::shared_ptr<const VC> volumeConductor_;
std::shared_ptr<const typename Traits::ElementSearch> search_;
typename Traits::Problem problem_;
std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::DirichletExtension dirichletExtension_;
typename Traits::BoundaryCondition boundaryCondition_;
typename Traits::FunctionSpace functionSpace_;
......
......@@ -113,6 +113,24 @@ namespace duneuro
}
}
}
template <class T>
void set_matrix_row(DenseMatrix<T>& matrix, std::size_t row,
const std::vector<T>& vector)
{
if (row >= matrix.rows()) {
DUNE_THROW(Dune::Exception, "tried to set row " << row << " but only " << matrix.rows()
<< " rows are present");
}
if (vector.size() != matrix.cols()) {
DUNE_THROW(Dune::Exception, "tried to set row with " << vector.size()
<< " entries, but row has actually "
<< matrix.cols() << " entries");
}
for (unsigned int i = 0; i < vector.size(); ++i) {
matrix(row, i) = vector[i];
}
}
}
#endif // DUNEURO_MATRIX_UTILITIES_HH
......@@ -4,6 +4,8 @@
#include <dune/common/parametertree.hh>
#include <duneuro/common/default_grids.hh>
#include <duneuro/common/cg_solver.hh>
#include <duneuro/common/cg_solver_backend.hh>
#include <duneuro/common/dg_solver.hh>
#include <duneuro/common/dg_solver_backend.hh>
#include <duneuro/common/fitted_driver_data.hh>
......@@ -14,6 +16,7 @@
#include <duneuro/io/volume_conductor_reader.hh>
#include <duneuro/io/vtk_writer.hh>
#include <duneuro/tes/tdcs_driver_interface.hh>
#include <duneuro/tes/tdcs_patch_cg_parameter.hh>
#include <duneuro/tes/tdcs_patch_dg_parameter.hh>
#if HAVE_DUNE_SUBGRID
#include <duneuro/common/geometry_adaption.hh>
......@@ -24,6 +27,16 @@ namespace duneuro
template <FittedSolverType solverType, class VC, ElementType et, int degree>
struct TDCSSelectFittedSolver;
template <class VC, ElementType et, int degree>
struct TDCSSelectFittedSolver<FittedSolverType::cg, VC, et, degree> {
using ProblemType = TDCSPatchCGParameter<VC>;
using SolverType = CGSolver<VC, et, degree, ProblemType>;
using SolverBackendType = CGSolverBackend<SolverType, et>;
};
template <FittedSolverType solverType, class VC, ElementType et, int degree>
struct TDCSSelectFittedSolver;
template <class VC, ElementType et, int degree>
struct TDCSSelectFittedSolver<FittedSolverType::dg, VC, et, degree> {
using ProblemType = TDCSPatchDGParameter<VC>;
......
......@@ -3,6 +3,10 @@
#include <duneuro/tes/fitted_tdcs_driver.hh>
template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::hexahedron,
......
......@@ -3,11 +3,17 @@
#include <duneuro/tes/fitted_tdcs_driver.hh>
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1>;
#if HAVE_DUNE_SUBGRID
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, true>;
template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, true>;
#endif
This diff is collapsed.
#include <config.h>
#include <duneuro/tes/fitted_tdcs_point_driver.hh>
template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1>;
#include <config.h>
#include <duneuro/tes/fitted_tdcs_point_driver.hh>
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1>;
#if HAVE_DUNE_SUBGRID
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, true>;
template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, true>;
#endif
#ifndef DUNEURO_TDCS_POINT_RHS_HH
#define DUNEURO_TDCS_POINT_RHS_HH
#include <dune/common/fvector.hh>
#include <dune/pdelab/backend/interface.hh>
#include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
#include <duneuro/tes/tdcs_rhs_interface.hh>
namespace duneuro
{
template <class GFS, class V>
class FittedTdcsRHS
: public TdcsRHSInterface<typename GFS::Traits::GridViewType, V>
{
public:
using GV = typename GFS::Traits::GridViewType;
enum { dim = GV::dimension };
using Real = typename GV::ctype;
using BaseT = TdcsRHSInterface<GV, V>;
using LocalCoordinate = typename BaseT::LocalCoordinate;
using Element = typename BaseT::Element;
using VectorType = typename BaseT::VectorType;
using LFS = Dune::PDELab::LocalFunctionSpace<GFS>;
using Cache = Dune::PDELab::LFSIndexCache<LFS>;
using FESwitch = Dune::FiniteElementInterfaceSwitch<typename LFS::Traits::FiniteElementType>;
using BasisSwitch = Dune::BasisInterfaceSwitch<typename FESwitch::Basis>;
using RangeType = typename BasisSwitch::Range;
FittedTdcsRHS(const GFS& gfs)
: gfs_(gfs)
, lfsReference_(gfs_)
, cacheReference_(lfsReference_)
, lfsElectrode_(gfs_)
, cacheElectrode_(lfsElectrode_)
{
}
virtual void bind(const Element& referenceElement, const LocalCoordinate& referenceLocal,
const Element& electrodeElement,
const LocalCoordinate& electrodeLocal) override
{
bind(lfsReference_, cacheReference_, phiReference_, referenceElement, referenceLocal);
bind(lfsElectrode_, cacheElectrode_, phiElectrode_, electrodeElement, electrodeLocal);
}
virtual void assembleRightHandSide(VectorType& output) const override
{
for (unsigned int i = 0; i < cacheElectrode_.size(); ++i) {
output[cacheElectrode_.containerIndex(i)] = phiElectrode_[i];
}
for (unsigned int i = 0; i < cacheReference_.size(); ++i) {
output[cacheReference_.containerIndex(i)] -= phiReference_[i];
}
}
private:
const GFS& gfs_;
LFS lfsReference_;
Cache cacheReference_;
std::vector<RangeType> phiReference_;
LFS lfsElectrode_;
Cache cacheElectrode_;
std::vector<RangeType> phiElectrode_;
void bind(LFS& lfs, Cache& cache, std::vector<RangeType>& phi, const Element& element,
const LocalCoordinate& local)
{
lfs.bind(element);
cache.update();
phi.resize(lfs.size());
FESwitch::basis(lfs.finiteElement()).evaluateFunction(local, phi);
}
};
}
#endif // DUNEURO_TDCS_POINT_RHS_HH
#ifndef DUNEURO_FITTED_TDCS_RHS_FACTORY_HH
#define DUNEURO_FITTED_TDCS_RHS_FACTORY_HH
#include <dune/common/parametertree.hh>
#include <duneuro/tes/fitted_tdcs_point_rhs.hh>
#include <duneuro/tes/tdcs_rhs_interface.hh>
namespace duneuro
{
struct FittedTdcsRHSFactory {
template <class Vector, class Solver>
static std::unique_ptr<TdcsRHSInterface<
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<FittedTdcsRHS<
typename Solver::Traits::FunctionSpace::GFS, Vector>>(solver.functionSpace().getGFS());
} else {
DUNE_THROW(Dune::Exception, "unknown transfer matrix type \"" << type << "\"");
}
}
};
}
#endif // DUNEURO_FITTED_TDCS_RHS_FACTORY_HH
......@@ -3,6 +3,14 @@
#include <duneuro/tes/udg_tdcs_driver.hh>
#endif
extern template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
extern template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1>;
extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1>;
extern template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1>;
extern template class duneuro::FittedTDCSDriver<2, duneuro::ElementType::hexahedron,
......@@ -12,6 +20,9 @@ extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::tetrahe
extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1>;
#if HAVE_DUNE_SUBGRID
extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, true>;
#
extern template class duneuro::FittedTDCSDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, true>;
#endif
......@@ -44,7 +55,30 @@ namespace duneuro
if (type == "fitted") {
auto solverType = config.get<std::string>("solver_type");
auto elementType = config.get<std::string>("element_type");
if (solverType == "dg") {
if (solverType == "cg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSDriver<3, ElementType::tetrahedron,
FittedSolverType::cg, 1>>(
data.fittedData, patchSet, config, dataTree);
} else if (elementType == "hexahedron") {
auto geometryAdapted = config.get<bool>("geometry_adapted", false);
if (geometryAdapted) {
#if HAVE_DUNE_SUBGRID
return Dune::Std::make_unique<FittedTDCSDriver<3, ElementType::hexahedron,
FittedSolverType::cg, 1, true>>(
data.fittedData, patchSet, config, dataTree);
#else
DUNE_THROW(Dune::Exception, "geometry adaption needs dune-subgrid");
#endif
} else {
return Dune::Std::make_unique<FittedTDCSDriver<3, ElementType::hexahedron,
FittedSolverType::cg, 1, false>>(
data.fittedData, patchSet, config, dataTree);
}
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "dg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSDriver<3, ElementType::tetrahedron,
FittedSolverType::dg, 1>>(
......@@ -104,7 +138,19 @@ namespace duneuro
if (type == "fitted") {
auto solverType = config.get<std::string>("solver_type");
auto elementType = config.get<std::string>("element_type");
if (solverType == "dg") {
if (solverType == "cg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSDriver<2, ElementType::tetrahedron,
FittedSolverType::cg, 1>>(
data.fittedData, patchSet, config, dataTree);
} else if (elementType == "hexahedron") {
return Dune::Std::make_unique<FittedTDCSDriver<2, ElementType::hexahedron,
FittedSolverType::cg, 1, false>>(
data.fittedData, patchSet, config, dataTree);
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "dg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSDriver<2, ElementType::tetrahedron,
FittedSolverType::dg, 1>>(
......
#ifndef DUNEURO_TDCS_MATRIX_SOLVER_HH
#define DUNEURO_TDCS_MATRIX_SOLVER_HH
#include <dune/common/parametertree.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh>
#include <duneuro/common/make_dof_vector.hh>
#include <duneuro/eeg/electrode_projection_interface.hh>
#include <duneuro/io/data_tree.hh>
namespace duneuro {
template <class S> struct TDCSMatrixSolverTraits {
using Solver = S;
static const unsigned int dimension = Solver::Traits::dimension;
using FunctionSpace = typename Solver::Traits::FunctionSpace;
using DomainDOFVector = typename Solver::Traits::DomainDOFVector;
using RangeDOFVector = typename Solver::Traits::RangeDOFVector;
using CoordinateFieldType = typename Solver::Traits::CoordinateFieldType;
using Coordinate = Dune::FieldVector<CoordinateFieldType, dimension>;
using ProjectedPosition =
duneuro::ProjectedElectrode<typename Solver::Traits::GridView>;
};
template <class S, class RHSFactory> class TDCSMatrixSolver {
public:
using Traits = TDCSMatrixSolverTraits<S>;
TDCSMatrixSolver(std::shared_ptr<typename Traits::Solver> solver,
const Dune::ParameterTree &config)
: solver_(solver),
rightHandSideVector_(solver_->functionSpace().getGFS(), 0.0),
config_(config) {}
template <class SolverBackend>
std::unique_ptr<DenseMatrix<double>>
solve(SolverBackend &solverBackend,
const ElectrodeProjectionInterface<
typename Traits::Solver::Traits::GridView> &projectedElectrodes,
const Dune::ParameterTree &config, DataTree dataTree = DataTree()) {
using GV = typename S::Traits::FunctionSpace::GFS::Traits::GridViewType;
auto tdcsMatrix = Dune::Std::make_unique<DenseMatrix<double>>(
projectedElectrodes.size(),
GV::dimension * solver_->volumeConductor()->gridView().size(0));
auto solver_config = config.sub("solver");
typename Traits::DomainDOFVector solution(solver_->functionSpace().getGFS(),
0.0);
for (std::size_t index = 1; index < projectedElectrodes.size(); ++index) {
solve(solverBackend.get(), projectedElectrodes.getProjection(index),
projectedElectrodes.getProjection(0), solution,
rightHandSideVector_, solver_config,
dataTree.sub("solver.electrode_" + std::to_string(index)));
auto result = evaluateElectricField(*solver_, solution);
set_matrix_row(*tdcsMatrix, index, result);
}
return tdcsMatrix;
}
private:
std::shared_ptr<typename Traits::Solver> solver_;
typename Traits::RangeDOFVector rightHandSideVector_;
Dune::ParameterTree config_;
template <class SolverBackend>
void solve(SolverBackend &solverBackend,
const typename Traits::ProjectedPosition &reference,
const typename Traits::ProjectedPosition &electrode,
typename Traits::DomainDOFVector &solution,
typename Traits::RangeDOFVector &rightHandSideVector,
const Dune::ParameterTree &config,
DataTree dataTree = DataTree()) const {
Dune::Timer timer;
rightHandSideVector = 0.0;
// assemble right hand side
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();
// solve system
solver_->solve(solverBackend, rightHandSideVector, solution, config,
dataTree.sub("linear_system_solver"));
timer.stop();
dataTree.set("time_solution", timer.lastElapsed());
dataTree.set("time", timer.elapsed());
}
template <class Solver>
std::vector<double> evaluateElectricField(
const Solver &solver,
const typename Solver::Traits::DomainDOFVector &function) {
using GV =
typename Solver::Traits::FunctionSpace::GFS::Traits::GridViewType;
using DGFG = Dune::PDELab::DiscreteGridFunctionGradient<
typename Solver::Traits::FunctionSpace::GFS,
typename Solver::Traits::DomainDOFVector>;
DGFG dgfg(solver_->functionSpace().getGFS(), function);
// evaluate discrete grid function
std::vector<double> result(GV::dimension *
solver_->volumeConductor()->gridView().size(0));
std::size_t offset = 0;
for (const auto &element :
Dune::elements(solver_->volumeConductor()->gridView())) {
typename DGFG::Traits::RangeType y;
typename DGFG::Traits::DomainType xglobal = element.geometry().center();
typename DGFG::Traits::DomainType xlocal =
element.geometry().local(xglobal);
dgfg.evaluate(element, xlocal, y);
const auto &conductivity = solver_->volumeConductor()->tensor(element);
std::vector<double> z(GV::dimension);
for (unsigned int i = 0; i < GV::dimension; ++i) {
for (unsigned int j = 0; j < GV::dimension; ++j) {
z[i] += conductivity[i][j] * y[j];
}
result[offset + i] = z[i];
}
offset += GV::dimension;
}
return result;
}
};
}
#endif // DUNEURO_TDCS_MATRIX_SOLVER_HH
#ifndef DUNEURO_TDCS_POINT_DRIVER_FACTORY_HH
#define DUNEURO_TDCS_POINT_DRIVER_FACTORY_HH
#include <memory>
#include <dune/common/parametertree.hh>
#include <duneuro/common/fitted_driver_data.hh>
#include <duneuro/io/data_tree.hh>
#include <duneuro/tes/tdcs_point_driver_interface.hh>
namespace duneuro
{
template <int dim>
struct TDCSPointDriverData {
FittedDriverData<dim> fittedData;
};
template <int dim>
struct TDCSPointDriverFactory {
static std::unique_ptr<TDCSPointDriverInterface<dim>>
make_tdcs_point_driver(const Dune::ParameterTree& config,
const TDCSPointDriverData<dim>& data = TDCSPointDriverData<dim>(),
DataTree dataTree = DataTree());
};
}
#include <duneuro/tes/tdcs_point_driver_factory_impl.hh>
#endif // DUNEURO_TDCS_POINT_DRIVER_FACTORY_HH
#include <dune/common/std/memory.hh>
#include <duneuro/tes/fitted_tdcs_point_driver.hh>
extern template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<2, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::cg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::tetrahedron,
duneuro::FittedSolverType::dg, 1, false>;
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, false>;
#if HAVE_DUNE_SUBGRID
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::cg, 1, true>;
extern template class duneuro::FittedTDCSPointDriver<3, duneuro::ElementType::hexahedron,
duneuro::FittedSolverType::dg, 1, true>;
#endif
namespace duneuro
{
template <>
std::unique_ptr<TDCSPointDriverInterface<2>>
TDCSPointDriverFactory<2>::make_tdcs_point_driver(const Dune::ParameterTree& config,
const TDCSPointDriverData<2>& data, DataTree dataTree)
{
auto type = config.get<std::string>("type");
if (type == "fitted") {
auto solverType = config.get<std::string>("solver_type");
auto elementType = config.get<std::string>("element_type");
if (solverType == "cg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::tetrahedron,
FittedSolverType::cg, 1>>(
data.fittedData, config, dataTree);
} else if (elementType == "hexahedron") {
auto geometryAdapted = config.get<bool>("geometry_adapted", false);
if (geometryAdapted) {
#if HAVE_DUNE_SUBGRID
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::hexahedron,
FittedSolverType::cg, 1, true>>(
data.fittedData, config, dataTree);
#else
DUNE_THROW(Dune::Exception, "geometry adaption needs dune-subgrid");
#endif
} else {
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::hexahedron,
FittedSolverType::cg, 1, false>>(
data.fittedData, config, dataTree);
}
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "dg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::tetrahedron,
FittedSolverType::dg, 1>>(
data.fittedData, config, dataTree);
} else if (elementType == "hexahedron") {
auto geometryAdapted = config.get<bool>("geometry_adapted", false);
if (geometryAdapted) {
#if HAVE_DUNE_SUBGRID
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::hexahedron,
FittedSolverType::dg, 1, true>>(
data.fittedData, config, dataTree);
#else
DUNE_THROW(Dune::Exception, "geometry adaption needs dune-subgrid");
#endif
} else {
return Dune::Std::make_unique<FittedTDCSPointDriver<2, ElementType::hexahedron,
FittedSolverType::dg, 1, false>>(
data.fittedData, config, dataTree);
}
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown solver type \"" << solverType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown type \"" << type << "\"");
}
}
template <>
std::unique_ptr<TDCSPointDriverInterface<3>>
TDCSPointDriverFactory<3>::make_tdcs_point_driver(const Dune::ParameterTree& config,
const TDCSPointDriverData<3>& data, DataTree dataTree)
{
auto type = config.get<std::string>("type");
if (type == "fitted") {
auto solverType = config.get<std::string>("solver_type");
auto elementType = config.get<std::string>("element_type");
if (solverType == "cg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSPointDriver<3, ElementType::tetrahedron,
FittedSolverType::cg, 1>>(
data.fittedData, config, dataTree);
} else if (elementType == "hexahedron") {
auto geometryAdapted = config.get<bool>("geometry_adapted", false);
if (geometryAdapted) {
#if HAVE_DUNE_SUBGRID
return Dune::Std::make_unique<FittedTDCSPointDriver<3, ElementType::hexahedron,
FittedSolverType::cg, 1, true>>(
data.fittedData, config, dataTree);
#else
DUNE_THROW(Dune::Exception, "geometry adaption needs dune-subgrid");
#endif
} else {
return Dune::Std::make_unique<FittedTDCSPointDriver<3, ElementType::hexahedron,
FittedSolverType::cg, 1, false>>(
data.fittedData, config, dataTree);
}
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else if (solverType == "dg") {
if (elementType == "tetrahedron") {
return Dune::Std::make_unique<FittedTDCSPointDriver<3, ElementType::tetrahedron,
FittedSolverType::dg, 1>>(
data.fittedData, config, dataTree);
} else if (elementType == "hexahedron") {
return Dune::Std::make_unique<FittedTDCSPointDriver<3, ElementType::hexahedron,
FittedSolverType::dg, 1, false>>(
data.fittedData, config, dataTree);
} else {
DUNE_THROW(Dune::Exception, "unknown element type \"" << elementType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown solver type \"" << solverType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown type \"" << type << "\"");
}
}
}
#ifndef DUNEURO_TDCS_POINT_DRIVER_INTERFACE_HH
#define DUNEURO_TDCS_POINT_DRIVER_INTERFACE_HH
#include <vector>
#include <duneuro/common/dense_matrix.hh>
#include <duneuro/common/function.hh>
#include <duneuro/io/data_tree.hh>
namespace duneuro
{
template <int dim>
struct TDCSPointDriverInterface {
static const int dimension = dim;
using FieldType = double;
using CoordinateType = Dune::FieldVector<FieldType, dimension>;
/**
* \brief create a domain function for the given interface
*
* This domain function mainly serves as data storage, as the internal data structure is hidden
* through type erasure. It can be passed back to the driver which knows how to treat it.
*/
virtual std::unique_ptr<Function> makeDomainFunction() const = 0;
/**
* \brief solve the tdcs point forward problem for given elctrodes
*
* solve the tdcs forward problem for two point electrodes
*/
virtual void solveTDCSPointForward(Function& solution, const Dune::ParameterTree& config,
DataTree dataTree = DataTree()) = 0;
/**
* \brief set the tdcs electrodes of this driver
*
* Subsequent calls to evaluateAtElectrodes will use the first two given electrodes.
* Note that the electrodes might be projected onto the drivers domain.
*/
virtual void setElectrodes(const std::vector<CoordinateType>& electrodes,
const Dune::ParameterTree& config) = 0;
/**
* \brief write the given solution to a file
*/
virtual void write(const Function& solution, const Dune::ParameterTree& config,
DataTree dataTree = DataTree()) const = 0;
/**
* \brief write the model without a solution to a file
*/
virtual void write(const Dune::ParameterTree& config, DataTree dataTree = DataTree()) const = 0;
/**
* \brief compute the tdcs matrix
*
* compute the electric field \sigma \nabla u in each element center
* for each pair of reference electrode and every other electrode.
* Note that setElectrodes has to be called before using this method.
*/
virtual std::unique_ptr<DenseMatrix<FieldType>>
computeTdcsMatrix(const Dune::ParameterTree& config, DataTree dataTree = DataTree()) = 0;
virtual std::vector<CoordinateType> getProjectedElectrodes() const =0;
/**
* \brief write the each element center into a matrix
*
* gives the corresponding ordering of the element center used in for the tdcsMatrix
*/
virtual std::unique_ptr<DenseMatrix<FieldType>> elementCenter() = 0;
virtual ~TDCSPointDriverInterface()
{
}
};
}
#endif // DUNEURO_TDCS_POINT_DRIVER_INTERFACE_HH
#ifndef DUNEURO_TDCS_RHS_INTERFACE_HH
#define DUNEURO_TDCS_RHS_INTERFACE_HH
#include <dune/common/fvector.hh>
namespace duneuro
{
template <class GV, class V>
struct TdcsRHSInterface {
public:
using Element = typename GV::template Codim<0>::Entity;
using LocalCoordinate = Dune::FieldVector<typename GV::ctype, GV::dimension>;
using VectorType = V;
virtual void bind(const Element& referenceElement, const LocalCoordinate& referenceLocal,
const Element& electrodeElement, const LocalCoordinate& electrodeLocal) = 0;
virtual void assembleRightHandSide(VectorType& vector) const = 0;
virtual ~TdcsRHSInterface(){};
};
}
#endif // DUNEURO_TDCS_RHS_INTERFACE_HH
set(duneuro_SOURCES ${CMAKE_SOURCE_DIR}/duneuro/meeg/fitted_meeg_driver_2d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/tes/fitted_tdcs_driver_2d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/tes/fitted_tdcs_point_driver_2d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/meeg/fitted_meeg_driver_3d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/tes/fitted_tdcs_driver_3d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/tes/fitted_tdcs_point_driver_3d.cc)
set(duneuro_SOURCES ${duneuro_SOURCES} ${CMAKE_SOURCE_DIR}/duneuro/common/deprecated.cc)
if (dune-udg_FOUND)
......