Commit 1758a0c6 authored by Andreas Nüßing's avatar Andreas Nüßing

[PenaltyFluxWeighting] extract weight computation from LOP

For the DG, CutFEM and UDG methods, there were several places where the
weights for the penalty and average terms were computed. Besides being
just copy and pasted and thus error prone, it is also not very extensible
when one wants to try different weighting strategies.
This commit extracts the penalty and average weights from the local
operators and transfers them into an additional class, similar to the
edge norm provider.
parent fd23afb1
......@@ -15,6 +15,7 @@
#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/penalty_flux_weighting.hh>
#include <duneuro/common/random.hh>
#include <duneuro/common/vector_initialization.hh>
......@@ -33,7 +34,9 @@ namespace duneuro
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 PenaltyFluxWeighting = UnfittedDynamicPenaltyFluxWeights;
using LocalOperator =
ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
using WrappedLocalOperator = Dune::UDG::CutFEMMultiPhaseLocalOperatorWrapper<LocalOperator>;
// using WrappedLocalOperator = Dune::UDG::MultiPhaseLocalOperatorWrapper<LocalOperator>;
using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<FundamentalGridView>;
......@@ -70,9 +73,11 @@ namespace duneuro
, 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"))
, weighting_(config.get<std::string>("weighting"))
, localOperator_(
*problem_, edgeNormProvider_, weighting_,
ConvectionDiffusion_DG_Scheme::fromString(config.get<std::string>("scheme")),
config.get<RF>("penalty"))
, wrappedLocalOperator_(localOperator_)
, unfittedSubTriangulation_(subTriangulation_->gridView(), *subTriangulation_)
, rawGridOperator_(functionSpace_.getGFS(), functionSpace_.getGFS(),
......@@ -133,6 +138,7 @@ namespace duneuro
std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_;
typename Traits::PenaltyFluxWeighting weighting_;
typename Traits::LocalOperator localOperator_;
typename Traits::WrappedLocalOperator wrappedLocalOperator_;
typename Traits::UnfittedSubTriangulation unfittedSubTriangulation_;
......
......@@ -11,6 +11,7 @@
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/flags.hh>
#include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/random.hh>
#include <duneuro/io/data_tree.hh>
......@@ -44,7 +45,9 @@ namespace duneuro
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 PenaltyFluxWeighting = FittedDynamicPenaltyFluxWeights;
using LocalOperator =
ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
using Assembler = GalerkinGlobalAssembler<FunctionSpace, LocalOperator, DF, RF, JF>;
using LinearSolver =
LinearProblemSolver<typename Assembler::GO, DomainDOFVector, RangeDOFVector>;
......@@ -64,11 +67,10 @@ namespace duneuro
, problem_(problem)
, functionSpace_(volumeConductor_->gridView())
, edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
, weighting_(config.get<std::string>("weighting"))
, localOperator_(
*problem_, edgeNormProvider_,
*problem_, edgeNormProvider_, weighting_,
ConvectionDiffusion_DG_Scheme::fromString(config.get<std::string>("scheme")),
config.get<bool>("weights") ? ConvectionDiffusion_DG_Weights::weightsOn :
ConvectionDiffusion_DG_Weights::weightsOff,
config.get<DF>("penalty"), false, config.get<unsigned int>("intorderadd", 0))
, assembler_(functionSpace_, localOperator_, 2 * VC::dim + 1)
, linearSolver_(*assembler_, config)
......@@ -146,6 +148,7 @@ namespace duneuro
std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_;
typename Traits::PenaltyFluxWeighting weighting_;
typename Traits::LocalOperator localOperator_;
typename Traits::Assembler assembler_;
typename Traits::LinearSolver linearSolver_;
......
#ifndef DUNEURO_PENALTY_FLUX_WEIGHTING_HH
#define DUNEURO_PENALTY_FLUX_WEIGHTING_HH
#include <dune/common/fmatrix.hh>
#include <dune/geometry/referenceelements.hh>
namespace duneuro
{
template <class T>
struct PenaltyFluxWeights {
T penaltyWeight;
T fluxInsideWeight;
T fluxOutsideWeight;
};
/**
* harmonic average for penalty weight and weighting of the flux only based on the tensors of
* inside and outside
*/
class TensorOnlyPenaltyFluxWeights
{
public:
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& intersection, const T& insideTensor,
const T& outsideTensor) const
{
using DF = typename IG::ctype;
static const int dim = IG::dimension;
using GlobalCoordinate = Dune::FieldVector<DF, dim>;
using LocalCoordinate = Dune::FieldVector<DF, dim - 1>;
// compute normal
// workaround for intersections that do not provide the centerUnitOutderNormal function
const LocalCoordinate localCenter =
Dune::ReferenceElements<DF, dim - 1>::general(intersection.geometry().type())
.position(0, 0);
const GlobalCoordinate normal = intersection.unitOuterNormal(localCenter);
// compute tensor in normal direction
GlobalCoordinate insideTensorNormal, outsideTensorNormal;
insideTensor.mv(normal, insideTensorNormal);
outsideTensor.mv(normal, outsideTensorNormal);
const auto deltaInside = insideTensorNormal * normal;
const auto deltaOutside = outsideTensorNormal * normal;
// compute result
PenaltyFluxWeights<DF> result;
result.penaltyWeight =
2.0 * deltaInside * deltaOutside / (deltaInside + deltaOutside + 1e-20);
result.fluxInsideWeight = deltaOutside / (deltaInside + deltaOutside + 1e-20);
result.fluxOutsideWeight = deltaInside / (deltaInside + deltaOutside + 1e-20);
return result;
}
};
/**
* no weighting for the penalty and a .5 average of the flux
*/
class ConstantPenaltyFluxWeights
{
public:
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& DUNE_UNUSED(intersection),
const T& DUNE_UNUSED(insideTensor),
const T& DUNE_UNUSED(outsideTensor)) const
{
using DF = typename IG::ctype;
PenaltyFluxWeights<DF> result;
result.penaltyWeight = 1.0;
result.fluxInsideWeight = .5;
result.fluxOutsideWeight = .5;
return result;
}
};
class AnnavarapuPenaltyFluxWeights
{
public:
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& intersection, const T& insideTensor,
const T& outsideTensor) const
{
using DF = typename IG::ctype;
static const int dim = IG::dimension;
using GlobalCoordinate = Dune::FieldVector<DF, dim>;
using LocalCoordinate = Dune::FieldVector<DF, dim - 1>;
// compute normal
// workaround for intersections that do not provide the centerUnitOutderNormal function
const LocalCoordinate localCenter =
Dune::ReferenceElements<DF, dim - 1>::general(intersection.geometry().type())
.position(0, 0);
const GlobalCoordinate normal = intersection.unitOuterNormal(localCenter);
// compute tensor in normal direction
GlobalCoordinate insideTensorNormal, outsideTensorNormal;
insideTensor.mv(normal, insideTensorNormal);
outsideTensor.mv(normal, outsideTensorNormal);
const auto deltaInside = insideTensorNormal * normal;
const auto deltaOutside = outsideTensorNormal * normal;
// compute volumes of the different cut parts
const auto iv = intersection.intersection().insideVolume();
const auto ov = intersection.intersection().outsideDomainIndex() == -1 ?
iv :
intersection.intersection().outsideVolume();
const auto fv = intersection.intersection().area();
// compute result
PenaltyFluxWeights<DF> result;
result.penaltyWeight = fv / (deltaInside / iv + deltaOutside / ov + 1e-20);
result.fluxInsideWeight = deltaOutside * iv / (deltaOutside * iv + deltaInside * ov + 1e-20);
result.fluxOutsideWeight = deltaInside * ov / (deltaOutside * iv + deltaInside * ov + 1e-20);
return result;
}
};
class BarrauPenaltyFluxWeights
{
public:
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& intersection, const T& insideTensor,
const T& outsideTensor) const
{
using DF = typename IG::ctype;
static const int dim = IG::dimension;
using GlobalCoordinate = Dune::FieldVector<DF, dim>;
using LocalCoordinate = Dune::FieldVector<DF, dim - 1>;
// compute normal
// workaround for intersections that do not provide the centerUnitOutderNormal function
const LocalCoordinate localCenter =
Dune::ReferenceElements<DF, dim - 1>::general(intersection.geometry().type())
.position(0, 0);
const GlobalCoordinate normal = intersection.unitOuterNormal(localCenter);
// compute tensor in normal direction
GlobalCoordinate insideTensorNormal, outsideTensorNormal;
insideTensor.mv(normal, insideTensorNormal);
outsideTensor.mv(normal, outsideTensorNormal);
const auto deltaInside = insideTensorNormal * normal;
const auto deltaOutside = outsideTensorNormal * normal;
// compute volumes of the different cut parts
const auto iv = intersection.intersection().insideVolume();
const auto ov = intersection.intersection().outsideDomainIndex() == -1 ?
iv :
intersection.intersection().outsideVolume();
const auto fv = intersection.intersection().area();
// compute result
PenaltyFluxWeights<DF> result;
result.penaltyWeight =
std::max(deltaInside * iv / (iv + ov), deltaOutside * ov / (iv + ov)) * fv / (iv + ov);
result.fluxInsideWeight = deltaOutside * iv / (deltaOutside * iv + deltaInside * ov + 1e-20);
result.fluxOutsideWeight = deltaInside * ov / (deltaOutside * iv + deltaInside * ov + 1e-20);
return result;
}
};
enum class PenaltyFluxWeightsTypes { constant, tensorOnly, annavarapu, barrau };
static inline PenaltyFluxWeightsTypes penaltyFluxWeightingFromString(const std::string& name)
{
if (name == "constant")
return PenaltyFluxWeightsTypes::constant;
else if (name == "tensorOnly")
return PenaltyFluxWeightsTypes::tensorOnly;
else if (name == "annavarapu")
return PenaltyFluxWeightsTypes::annavarapu;
else if (name == "barrau")
return PenaltyFluxWeightsTypes::barrau;
else
DUNE_THROW(Dune::Exception, "unknown weighting type \"" << name << "\"");
}
class FittedDynamicPenaltyFluxWeights
{
public:
explicit FittedDynamicPenaltyFluxWeights(PenaltyFluxWeightsTypes type) : type_(type)
{
}
explicit FittedDynamicPenaltyFluxWeights(const std::string type)
: FittedDynamicPenaltyFluxWeights(penaltyFluxWeightingFromString(type))
{
}
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& intersection, const T& insideTensor,
const T& outsideTensor) const
{
switch (type_) {
case PenaltyFluxWeightsTypes::constant:
return ConstantPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
case PenaltyFluxWeightsTypes::tensorOnly:
return TensorOnlyPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
default: DUNE_THROW(Dune::Exception, "illegal weighting for a DG method");
}
}
private:
PenaltyFluxWeightsTypes type_;
};
class UnfittedDynamicPenaltyFluxWeights
{
public:
explicit UnfittedDynamicPenaltyFluxWeights(PenaltyFluxWeightsTypes type) : type_(type)
{
}
explicit UnfittedDynamicPenaltyFluxWeights(const std::string& type)
: UnfittedDynamicPenaltyFluxWeights(penaltyFluxWeightingFromString(type))
{
}
template <class IG, class T>
PenaltyFluxWeights<typename IG::ctype> operator()(const IG& intersection, const T& insideTensor,
const T& outsideTensor) const
{
switch (type_) {
case PenaltyFluxWeightsTypes::constant:
return ConstantPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
case PenaltyFluxWeightsTypes::tensorOnly:
return TensorOnlyPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
case PenaltyFluxWeightsTypes::annavarapu:
return AnnavarapuPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
case PenaltyFluxWeightsTypes::barrau:
return BarrauPenaltyFluxWeights()(intersection, insideTensor, outsideTensor);
default: DUNE_THROW(Dune::Exception, "illegal weighting for a DG method");
}
}
private:
PenaltyFluxWeightsTypes type_;
};
}
#endif // DUNEURO_PENALTY_FLUX_WEIGHTING_HH
......@@ -16,7 +16,7 @@
#include <duneuro/common/convection_diffusion_udg_default_parameter.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/linear_problem_solver.hh>
#include <duneuro/common/random.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/udg_multi_phase_space.hh>
#include <duneuro/common/vector_initialization.hh>
......@@ -35,7 +35,9 @@ namespace duneuro
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 PenaltyFluxWeighting = UnfittedDynamicPenaltyFluxWeights;
using LocalOperator =
ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
using WrappedLocalOperator = Dune::UDG::MultiPhaseLocalOperatorWrapper<LocalOperator>;
using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<FundamentalGridView>;
using MatrixBackend = Dune::PDELab::istl::BCRSMatrixBackend<>;
......@@ -68,9 +70,11 @@ namespace duneuro
, 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"))
, weighting_(config.get<std::string>("weighting"))
, localOperator_(
*problem_, edgeNormProvider_, weighting_,
ConvectionDiffusion_DG_Scheme::fromString(config.get<std::string>("scheme")),
config.get<RF>("penalty"), false, config.get<RF>("intorderadd"))
, wrappedLocalOperator_(localOperator_)
, unfittedSubTriangulation_(subTriangulation_->gridView(), *subTriangulation_)
, gridOperator_(functionSpace_.getGFS(), functionSpace_.getGFS(), unfittedSubTriangulation_,
......@@ -130,6 +134,7 @@ namespace duneuro
std::shared_ptr<typename Traits::Problem> problem_;
typename Traits::FunctionSpace functionSpace_;
typename Traits::EdgeNormProvider edgeNormProvider_;
typename Traits::PenaltyFluxWeighting weighting_;
typename Traits::LocalOperator localOperator_;
typename Traits::WrappedLocalOperator wrappedLocalOperator_;
typename Traits::UnfittedSubTriangulation unfittedSubTriangulation_;
......
......@@ -13,6 +13,7 @@
#include <duneuro/common/element_patch.hh>
#include <duneuro/common/entityset_volume_conductor.hh>
#include <duneuro/common/logged_timer.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/common/sub_function_space.hh>
#include <duneuro/common/subset_entityset.hh>
#include <duneuro/eeg/source_model_interface.hh>
......@@ -39,7 +40,8 @@ namespace duneuro
using Problem =
SubtractionDGDefaultParameter<SubEntitySet, typename V::field_type, SubVolumeConductor>;
using EdgeNormProvider = MultiEdgeNormProvider;
using LOP = SubtractionDG<Problem, EdgeNormProvider>;
using PenaltyFluxWeighting = FittedDynamicPenaltyFluxWeights;
using LOP = SubtractionDG<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
using DOF = typename SUBFS::DOF;
using AS = Dune::PDELab::GalerkinGlobalAssembler<SUBFS, LOP, Dune::SolverCategory::sequential>;
using SubLFS = Dune::PDELab::LocalFunctionSpace<typename SUBFS::GFS>;
......@@ -59,12 +61,11 @@ namespace duneuro
, elementNeighborhoodMap_(std::make_shared<ElementNeighborhoodMap<typename VC::GridView>>(
volumeConductor_->gridView()))
, edgeNormProvider_(solverConfig.get<std::string>("edge_norm_type"), 1.0)
, weighting_(solverConfig.get<std::string>("weighting"))
, config_(config)
, intorderadd_lb_(config.get<unsigned int>("intorderadd_lb"))
, scheme_(
ConvectionDiffusion_DG_Scheme::fromString(solverConfig.get<std::string>("scheme")))
, weights_(solverConfig.get<bool>("weights") ? ConvectionDiffusion_DG_Weights::weightsOn :
ConvectionDiffusion_DG_Weights::weightsOff)
, penalty_(solverConfig.get<double>("penalty"))
{
}
......@@ -107,7 +108,7 @@ namespace duneuro
this->dipole().moment());
problem_ = std::make_shared<Problem>(subVolumeConductor_->entitySet(), subVolumeConductor_);
problem_->bind(this->dipoleElement(), this->localDipolePosition(), this->dipole().moment());
lop_ = std::make_shared<LOP>(*problem_, weights_, config_.get<unsigned int>("intorderadd"),
lop_ = std::make_shared<LOP>(*problem_, weighting_, config_.get<unsigned int>("intorderadd"),
config_.get<unsigned int>("intorderadd_lb"));
subFS_ = std::make_shared<SUBFS>(subVolumeConductor_);
dataTree.set("sub_dofs", subFS_->getGFS().size());
......@@ -162,6 +163,7 @@ namespace duneuro
std::shared_ptr<Problem> problem_;
std::shared_ptr<HostProblem> hostProblem_;
EdgeNormProvider edgeNormProvider_;
PenaltyFluxWeighting weighting_;
std::shared_ptr<LOP> lop_;
std::shared_ptr<SUBFS> subFS_;
std::shared_ptr<AS> assembler_;
......@@ -171,7 +173,6 @@ namespace duneuro
std::vector<typename HostGridView::Intersection> patchBoundaryIntersections_;
unsigned int intorderadd_lb_;
ConvectionDiffusion_DG_Scheme::Type scheme_;
ConvectionDiffusion_DG_Weights::Type weights_;
double penalty_;
void assembleLocalDefaultSubtraction(VectorType& vector) const
......@@ -227,23 +228,8 @@ namespace duneuro
const auto& A_n = volumeConductor_->tensor(outside);
const auto& n_F = is.centerUnitOuterNormal();
RF omega_s;
RF omega_n;
RF harmonic_average;
if (weights_ == ConvectionDiffusion_DG_Weights::weightsOn) {
Dune::FieldVector<RF, dim> An_F_s;
A_s.mv(n_F, An_F_s);
Dune::FieldVector<RF, dim> An_F_n;
A_n.mv(n_F, An_F_n);
RF delta_s = (An_F_s * n_F);
RF delta_n = (An_F_n * n_F);
omega_s = delta_n / (delta_s + delta_n + 1e-20);
omega_n = delta_s / (delta_s + delta_n + 1e-20);
harmonic_average = 2.0 * delta_s * delta_n / (delta_s + delta_n + 1e-20);
} else {
omega_s = omega_n = 0.5;
harmonic_average = 1.0;
}
auto weights = weighting_(is, A_s, A_n);
// note: edgenorm provider needs the intersectiongeometry interface
RF h_F;
......@@ -255,7 +241,7 @@ namespace duneuro
const int degree = std::max(order_s, order_n);
const RF penalty_factor =
(penalty_ / h_F) * harmonic_average * degree * (degree + VC::dim - 1);
(penalty_ / h_F) * weights.penaltyWeight * degree * (degree + VC::dim - 1);
const int intorder = intorderadd_lb_ + 2 * degree;
......@@ -298,14 +284,16 @@ namespace duneuro
for (unsigned int i = 0; i < hostcache_inside.size(); i++) {
auto index = hostcache_inside.containerIndex(i);
vector[index] += phi_s[i] * term1;
vector[index] += -phi_s[i] * omega_n * term1;
vector[index] += (Agradpsi_s[i] * n_F) * term2 * omega_s; // symmetry term
vector[index] += -phi_s[i] * weights.fluxOutsideWeight * term1;
vector[index] +=
(Agradpsi_s[i] * n_F) * term2 * weights.fluxInsideWeight; // symmetry term
vector[index] += -phi_s[i] * term3; // penalty term
}
for (unsigned int i = 0; i < hostcache_outside.size(); i++) {
auto index = hostcache_outside.containerIndex(i);
vector[index] += -phi_n[i] * omega_s * term1;
vector[index] += (Agradpsi_n[i] * n_F) * term2 * omega_n; // symmetry term
vector[index] += -phi_n[i] * weights.fluxInsideWeight * term1;
vector[index] +=
(Agradpsi_n[i] * n_F) * term2 * weights.fluxOutsideWeight; // symmetry term
vector[index] += phi_n[i] * term3; // penalty term
}
}
......
......@@ -25,19 +25,17 @@
namespace duneuro
{
/**** class definition ****/
template <class PROBLEMDATA>
template <class PROBLEMDATA, class PenaltyFluxWeighting>
class SubtractionDGLambda
{
public:
/*** constructor ***/
SubtractionDGLambda(
PROBLEMDATA& problem_, unsigned int intorderadd_ = 0, unsigned int intorderadd_lb_ = 0,
ConvectionDiffusion_DG_Weights::Type weights_ = ConvectionDiffusion_DG_Weights::weightsOn)
: param(problem_)
SubtractionDGLambda(PROBLEMDATA& problem_, const PenaltyFluxWeighting& weighting_,
unsigned int intorderadd_ = 0, unsigned int intorderadd_lb_ = 0)
: param(problem_), weighting(weighting_)
{
intorderadd = intorderadd_;
intorderadd_lb = intorderadd_lb_;
weights = weights_;
}
/*** lambda_boundary method ***/
......@@ -203,32 +201,13 @@ namespace duneuro
A_s = param.A(insideEntity, inside_local);
A_n = param.A(outsideEntity, outside_local);
/* evaluate normal at intersectino center */
const Dune::FieldVector<DF, dim> n_F = ig.centerUnitOuterNormal();
/* tensor times normal */
Dune::FieldVector<RF, dim> An_F_s;
A_s.mv(n_F, An_F_s);
Dune::FieldVector<RF, dim> An_F_n;
A_n.mv(n_F, An_F_n);
RF omega_s;
RF omega_n;
RF harmonic_average(0.0);
if (weights == ConvectionDiffusion_DG_Weights::weightsOn) {
RF delta_s = (An_F_s * n_F);
RF delta_n = (An_F_n * n_F);
omega_s = delta_n / (delta_s + delta_n + 1e-20);
omega_n = delta_s / (delta_s + delta_n + 1e-20);
harmonic_average = 2.0 * delta_s * delta_n / (delta_s + delta_n + 1e-20);
} else {
omega_s = omega_n = 0.5;
harmonic_average = 1.0;
}
auto weights = weighting(ig, A_s, A_n);
std::vector<RangeType> phi_s(lfsv_s.size());
std::vector<RangeType> phi_n(lfsv_n.size());
/** loop over quadrature points **/
for (const auto& qp: rule) {
for (const auto& qp : rule) {
/* unit outer normal */
const Dune::FieldVector<DF, dim> n_F_local = ig.unitOuterNormal(qp.position());
......@@ -264,8 +243,8 @@ namespace duneuro
}
/* combine the terms */
RF termls = (omega_s * (sigma_corr_grad_u_infty_s * n_F_local)
+ omega_n * (sigma_corr_grad_u_infty_n * n_F_local))
RF termls = (weights.fluxInsideWeight * (sigma_corr_grad_u_infty_s * n_F_local)
+ weights.fluxOutsideWeight * (sigma_corr_grad_u_infty_n * n_F_local))
* factor;
if (std::isnan(termls)) {
......@@ -282,9 +261,9 @@ namespace duneuro
private:
PROBLEMDATA& param;
PenaltyFluxWeighting weighting;
unsigned int intorderadd;
unsigned int intorderadd_lb;
ConvectionDiffusion_DG_Weights::Type weights;
};
}
......
......@@ -9,9 +9,9 @@ namespace duneuro
enum class SubtractionContinuityType { continuous, discontinuous };
/**** class definition of the operator ****/
template <typename PROBLEMDATA, typename EdgeNormProvider,
template <typename PROBLEMDATA, typename EdgeNormProvider, typename PenaltyFluxWeighting,
SubtractionContinuityType continuityType = SubtractionContinuityType::discontinuous>
class SubtractionDG : public SubtractionDGLambda<PROBLEMDATA>,
class SubtractionDG : public SubtractionDGLambda<PROBLEMDATA, PenaltyFluxWeighting>,
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::FullSkeletonPattern,
public Dune::PDELab::LocalOperatorDefaultFlags
......@@ -22,15 +22,15 @@ namespace duneuro
enum { doLambdaVolume = true };
enum { doLambdaSkeleton = continuityType == SubtractionContinuityType::discontinuous };
using SubtractionDGLambda<PROBLEMDATA>::lambda_volume;
using SubtractionDGLambda<PROBLEMDATA>::lambda_boundary;
using SubtractionDGLambda<PROBLEMDATA>::lambda_skeleton;
using SubtractionDGLambda<PROBLEMDATA, PenaltyFluxWeighting>::lambda_volume;
using SubtractionDGLambda<PROBLEMDATA, PenaltyFluxWeighting>::lambda_boundary;
using SubtractionDGLambda<PROBLEMDATA, PenaltyFluxWeighting>::lambda_skeleton;
/*** Constructor ***/
SubtractionDG(PROBLEMDATA& problem_, const ConvectionDiffusion_DG_Weights::Type weights_ =
ConvectionDiffusion_DG_Weights::weightsOn,
SubtractionDG(PROBLEMDATA& problem_, const PenaltyFluxWeighting& weighting_,
unsigned int intorderadd_ = 0, unsigned int intorderadd_lb_ = 0)
: SubtractionDGLambda<PROBLEMDATA>(problem_, intorderadd_, intorderadd_lb_, weights_)
: SubtractionDGLambda<PROBLEMDATA, PenaltyFluxWeighting>(problem_, weighting_, intorderadd_,
intorderadd_lb_)
{
}
};
......
......@@ -7,6 +7,7 @@
#include <dune/pdelab/boilerplate/pdelab.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/eeg/source_model_interface.hh>
#include <duneuro/eeg/subtraction_dg_default_parameter.hh>
#include <duneuro/eeg/subtraction_dg_operator.hh>
......@@ -22,7 +23,8 @@ namespace duneuro
using Problem = SubtractionDGDefaultParameter<typename FS::GFS::Traits::GridViewType,
typename V::field_type, VC>;
using EdgeNormProvider = MultiEdgeNormProvider;
using LOP = SubtractionDG<Problem, EdgeNormProvider, continuityType>;
using PenaltyFluxWeighting = FittedDynamicPenaltyFluxWeights;
using LOP = SubtractionDG<Problem, EdgeNormProvider, PenaltyFluxWeighting, continuityType>;
using DOF = typename FS::DOF;
using AS = Dune::PDELab::GalerkinGlobalAssembler<FS, LOP, Dune::SolverCategory::sequential>;
using ElementType = typename BaseT::ElementType;
......@@ -36,10 +38,9 @@ namespace duneuro
: BaseT(search)
, problem_(volumeConductor->gridView(), volumeConductor)
, 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"))
, weighting_(solverConfig.get<std::string>("weighting", "tensorOnly"))
, lop_(problem_, weighting_, config.get<unsigned int>("intorderadd"),
config.get<unsigned int>("intorderadd_lb"))
, x_(fs.getGFS(), 0.0)
, res_(fs.getGFS(), 0.0)
, interp_(fs.getGFS(), 0.0)
......@@ -84,6 +85,7 @@ namespace duneuro
private:
Problem problem_;
EdgeNormProvider edgeNormProvider_;
PenaltyFluxWeighting weighting_;
LOP lop_;
mutable DOF x_;
mutable DOF res_;
......
......@@ -10,6 +10,7 @@
#include <duneuro/common/convection_diffusion_dg_operator.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/penalty_flux_weighting.hh>
#include <duneuro/tes/tdcs_patch_udg_parameter.hh>
namespace duneuro
......@@ -25,7 +26,9 @@ namespace duneuro
static const int dimension = FundamentalGridView::dimension;
using Problem = TDCSPatchUDGParameter<FundamentalGridView>;
using EdgeNormProvider = MultiEdgeNormProvider;
using LocalOperator = ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider>;
using PenaltyFluxWeighting = UnfittedDynamicPenaltyFluxWeights;
using LocalOperator =
ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider, PenaltyFluxWeighting>;
using WrappedLocalOperator = Dune::UDG::MultiPhaseLocalOperatorWrapper<LocalOperator>;
using UnfittedSubTriangulation = Dune::PDELab::UnfittedSubTriangulation<FundamentalGridView>;
using MatrixBackend = Dune::PDELab::istl::BCRSMatrixBackend<>;
......@@ -37,9 +40,10 @@ namespace duneuro
TDCSPatchUDGAssembler(Problem& problem, std::shared_ptr<ST> subTriangulation, const FS& fs,
const Dune::ParameterTree& config)
: edgeNormProvider_(config.get<std::string>("edge_norm_type"), 1.0)
, weighting_(config.get<std::string>("weighting"))
, lop_(problem, edgeNormProvider_,
ConvectionDiffusion_DG_Scheme::fromString(config.get<std::string>("scheme")),
ConvectionDiffusion_DG_Weights::weightsOn, config.get<double>("penalty"))
weighting_, config.get<double>("penalty"))
, wlop_(lop_)
, ust_(subTriangulation->gridView(), *subTriangulation)
, go_(fs.getGFS(), fs.getGFS(), ust_, wlop_, MatrixBackend(2 * dimension + 1))
......@@ -56,6 +60,7 @@ namespace duneuro
private:
EdgeNormProvider edgeNormProvider_;
PenaltyFluxWeighting weighting_;
LocalOperator lop_;
WrappedLocalOperator wlop_;
UnfittedSubTriangulation ust_;
......
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