Commit 175f8007 authored by Andreas Nüßing's avatar Andreas Nüßing

Merge branch 'feature/cutfem-eeg' into 'master'

[CutFEM] add first implementation of cutfem

See merge request duneuro/duneuro!47
parents 93a3472f e1055324
......@@ -48,13 +48,6 @@ namespace duneuro
}
};
/**
* \brief A struct for choosing the weighted/non-weighted variant of the DG scheme.
*/
struct ConvectionDiffusion_DG_Weights {
enum Type { weightsOn, weightsOff };
};
/**
* \brief Used when evaluating the diffusion tensor on an intersection
*
......@@ -162,16 +155,17 @@ namespace duneuro
#warning Assuming piecewise constant diffusion tensor at the spots marked by (*).
#warning Assuming continuity of the velocity field at the spots marked by (**) such that we may choose any side for its evaluation.
#warning UDG assembler: Skeleton/boundary integrals: We evaluate data functions on the inside/outside "host entities" not on the "fundamental mesh home entities".
template <typename T, typename EdgeNormProvider>
template <typename T, typename EdgeNormProvider, typename PenaltyFluxWeighting>
class ConvectionDiffusion_DG_LocalOperator
: public Dune::PDELab::
NumericalJacobianApplyVolume<ConvectionDiffusion_DG_LocalOperator<T, EdgeNormProvider>>,
NumericalJacobianApplyVolume<ConvectionDiffusion_DG_LocalOperator<T, EdgeNormProvider,
PenaltyFluxWeighting>>,
public Dune::PDELab::
NumericalJacobianApplySkeleton<ConvectionDiffusion_DG_LocalOperator<T,
EdgeNormProvider>>,
NumericalJacobianApplySkeleton<ConvectionDiffusion_DG_LocalOperator<T, EdgeNormProvider,
PenaltyFluxWeighting>>,
public Dune::PDELab::
NumericalJacobianApplyBoundary<ConvectionDiffusion_DG_LocalOperator<T,
EdgeNormProvider>>,
NumericalJacobianApplyBoundary<ConvectionDiffusion_DG_LocalOperator<T, EdgeNormProvider,
PenaltyFluxWeighting>>,
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::FullSkeletonPattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
......@@ -199,29 +193,30 @@ namespace duneuro
//! mesh
ConvectionDiffusion_DG_LocalOperator(
T& param_, const EdgeNormProvider& edgenormprovider_,
const PenaltyFluxWeighting& weighting_,
const ConvectionDiffusion_DG_Scheme::Type scheme_ = ConvectionDiffusion_DG_Scheme::NIPG,
const ConvectionDiffusion_DG_Weights::Type weights_ =
ConvectionDiffusion_DG_Weights::weightsOff,
Real alpha_ = 0.0, const bool useOutflowBoundaryConditionAndItsFluxOnInflow_ = false,
const int intorderadd_ = 0)
: Dune::PDELab::
NumericalJacobianApplyVolume<ConvectionDiffusion_DG_LocalOperator<T,
EdgeNormProvider>>(
NumericalJacobianApplyVolume<ConvectionDiffusion_DG_LocalOperator<T, EdgeNormProvider,
PenaltyFluxWeighting>>(
1.0e-7)
, Dune::PDELab::
NumericalJacobianApplySkeleton<ConvectionDiffusion_DG_LocalOperator<T,
EdgeNormProvider>>(
EdgeNormProvider,
PenaltyFluxWeighting>>(
1.0e-7)
, Dune::PDELab::
NumericalJacobianApplyBoundary<ConvectionDiffusion_DG_LocalOperator<T,
EdgeNormProvider>>(
EdgeNormProvider,
PenaltyFluxWeighting>>(
1.0e-7)
, param(param_)
, useOutflowBoundaryConditionAndItsFluxOnInflow(
useOutflowBoundaryConditionAndItsFluxOnInflow_)
, edgenormprovider(edgenormprovider_)
, weighting(weighting_)
, scheme(scheme_)
, weights(weights_)
, alpha(alpha_)
, intorderadd(intorderadd_)
, quadrature_factor(2)
......@@ -410,27 +405,14 @@ namespace duneuro
maxH = std::max(h_F, maxH);
assert(h_F > 1e-20);
// compute weights
RF omega_s;
RF omega_n;
RF harmonic_average(0.0);
if (weights == ConvectionDiffusion_DG_Weights::weightsOn) {
const RF delta_s = (An_F_s * n_F);
const 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);
// get polynomial degree
const int degree = std::max(order_s, order_n);
// penalty factor
const RF penalty_factor = (alpha / h_F) * harmonic_average * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * harmonic_average;
const RF penalty_factor = (alpha / h_F) * weights.penaltyWeight * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * weights.penaltyWeight;
// create copies of inside and outside entities
const auto& outsideEntity = ig.outside();
......@@ -484,7 +466,9 @@ namespace duneuro
const RF factor = qp.weight() * ig.geometry().integrationElement(qp.position());
// diffusion term
const RF term2 = -(omega_s * (An_F_s * gradu_s) + omega_n * (An_F_n * gradu_n)) * factor;
const RF term2 = -(weights.fluxInsideWeight * (An_F_s * gradu_s)
+ weights.fluxOutsideWeight * (An_F_n * gradu_n))
* factor;
for (size_type i = 0; i < lfsu_s.size(); i++)
r_s.accumulate(lfsu_s, i, term2 * phi_s[i]);
for (size_type i = 0; i < lfsu_n.size(); i++)
......@@ -493,11 +477,15 @@ namespace duneuro
// (non-)symmetric IP term
const RF term3 = (u_s - u_n) * factor;
for (size_type i = 0; i < lfsu_s.size(); i++)
// r_s.accumulate(lfsu_s,i,term3 * theta * omega_s * (An_F_s*tgradphi_s[i]));
r_s.accumulate(lfsu_s, i, term3 * theta * omega_s * (An_F_s * gradphi_s[i][0]));
// r_s.accumulate(lfsu_s,i,term3 * theta * weights.fluxInsideWeight *
// (An_F_s*tgradphi_s[i]));
r_s.accumulate(lfsu_s, i,
term3 * theta * weights.fluxInsideWeight * (An_F_s * gradphi_s[i][0]));
for (size_type i = 0; i < lfsu_n.size(); i++)
// r_n.accumulate(lfsu_n,i,term3 * theta * omega_n * (An_F_n*tgradphi_n[i]));
r_n.accumulate(lfsu_n, i, term3 * theta * omega_n * (An_F_n * gradphi_n[i][0]));
// r_n.accumulate(lfsu_n,i,term3 * theta * weights.fluxOutsideWeight *
// (An_F_n*tgradphi_n[i]));
r_n.accumulate(lfsu_n, i,
term3 * theta * weights.fluxOutsideWeight * (An_F_n * gradphi_n[i][0]));
// standard IP term integral
const RF term4 = penalty_factor * (u_s - u_n) * factor;
......@@ -562,26 +550,14 @@ namespace duneuro
assert(h_F > 1e-20);
// compute weights
RF omega_s;
RF omega_n;
RF harmonic_average(0.0);
if (weights == ConvectionDiffusion_DG_Weights::weightsOn) {
const RF delta_s = (An_F_s * n_F);
const 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);
// get polynomial degree
const int degree = std::max(order_s, order_n);
// penalty factor
const RF penalty_factor = (alpha / h_F) * harmonic_average * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * harmonic_average;
const RF penalty_factor = (alpha / h_F) * weights.penaltyWeight * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * weights.penaltyWeight;
// create copies of inside and outside entities
const auto& insideEntity = ig.inside();
......@@ -619,12 +595,13 @@ namespace duneuro
// do all terms in the order: I convection, II diffusion, III consistency, IV ip
for (size_type j = 0; j < lfsu_s.size(); j++) {
// const RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
const RF temp1 = -(An_F_s * gradphi_s[j][0]) * omega_s * factor;
// const RF temp1 = -(An_F_s*tgradphi_s[j])*weights.fluxInsideWeight*factor;
const RF temp1 = -(An_F_s * gradphi_s[j][0]) * weights.fluxInsideWeight * factor;
for (size_type i = 0; i < lfsu_s.size(); i++) {
mat_ss.accumulate(lfsu_s, i, lfsu_s, j, temp1 * phi_s[i]);
mat_ss.accumulate(lfsu_s, i, lfsu_s, j,
phi_s[j] * factor * theta * omega_s * (An_F_s * gradphi_s[i][0]));
mat_ss.accumulate(lfsu_s, i, lfsu_s, j, phi_s[j] * factor * theta
* weights.fluxInsideWeight
* (An_F_s * gradphi_s[i][0]));
mat_ss.accumulate(lfsu_s, i, lfsu_s, j, phi_s[j] * ipfactor * phi_s[i]);
if (std::isnan(mat_ss.container()(lfsu_s, i, lfsu_s, j))) {
for (int k = 0; k < ig.geometry().corners(); ++k) {
......@@ -658,11 +635,12 @@ namespace duneuro
}
}
for (size_type j = 0; j < lfsu_n.size(); j++) {
const RF temp1 = -(An_F_n * gradphi_n[j][0]) * omega_n * factor;
const RF temp1 = -(An_F_n * gradphi_n[j][0]) * weights.fluxOutsideWeight * factor;
for (size_type i = 0; i < lfsu_s.size(); i++) {
mat_sn.accumulate(lfsu_s, i, lfsu_n, j, temp1 * phi_s[i]);
mat_sn.accumulate(lfsu_s, i, lfsu_n, j,
-phi_n[j] * factor * theta * omega_s * (An_F_s * gradphi_s[i][0]));
mat_sn.accumulate(lfsu_s, i, lfsu_n, j, -phi_n[j] * factor * theta
* weights.fluxInsideWeight
* (An_F_s * gradphi_s[i][0]));
mat_sn.accumulate(lfsu_s, i, lfsu_n, j, -phi_n[j] * ipfactor * phi_s[i]);
if (std::isnan(mat_sn.container()(lfsu_s, i, lfsu_n, j))) {
DUNE_THROW(Dune::Exception, "NAN found");
......@@ -670,11 +648,12 @@ namespace duneuro
}
}
for (size_type j = 0; j < lfsu_s.size(); j++) {
const RF temp1 = -(An_F_s * gradphi_s[j][0]) * omega_s * factor;
const RF temp1 = -(An_F_s * gradphi_s[j][0]) * weights.fluxInsideWeight * factor;
for (size_type i = 0; i < lfsu_n.size(); i++) {
mat_ns.accumulate(lfsu_n, i, lfsu_s, j, -temp1 * phi_n[i]);
mat_ns.accumulate(lfsu_n, i, lfsu_s, j,
phi_s[j] * factor * theta * omega_n * (An_F_n * gradphi_n[i][0]));
mat_ns.accumulate(lfsu_n, i, lfsu_s, j, phi_s[j] * factor * theta
* weights.fluxOutsideWeight
* (An_F_n * gradphi_n[i][0]));
mat_ns.accumulate(lfsu_n, i, lfsu_s, j, -phi_s[j] * ipfactor * phi_n[i]);
if (std::isnan(mat_ns.container()(lfsu_n, i, lfsu_s, j))) {
DUNE_THROW(Dune::Exception, "NAN found");
......@@ -682,11 +661,12 @@ namespace duneuro
}
}
for (size_type j = 0; j < lfsu_n.size(); j++) {
const RF temp1 = -(An_F_n * gradphi_n[j][0]) * omega_n * factor;
const RF temp1 = -(An_F_n * gradphi_n[j][0]) * weights.fluxOutsideWeight * factor;
for (size_type i = 0; i < lfsu_n.size(); i++) {
mat_nn.accumulate(lfsu_n, i, lfsu_n, j, -temp1 * phi_n[i]);
mat_nn.accumulate(lfsu_n, i, lfsu_n, j,
-phi_n[j] * factor * theta * omega_n * (An_F_n * gradphi_n[i][0]));
mat_nn.accumulate(lfsu_n, i, lfsu_n, j, -phi_n[j] * factor * theta
* weights.fluxOutsideWeight
* (An_F_n * gradphi_n[i][0]));
mat_nn.accumulate(lfsu_n, i, lfsu_n, j, phi_n[j] * ipfactor * phi_n[i]);
if (std::isnan(mat_nn.container()(lfsu_n, i, lfsu_n, j))) {
DUNE_THROW(Dune::Exception, "NAN found");
......@@ -736,10 +716,6 @@ namespace duneuro
Dune::FieldVector<RF, dim> An_F_s;
A_s.mv(n_F, An_F_s);
// evaluate boundary condition type (see also (***))
const Dune::FieldVector<DF, dim - 1> face_local =
Dune::ReferenceElements<DF, dim - 1>::general(gtface).position(0, 0);
// face diameter
RF h_F;
edgenormprovider.edgeNorm(ig, h_F, true);
......@@ -748,18 +724,14 @@ namespace duneuro
assert(h_F > 1e-20);
// compute weights
RF harmonic_average;
if (weights == ConvectionDiffusion_DG_Weights::weightsOn)
harmonic_average = An_F_s * n_F;
else
harmonic_average = 1.0;
auto weights = weighting(ig, A_s, A_s);
// get polynomial degree
const int degree = order_s;
// penalty factor
const RF penalty_factor = (alpha / h_F) * harmonic_average * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * harmonic_average;
const RF penalty_factor = (alpha / h_F) * weights.penaltyWeight * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * weights.penaltyWeight;
// create copy of inside Entity
const auto& insideEntity = ig.inside();
......@@ -887,18 +859,14 @@ namespace duneuro
assert(h_F > 1e-20);
// compute weights
RF harmonic_average;
if (weights == ConvectionDiffusion_DG_Weights::weightsOn)
harmonic_average = An_F_s * n_F;
else
harmonic_average = 1.0;
auto weights = weighting(ig, A_s, A_s);
// get polynomial degree
const int degree = order_s;
// penalty factor
const RF penalty_factor = (alpha / h_F) * harmonic_average * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * harmonic_average;
const RF penalty_factor = (alpha / h_F) * weights.penaltyWeight * degree * (degree + dim - 1);
// const RF penalty_factor = (alpha/h_F) * weights.penaltyWeight;
// create copy of inside entity
const auto& insideEntity = ig.inside();
......@@ -1026,8 +994,8 @@ namespace duneuro
// DG scheme related parameters
EdgeNormProvider edgenormprovider;
const PenaltyFluxWeighting weighting;
const ConvectionDiffusion_DG_Scheme::Type scheme;
const ConvectionDiffusion_DG_Weights::Type weights;
Real alpha;
Real theta;
......
// -*- 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/penalty_flux_weighting.hh>
#include <duneuro/common/random.hh>
#include <duneuro/common/vector_initialization.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 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>;
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)
, weighting_(config.get<std::string>("weights"))
, 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(),
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;
initialize(Dune::PDELab::Backend::native(solution), config.hasSub("initialization") ?
config.sub("initialization") :
Dune::ParameterTree());
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_;
}
#if HAVE_TBB
tbb::mutex& functionSpaceMutex()
{
return fsMutex_;
}
#endif
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::PenaltyFluxWeighting weighting_;
typename Traits::LocalOperator localOperator_;
typename Traits::WrappedLocalOperator wrappedLocalOperator_;
typename Traits::UnfittedSubTriangulation unfittedSubTriangulation_;
typename Traits::RawGridOperator rawGridOperator_;
typename Traits::GridOperator gridOperator_;
typename Traits::LinearSolver linearSolver_;
#if HAVE_TBB
tbb::mutex fsMutex_;
#endif
};
}
#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
#include <duneuro/common/deprecated.hh>
#include <iostream>
namespace duneuro
{
void issueDeprecationWarning(const std::string& msg)
{
std::cout << "DEPRECATION WARNING:\n" << msg << std::endl;
}
}
#ifndef DUNEURO_DEPRECATED_HH
#define DUNEURO_DEPRECATED_HH
#include <string>
namespace duneuro
{
void issueDeprecationWarning(const std::string& msg);
}
#endif // DUNEURO_DEPRECATED_HH
......@@ -11,6 +11,7 @@
#include <duneuro/common/edge_norm_provider.hh>