Commit 16846497 authored by Andreas Nüßing's avatar Andreas Nüßing

[CutFEM] add first implementation of cutfem

We add a first implementation of a cutfem MEEG driver. This commit basically
adds new necessary files and modifies the existing UDG files. They now
behave in a similar way as the FittedMEEGDriver for the distinction between
CG and DG.
parent 93a3472f
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNEURO_CUTFEM_GRIDOPERATOR_HH
#define DUNEURO_CUTFEM_GRIDOPERATOR_HH
#include <dune/istl/io.hh>
#include <dune/pdelab/gridfunctionspace/interpolate.hh>
#include <dune/pdelab/gridoperator/common/borderdofexchanger.hh>
#include <dune/pdelab/gridoperator/common/gridoperatorutilities.hh>
#include <dune/udg/pdelab/assembler/assembler.hh>
#include <dune/udg/pdelab/assembler/localassembler.hh>
#include <dune/udg/pdelab/ghostpenalty.hh>
namespace duneuro
{
template <class GO, class ST, class ENP>
class CutFEMGridOperator
{
public:
using Traits = typename GO::Traits;
//! Constructor for non trivial constraints
explicit CutFEMGridOperator(const GO& go, std::shared_ptr<const ST> subTriangulation,
const ENP& edgeNormProvider, const Dune::ParameterTree& config)
: go_(go)
, subTriangulation_(subTriangulation)
, edgeNormProvider_(edgeNormProvider)
, ghost_penalty_(config.get<double>("ghost_penalty"))
, conductivities_(config.get<std::vector<double>>("conductivities"))
{
}
void fill_pattern(typename GO::Pattern& p) const
{
go_.fill_pattern(p);
}
void residual(const typename GO::Domain& x, typename GO::Range& r) const
{
go_.residual(x, r);
}
//! Assembler jacobian
void jacobian(const typename GO::Domain& x, typename GO::Jacobian& a) const
{
go_.jacobian(x, a);
Dune::UDG::add_ghost_penalty_cutfem(*subTriangulation_, go_.trialGridFunctionSpace(),
edgeNormProvider_, ghost_penalty_, conductivities_, 0,
Dune::PDELab::Backend::native(a));
}
const typename GO::Traits::TrialGridFunctionSpace& trialGridFunctionSpace() const
{
return go_.trialGridFunctionSpace();
}
const typename GO::Traits::TestGridFunctionSpace& testGridFunctionSpace() const
{
return go_.testGridFunctionSpace();
}
const typename GO::Traits::MatrixBackend& matrixBackend() const
{
return go_.matrixBackend();
}
private:
const GO& go_;
std::shared_ptr<const ST> subTriangulation_;
ENP edgeNormProvider_;
double ghost_penalty_;
std::vector<double> conductivities_;
};
}
#endif
#ifndef DUNEURO_CUTFEM_MULTI_PHASE_SPACE_HH
#define DUNEURO_CUTFEM_MULTI_PHASE_SPACE_HH
#include <dune/localfunctions/lagrange.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/udg/pdelab/finiteelementmap.hh>
#include <dune/udg/vtriangulation.hh>
#include <duneuro/common/grid_function_space_utilities.hh>
namespace duneuro
{
struct CutFEMLeafOrderingParams : public Dune::PDELab::NoConstOrderingSize<true> {
};
template <class TGV, class N, int degree, int phases>
class CutFEMMultiPhaseSpace
{
public:
typedef TGV GV;
enum { dim = GV::dimension };
typedef typename GV::ctype ctype;
typedef N NT;
static const int blockSize = 1;
typedef Dune::LagrangeLocalFiniteElement<Dune::EquidistantPointSet, dim, double, double> LFE;
typedef VirtualSubTriangulation<GV> SubTriangulation;
typedef Dune::PDELab::istl::VectorBackend<> VBE;
typedef Dune::PDELab::UnfittedFiniteElementMapTraits<LFE, typename SubTriangulation::EntityPart>
UFEMTraits;
typedef Dune::PDELab::UnfittedFiniteElementMap<UFEMTraits, SubTriangulation> FEM;
typedef Dune::PDELab::LeafOrderingTag<CutFEMLeafOrderingParams> LeafOrderingTag;
typedef Dune::PDELab::GridFunctionSpace<GV, FEM, Dune::PDELab::NoConstraints, VBE,
LeafOrderingTag>
DomainGFS;
typedef Dune::PDELab::istl::VectorBackend<> PVBE;
typedef Dune::PDELab::PowerGridFunctionSpace<DomainGFS, phases, PVBE,
Dune::PDELab::EntityBlockedOrderingTag>
GFS;
typedef typename Dune::PDELab::Backend::Vector<GFS, N> DOF;
CutFEMMultiPhaseSpace(const GV& gv, std::shared_ptr<SubTriangulation> subTriangulation)
: gridView_(gv)
, entitySet_(gridView_)
, subTriangulation_(subTriangulation)
, lfe_(Dune::GeometryType(Dune::GeometryType::BasicType::cube, dim), degree)
{
for (unsigned int i = 0; i < phases; ++i) {
fems_[i] = std::make_shared<FEM>(lfe_, *subTriangulation_, i, false);
domainGfss_[i] = std::make_shared<DomainGFS>(entitySet_, *(fems_[i]));
}
gfs_ = make_power_gfs<DomainGFS, PVBE, typename GFS::Traits::OrderingTag>(domainGfss_);
gfs_->ordering();
}
// return gfs reference
GFS& getGFS()
{
return *gfs_;
}
// return gfs reference const version
const GFS& getGFS() const
{
return *gfs_;
}
CutFEMMultiPhaseSpace(const CutFEMMultiPhaseSpace&) = delete;
CutFEMMultiPhaseSpace& operator=(const CutFEMMultiPhaseSpace&) = delete;
private:
GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_;
std::shared_ptr<SubTriangulation> subTriangulation_;
LFE lfe_;
std::array<std::shared_ptr<FEM>, phases> fems_;
std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_;
std::shared_ptr<GFS> gfs_;
};
}
#endif // DUNEURO_CUTFEM_MULTI_PHASE_SPACE_HH
#ifndef DUNEURO_CUTFEM_SOLVER_HH
#define DUNEURO_CUTFEM_SOLVER_HH
#include <dune/udg/pdelab/cutfemmultiphaseoperator.hh>
#include <dune/udg/pdelab/multiphaseoperator.hh>
#include <dune/udg/pdelab/operator.hh>
#include <dune/udg/pdelab/subtriangulation.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/stationary/linearproblem.hh>
#include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/convection_diffusion_udg_default_parameter.hh>
#include <duneuro/common/cutfem_gridoperator.hh>
#include <duneuro/common/cutfem_multi_phase_space.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/random.hh>
namespace duneuro
{
template <class ST, int comps, int degree, class P, class DF, class RF, class JF>
struct CutFEMSolverTraits {
using SubTriangulation = ST;
using FundamentalGridView = typename ST::BaseT::GridView;
static const int dimension = FundamentalGridView::dimension;
static const int compartments = comps;
using Problem = P;
using FunctionSpace = CutFEMMultiPhaseSpace<FundamentalGridView, RF, degree, compartments>;
using DomainField = DF;
using RangeField = RF;
using DomainDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, DF>;
using RangeDOFVector = Dune::PDELab::Backend::Vector<typename FunctionSpace::GFS, RF>;
using EdgeNormProvider = MultiEdgeNormProvider;
using LocalOperator = ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider>;
using WrappedLocalOperator = Dune::UDG::CutFEMMultiPhaseLocalOperatorWrapper<LocalOperator>;
// using WrappedLocalOperator = Dune::UDG::MultiPhaseLocalOperatorWrapper<LocalOperator>;
using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<FundamentalGridView>;
using MatrixBackend = Dune::PDELab::istl::BCRSMatrixBackend<>;
using RawGridOperator =
Dune::UDG::UDGGridOperator<typename FunctionSpace::GFS, typename FunctionSpace::GFS,
WrappedLocalOperator, MatrixBackend, DF, RF, JF,
UnfittedSubTriangulation>;
using GridOperator = CutFEMGridOperator<RawGridOperator, ST, EdgeNormProvider>;
using SolverBackend = Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GridOperator>;
using LinearSolver = LinearProblemSolver<GridOperator, DomainDOFVector, RangeDOFVector>;
};
template <class ST, int compartments, int degree,
class P = ConvectionDiffusion_UDG_DefaultParameter<typename ST::BaseT::GridView>,
class DF = double, class RF = double, class JF = double>
class CutFEMSolver
{
public:
using Traits = CutFEMSolverTraits<ST, compartments, degree, P, DF, RF, JF>;
CutFEMSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation,
const Dune::ParameterTree& config)
: CutFEMSolver(subTriangulation, std::make_shared<typename Traits::Problem>(
config.get<std::vector<double>>("conductivities")),
config)
{
}
CutFEMSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<typename Traits::Problem> problem,
const Dune::ParameterTree& config)
: subTriangulation_(subTriangulation)
, problem_(problem)
, functionSpace_(subTriangulation_->gridView(), subTriangulation_)
, edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
, localOperator_(*problem_, edgeNormProvider_, ConvectionDiffusion_DG_Scheme::fromString(
config.get<std::string>("scheme")),
ConvectionDiffusion_DG_Weights::weightsOn, config.get<RF>("penalty"))
, wrappedLocalOperator_(localOperator_)
, unfittedSubTriangulation_(subTriangulation_->gridView(), *subTriangulation_)
, rawGridOperator_(functionSpace_.getGFS(), functionSpace_.getGFS(),
unfittedSubTriangulation_, wrappedLocalOperator_,
typename Traits::MatrixBackend(2 * Traits::dimension + 1))
, gridOperator_(rawGridOperator_, subTriangulation, edgeNormProvider_, config)
, linearSolver_(gridOperator_, config)
{
}
template <class SolverBackend>
void solve(SolverBackend& solverBackend, const typename Traits::RangeDOFVector& rightHandSide,
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));
linearSolver_.apply(solverBackend, solution, rightHandSide, config, dataTree);
dataTree.set("time", timer.elapsed());
}
template <class 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));
linearSolver_.apply(solverBackend, solution, config, dataTree);
dataTree.set("time", timer.elapsed());
}
const typename Traits::FunctionSpace& functionSpace() const
{
return functionSpace_;
}
const typename Traits::SubTriangulation& subTriangulation() const
{
return *subTriangulation_;
}
typename Traits::Problem& problem()
{
return *problem_;
}
private:
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_;
std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_;
typename Traits::LocalOperator localOperator_;
typename Traits::WrappedLocalOperator wrappedLocalOperator_;
typename Traits::UnfittedSubTriangulation unfittedSubTriangulation_;
typename Traits::RawGridOperator rawGridOperator_;
typename Traits::GridOperator gridOperator_;
typename Traits::LinearSolver linearSolver_;
};
}
#endif // DUNEURO_CUTFEM_SOLVER_HH
#ifndef DUNEURO_CUTFEM_SOLVER_BACKEND_HH
#define DUNEURO_CUTFEM_SOLVER_BACKEND_HH
#include <dune/pdelab/backend/istl.hh>
namespace duneuro
{
template <typename Solver>
struct CutFEMSolverBackendTraits {
using SolverBackend =
Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<typename Solver::Traits::GridOperator>;
};
template <typename Solver>
class CutFEMSolverBackend
{
public:
using Traits = CutFEMSolverBackendTraits<Solver>;
explicit CutFEMSolverBackend(std::shared_ptr<Solver> solver, const Dune::ParameterTree& config)
: solver_(solver)
, solverBackend_(config.get<unsigned int>("max_iterations", 5000),
config.get<unsigned int>("verbose", 0), true)
{
}
const typename Traits::SolverBackend& get() const
{
return solverBackend_;
}
typename Traits::SolverBackend& get()
{
return solverBackend_;
}
private:
std::shared_ptr<Solver> solver_;
typename Traits::SolverBackend solverBackend_;
};
}
#endif // DUNEURO_CUTFEM_SOLVER_BACKEND_HH
......@@ -31,6 +31,7 @@ namespace duneuro
};
enum class FittedSolverType { cg, dg };
enum class UnfittedSolverType { cutfem, udg};
}
#endif // DUNEURO_FLAGS_HH
......@@ -7,11 +7,13 @@ namespace duneuro
{
template <class GFS, int k, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, k, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, k>& gfss, Backend backend, OrderingTag tag);
make_power_gfs(std::array<std::shared_ptr<GFS>, k>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag());
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 1, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 1>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 1>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 1, Backend, OrderingTag>>(
*(gfss[0]), backend, tag);
......@@ -19,7 +21,8 @@ namespace duneuro
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 2, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 2>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 2>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 2, Backend, OrderingTag>>(
*(gfss[0]), *(gfss[1]), backend, tag);
......@@ -27,7 +30,8 @@ namespace duneuro
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 3, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 3>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 3>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 3, Backend, OrderingTag>>(
*(gfss[0]), *(gfss[1]), *(gfss[2]), backend, tag);
......@@ -35,7 +39,8 @@ namespace duneuro
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 4, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 4>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 4>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 4, Backend, OrderingTag>>(
*(gfss[0]), *(gfss[1]), *(gfss[2]), *(gfss[3]), backend, tag);
......@@ -43,7 +48,8 @@ namespace duneuro
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 5, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 5>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 5>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 5, Backend, OrderingTag>>(
*(gfss[0]), *(gfss[1]), *(gfss[2]), *(gfss[3]), *(gfss[4]), backend, tag);
......@@ -51,7 +57,8 @@ namespace duneuro
template <class GFS, class Backend, class OrderingTag>
std::shared_ptr<Dune::PDELab::PowerGridFunctionSpace<GFS, 6, Backend, OrderingTag>>
make_power_gfs(std::array<std::shared_ptr<GFS>, 6>& gfss, Backend backend, OrderingTag tag)
make_power_gfs(std::array<std::shared_ptr<GFS>, 6>& gfss, Backend backend = Backend(),
OrderingTag tag = OrderingTag())
{
return std::make_shared<Dune::PDELab::PowerGridFunctionSpace<GFS, 6, Backend, OrderingTag>>(
*(gfss[0]), *(gfss[1]), *(gfss[2]), *(gfss[3]), *(gfss[4]), *(gfss[5]), backend, tag);
......
......@@ -49,7 +49,7 @@ namespace duneuro
fems_[i] = std::make_shared<FEM>(lfe_, *subTriangulation_, i);
domainGfss_[i] = std::make_shared<DomainGFS>(entitySet_, *(fems_[i]));
}
gfs_ = make_power_gfs(domainGfss_, PVBE(), OrderingTag(blockSize));
gfs_ = duneuro::make_power_gfs(domainGfss_, PVBE(), OrderingTag(blockSize));
gfs_->ordering();
}
......@@ -67,6 +67,7 @@ namespace duneuro
UDGQkMultiPhaseSpace(const UDGQkMultiPhaseSpace&) = delete;
UDGQkMultiPhaseSpace& operator=(const UDGQkMultiPhaseSpace&) = delete;
private:
GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_;
......
......@@ -17,11 +17,12 @@
namespace duneuro
{
template <class ST, int compartments, int degree, class P, class DF, class RF, class JF>
template <class ST, int comps, int degree, class P, class DF, class RF, class JF>
struct UDGSolverTraits {
using SubTriangulation = ST;
using FundamentalGridView = typename ST::BaseT::GridView;
static const int dimension = FundamentalGridView::dimension;
static const int compartments = comps;
using Problem = P;
using FunctionSpace = UDGQkMultiPhaseSpace<FundamentalGridView, RF, degree, compartments>;
using DomainField = DF;
......
#ifndef DUNEURO_CUTFEM_SOURCE_MODEL_FACTORY_HH
#define DUNEURO_CUTFEM_SOURCE_MODEL_FACTORY_HH
#include <dune/common/parametertree.hh>
#include <dune/common/std/memory.hh>
#include <duneuro/common/exceptions.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/udg_partial_integration_source_model.hh>
namespace duneuro
{
struct CutFEMSourceModelFactory {
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,
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, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, false);
} else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\"");
}
}
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,
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, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, false);
} else {
DUNE_THROW(duneuro::SourceModelException, "unknown source model of type \"" << type
<< "\"");
}
}
};
}
#endif // DUNEURO_CUTFEM_SOURCE_MODEL_FACTORY_HH
......@@ -13,29 +13,25 @@
namespace duneuro
{
template <class ST, int compartments, int degree, class DF, class RF, class JF>
template <class S, class SMF>
struct UDGEEGFowardSolverTraits {
static const unsigned int dimension = ST::dim;
using Solver =
UDGSolver<ST, compartments, degree,
ConvectionDiffusion_UDG_DefaultParameter<typename ST::GridView>, DF, RF, JF>;
using SubTriangulation = ST;
using Solver = S;
using SubTriangulation = typename Solver::Traits::SubTriangulation;
static const unsigned int dimension = SubTriangulation::dim;
using FunctionSpace = typename Solver::Traits::FunctionSpace;
using DomainDOFVector = typename Solver::Traits::DomainDOFVector;
using RangeDOFVector = typename Solver::Traits::RangeDOFVector;
using CoordinateFieldType = typename ST::ctype;
using CoordinateFieldType = typename SubTriangulation::ctype;
using DipoleType = Dipole<CoordinateFieldType, dimension>;
using ElementSearch = KDTreeElementSearch<typename ST::BaseT::GridView>;
using ElementSearch = KDTreeElementSearch<typename SubTriangulation::BaseT::GridView>;
};
template <class ST, int compartments, int degree, class DF = double, class RF = double,
class JF = double>
template <class S, class SMF>
class UDGEEGFowardSolver
: public EEGForwardSolver<UDGEEGFowardSolver<ST, compartments, degree, DF, RF, JF>,
UDGEEGFowardSolverTraits<ST, compartments, degree, DF, RF, JF>>
: public EEGForwardSolver<UDGEEGFowardSolver<S, SMF>, UDGEEGFowardSolverTraits<S, SMF>>
{
public:
using Traits = UDGEEGFowardSolverTraits<ST, compartments, degree, DF, RF, JF>;
using Traits = UDGEEGFowardSolverTraits<S, SMF>;
UDGEEGFowardSolver(std::shared_ptr<typename Traits::SubTriangulation> subTriangulation,
std::shared_ptr<typename Traits::Solver> solver,
......@@ -51,9 +47,8 @@ namespace duneuro
void setSourceModel(const Dune::ParameterTree& config, DataTree dataTree = DataTree())
{
denseSourceModel_ =
UDGSourceModelFactory::template createDense<typename Traits::RangeDOFVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
denseSourceModel_ = SMF::template createDense<typename Traits::RangeDOFVector>(
*solver_, subTriangulation_, search_, config.get<std::size_t>("compartment"), config);
}
void bind(const typename Traits::DipoleType& dipole)
......@@ -101,7 +96,7 @@ namespace duneuro
private:
std::shared_ptr<typename Traits::SubTriangulation> subTriangulation_;
std::shared_ptr<typename Traits::Solver> solver_;
std::shared_ptr<KDTreeElementSearch<typename ST::BaseT::GridView>> search_;
std::shared_ptr<typename Traits::ElementSearch> search_;
std::unique_ptr<typename Traits::RangeDOFVector> rightHandSideVector_;
std::shared_ptr<SourceModelInterface<typename Traits::CoordinateFieldType, Traits::dimension,
typename Traits::RangeDOFVector>>
......
......@@ -35,10 +35,11 @@ namespace duneuro
UDGPartialIntegrationSourceModel(const GFS& gfs, std::shared_ptr<ST> subTriangulation,
std::shared_ptr<typename BaseT::SearchType> search,
std::size_t child)
std::size_t child, bool scaleToBBox)
: BaseT(search)
, subTriangulation_(subTriangulation)
, child_(child)
, scaleToBBox_(scaleToBBox)
, ulfs_(gfs)
, ucache_(ulfs_)
{
......@@ -76,8 +77,11 @@ namespace duneuro
// reset basis transformations
FESwitch::basis(childLfs.finiteElement()).reset();
auto boundingBoxLocal = ep.subEntity().boundingBox().local(
this->dipoleElement().geometry().global(this->localDipolePosition()));
auto boundingBoxLocal =
scaleToBBox_ ?
ep.subEntity().boundingBox().local(
this->dipoleElement().geometry().global(this->localDipolePosition())) :
this->localDipolePosition();
// evaluate gradiant in local bounding box coordinates
std::vector<Dune::FieldMatrix<Real, 1, dim>> gradpsi(childLfs.size());
......@@ -86,11 +90,16 @@ namespace duneuro
// get transformation of the boundingbox
const auto& boundingBoxJacobian =
ep.subEntity().boundingBox().jacobianInverseTransposed(CoordinateType(0.5));
const auto& entityJacobian =
ep.entity().geometry().jacobianInverseTransposed(CoordinateType(.5));
// multiply local gradiant with bounding box transformation
Dune::FieldMatrix<Real, 1, dim> tmp;
for (unsigned int i = 0; i < gradpsi.size(); i++) {
boundingBoxJacobian.mv(gradpsi[i][0], tmp[0]);
if (scaleToBBox_)
boundingBoxJacobian.mv(gradpsi[i][0], tmp[0]);
else
entityJacobian.mv(gradpsi[i][0], tmp[0]);
gradpsi[i] = tmp;
}
......@@ -111,6 +120,7 @@ namespace duneuro
private:
std::shared_ptr<ST> subTriangulation_;
std::size_t child_;
bool scaleToBBox_;
mutable ULFS ulfs_;
mutable UCache ucache_;
};
......
......@@ -24,7 +24,7 @@ namespace duneuro
if (type == "partial_integration") {
return Dune::Std::make_unique<UDGPartialIntegrationSourceModel<
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment);
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, true);
} else if (type == "patch_based_venant") {
return Dune::Std::make_unique<UDGPatchBasedVenantSourceModel<
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
......@@ -47,7 +47,7 @@ namespace duneuro
if (type == "partial_integration") {
return Dune::Std::make_unique<UDGPartialIntegrationSourceModel<
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment);
solver.functionSpace().getGFS(), subTriangulation, search, dipoleCompartment, true);
} else if (type == "patch_based_venant") {
return Dune::Std::make_unique<UDGPatchBasedVenantSourceModel<
typename Solver::Traits::FunctionSpace::GFS, ST, Vector>>(
......
......@@ -27,8 +27,9 @@ namespace duneuro
using ULFS = Dune::PDELab::UnfittedLocalFunctionSpace<GFS>;
using UCache = Dune::PDELab::LFSIndexCache<ULFS>;
UDGTransferMatrixRHS(const GFS& gfs, std::shared_ptr<ST> st, std::size_t child)
: st_(st), gfs_(gfs), child_(child)
UDGTransferMatrixRHS(const GFS& gfs, std::shared_ptr<ST> st, std::size_t child,
bool scaleToBBox)
: st_(st), gfs_(gfs), child_(child), scaleToBBox_(scaleToBBox)
{
}
......@@ -66,9 +67,11 @@ namespace duneuro
}
std::vector<RangeType> phi(childLfs.size());
FESwitch::basis(childLfs.finiteElement())
.evaluateFunction(
ep.boundingBox().local(electrodeElement.geometry().global(electrodeLocal)), phi);
auto local =
scaleToBBox_ ?
ep.boundingBox().local(electrodeElement.geometry().global(electrodeLocal)) :
electrodeLocal;
FESwitch::basis(childLfs.finiteElement()).evaluateFunction(local, phi);
for (unsigned int i = 0; i < childLfs.size(); ++i) {
output[ucache.containerIndex(childLfs.localIndex(i))] = phi[i];
......@@ -88,9 +91,11 @@ namespace duneuro
assert(childLfs.size() > 0);
std::vector<RangeType> phi(childLfs.size());
FESwitch::basis(childLfs.finiteElement())
.evaluateFunction(
ep.boundingBox().local(referenceElement.geometry().global(referenceLocal)), phi);
auto local =
scaleToBBox_ ?
ep.boundingBox().local(referenceElement.geometry().global(referenceLocal)) :
referenceLocal;
FESwitch::basis(childLfs.finiteElement()).evaluateFunction(local, phi);
for (unsigned int i = 0; i < childLfs.size(); ++i) {
output[ucache.containerIndex(childLfs.localIndex(i))] -= phi[i];
......@@ -103,6 +108,7 @@ namespace duneuro
std::shared_ptr<ST> st_;
const GFS& gfs_;
std::size_t child_;
bool scaleToBBox_;
};
}
......
......@@ -11,34 +11,32 @@
namespace duneuro
{
template <class ST, int compartments, int degree, class DF, class RF, class JF>