Commit bdf77348 authored by Santiago Ospina's avatar Santiago Ospina

Allow non-conforming elements for raviart thomas reconstruction in the richards local operator

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 6d7303aa
......@@ -159,8 +159,8 @@ public:
//! Called immediately after binding of local function space in
//! global assembler.
//! @{
template<typename EG, typename LFSUC, typename LFSVC>
void onBindLFSUV(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsw_cache)
template<typename EG, typename LFSUC, typename LFSWC>
void onBindLFSUV(const EG & eg, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
{
// Check that local matrix will be conforming
assert(lfsu_cache.size() == lfsw_cache.size());
......@@ -170,8 +170,8 @@ public:
xl.resize(lfsu_cache.size());
}
template<typename EG, typename LFSVC>
void onBindLFSV(const EG & eg, const LFSVC & lfsw_cache)
template<typename EG, typename LFSWC>
void onBindLFSV(const EG & eg, const LFSWC & lfsw_cache)
{
global_pl_view.bind(lfsw_cache);
pl.assign(lfsw_cache.size(),0.0);
......@@ -184,36 +184,36 @@ public:
// TODO local matrix = 0
}
template<typename IG, typename LFSUC, typename LFSVC>
void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSVC & lfsw_cache)
template<typename IG, typename LFSUC, typename LFSWC>
void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
{
global_sl_view.bind(lfsu_cache);
xl.resize(lfsu_cache.size());
}
template<typename IG, typename LFSUC, typename LFSVC>
template<typename IG, typename LFSUC, typename LFSWC>
void onBindLFSUVOutside(const IG & ig,
const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache)
const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
{
global_sn_view.bind(lfsu_n_cache);
xn.resize(lfsu_n_cache.size());
}
template<typename IG, typename LFSVC>
void onBindLFSVInside(const IG & ig, const LFSVC & lfsw_cache)
template<typename IG, typename LFSWC>
void onBindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
{
global_pl_view.bind(lfsw_cache);
rl.assign(lfsw_cache.size(),0.0);
}
template<typename IG, typename LFSVC>
template<typename IG, typename LFSWC>
void onBindLFSVOutside(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache)
const LFSWC & lfsw_s_cache,
const LFSWC & lfsw_n_cache)
{
global_pn_view.bind(lfsv_n_cache);
rn.assign(lfsv_n_cache.size(),0.0);
global_pn_view.bind(lfsw_n_cache);
rn.assign(lfsw_n_cache.size(),0.0);
}
//! @}
......@@ -221,26 +221,26 @@ public:
//! Called when the local function space is about to be rebound or
//! discarded
//! @{
template<typename EG, typename LFSVC>
void onUnbindLFSV(const EG & eg, const LFSVC & lfsw_cache)
template<typename EG, typename LFSWC>
void onUnbindLFSV(const EG & eg, const LFSWC & lfsw_cache)
{
//FIXME
// global_pl_view.add(rl);
// global_pl_view.commit();
}
template<typename IG, typename LFSVC>
void onUnbindLFSVInside(const IG & ig, const LFSVC & lfsw_cache)
template<typename IG, typename LFSWC>
void onUnbindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
{
//FIXME
// global_pl_view.add(rl);
// global_pl_view.commit();
}
template<typename IG, typename LFSVC>
template<typename IG, typename LFSWC>
void onUnbindLFSVOutside(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache)
const LFSWC & lfsw_s_cache,
const LFSWC & lfsw_n_cache)
{
//FIXME
// global_pn_view.add(rn);
......@@ -295,93 +295,76 @@ public:
return LocalAssembler::isNonOverlapping && eg.entity().partitionType() != Dune::InteriorEntity;
}
template<typename EG, typename LFSUC, typename LFSVC>
void assembleUVVolume(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsw_cache)
template<typename EG, typename LFSUC, typename LFSWC>
void assembleUVVolume(const EG & eg, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
{
rl_view.setWeight(local_assembler.weight());
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolume>::
alpha_volume(lop,eg,lfsw_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
}
template<typename EG, typename LFSVC>
void assembleVVolume(const EG & eg, const LFSVC & lfsw_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaVolume>::
// lambda_volume(lop,eg,lfsw_cache.localFunctionSpace(),rl_view);
}
template<typename EG, typename LFSWC>
void assembleVVolume(const EG & eg, const LFSWC & lfsw_cache)
{}
template<typename IG, typename LFSUC, typename LFSVC>
void assembleUVSkeleton(const IG & ig, const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache)
template<typename IG, typename LFSUC, typename LFSWC>
void assembleUVSkeleton(const IG & ig, const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
{
// rl_view.setWeight(local_assembler.weight());
// rn_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
// alpha_skeleton(lop,ig,
// lfsu_s_cache.localFunctionSpace(),xl,lfsv_s_cache.localFunctionSpace(),
// lfsu_n_cache.localFunctionSpace(),xn,lfsv_n_cache.localFunctionSpace(),
// rl_view,rn_view);
rl_view.setWeight(local_assembler.weight());
rn_view.setWeight(local_assembler.weight());
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
alpha_skeleton(lop,ig,
lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,
lfsw_n_cache.localFunctionSpace(),pn,lfsv_skeleton_n,
rl_view,rn_view);
}
template<typename IG, typename LFSVC>
void assembleVSkeleton(const IG & ig, const LFSVC & lfsv_s_cache, const LFSVC & lfsv_n_cache)
{
// rl_view.setWeight(local_assembler.weight());
// rn_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaSkeleton>::
// lambda_skeleton(lop, ig, lfsv_s_cache.localFunctionSpace(), lfsv_n_cache.localFunctionSpace(), rl_view, rn_view);
}
template<typename IG, typename LFSWC>
void assembleVSkeleton(const IG & ig, const LFSWC & lfsw_s_cache, const LFSWC & lfsw_n_cache)
{}
template<typename IG, typename LFSUC, typename LFSVC>
void assembleUVBoundary(const IG & ig, const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache)
template<typename IG, typename LFSUC, typename LFSWC>
void assembleUVBoundary(const IG & ig, const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
// alpha_boundary(lop,ig,lfsu_s_cache.localFunctionSpace(),xl,lfsv_s_cache.localFunctionSpace(),rl_view);
rl_view.setWeight(local_assembler.weight());
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,rl_view);
}
template<typename IG, typename LFSVC>
void assembleVBoundary(const IG & ig, const LFSVC & lfsv_s_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaBoundary>::
// lambda_boundary(lop,ig,lfsv_s_cache.localFunctionSpace(),rl_view);
}
template<typename IG, typename LFSWC>
void assembleVBoundary(const IG & ig, const LFSWC & lfsw_s_cache)
{}
template<typename IG, typename LFSUC, typename LFSVC>
template<typename IG, typename LFSUC, typename LFSWC>
static void assembleUVEnrichedCoupling(const IG & ig,
const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache,
const LFSUC & lfsu_coupling_cache, const LFSVC & lfsv_coupling_cache)
const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache,
const LFSUC & lfsu_coupling_cache, const LFSWC & lfsw_coupling_cache)
{
DUNE_THROW(Dune::NotImplemented,"Assembling of coupling spaces is not implemented for ");
}
template<typename IG, typename LFSVC>
template<typename IG, typename LFSWC>
static void assembleVEnrichedCoupling(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache,
const LFSVC & lfsv_coupling_cache)
const LFSWC & lfsw_s_cache,
const LFSWC & lfsw_n_cache,
const LFSWC & lfsw_coupling_cache)
{
DUNE_THROW(Dune::NotImplemented,"Assembling of coupling spaces is not implemented for ");
}
template<typename EG, typename LFSUC, typename LFSVC>
void assembleUVVolumePostSkeleton(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsw_cache)
template<typename EG, typename LFSUC, typename LFSWC>
void assembleUVVolumePostSkeleton(const EG & eg, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
// alpha_volume_post_skeleton(lop,eg,lfsu_cache.localFunctionSpace(),xl,lfsw_cache.localFunctionSpace(),rl_view);
rl_view.setWeight(local_assembler.weight());
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
alpha_volume_post_skeleton(lop,eg,lfsw_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
}
template<typename EG, typename LFSVC>
void assembleVVolumePostSkeleton(const EG & eg, const LFSVC & lfsw_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaVolumePostSkeleton>::
// lambda_volume_post_skeleton(lop,eg,lfsw_cache.localFunctionSpace(),rl_view);
}
template<typename EG, typename LFSWC>
void assembleVVolumePostSkeleton(const EG & eg, const LFSWC & lfsw_cache)
{}
//! @}
......
......@@ -61,6 +61,9 @@ class RaviartThomasFluxReconstruction // : public GridFunctionBase<>
GV gv;
// const IndexSet& is;
GFSU _gfsu;
Range _x;
// // discretization related stuff
// DF time; // this is just for set_time
// int intorderadd;
......@@ -87,13 +90,8 @@ public:
RaviartThomasFluxReconstruction(GO& go, int intorderadd_=2)
: _go(go)
, gv(go.trialGridFunctionSpace().gridView())
// , is(gv.indexSet())
// , time(0)
// , intorderadd(intorderadd_)
// , quadrature_factor(2)
// , lfe(RTFactory::get())
// , rtdof(is.size(0))
// , lastindex(-1)
, _gfsu(_go.trialGridFunctionSpace().entitySet(),gv)
, _x(_gfsu,0.0)
{}
void update (const Domain& p)
......@@ -110,12 +108,10 @@ public:
const GFSW& gfsw = _go.trialGridFunctionSpace();
GFSU gfsu(entity_set,gv);
auto& local_assembler_go = _go.localAssembler();
auto& local_operator = local_assembler_go.localOperator();
MBE mbe(0);
GOP gop(gfsu,gfsw,local_operator,mbe);
GOP gop(_gfsu,gfsw,local_operator,mbe);
const auto& global_assembler_gop = gop.assembler();
const auto& local_assembler_gop = gop.localAssembler();
......@@ -124,7 +120,7 @@ public:
LocalRaviartThomasEngine local_raviart_thomas_engine(local_assembler_gop,gfsv_volume,gfsv_skeleton);
local_raviart_thomas_engine.setPrescription(p);
// local_raviart_thomas_engine.setSolution(x);
local_raviart_thomas_engine.setSolution(_x);
global_assembler_gop.assemble(local_raviart_thomas_engine);
}
};
......
......@@ -292,6 +292,7 @@ public:
auto detB = BTransposed.determinant();
for (unsigned int i=0; i<lfsv.size(); i++) {
Vector y(gradphiv[i]);
y /= detB;
BTransposed.mtv(y,gradphiv[i]);
}
} else {
......@@ -346,6 +347,17 @@ public:
const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
R& r_s, R& r_n) const
{
// Switches between local and global interface in test local space
using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType>;
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
/* If domain dimension of the test local function space is dim we assume
* that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
* assume that it is a skeleton finite element function space
* which is non conforming with respect to the ansatz space.*/
constexpr int dimDomainLocalTest = BasisSwitchTest::dimDomainLocal;
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
const int order_n = lfsu_n.finiteElement().localBasis().order();
......@@ -390,9 +402,15 @@ public:
std::vector<Scalar> phiv_n(lfsv_n.size());
// evaluate basis functions of test function
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,phiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(p_local_n,phiv_n);
if constexpr (dimDomainLocalTest == dim-1) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(it.position(),phiv_n);
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
phiv_n = phiu_n;
}
// evaluate u
RF u_s = 0., u_n = 0.;
......@@ -555,6 +573,17 @@ public:
const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
R& r_s) const
{
// Switches between local and global interface in test local space
using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType>;
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
/* If domain dimension of the test local function space is dim we assume
* that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
* assume that it is a skeleton finite element function space
* which is non conforming with respect to the ansatz space.*/
constexpr int dimDomainLocalTest = BasisSwitchTest::dimDomainLocal;
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
const int degree = order_s;
......@@ -591,9 +620,13 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
// evaluate basis functions of test function
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,phiv_s);
if constexpr (dimDomainLocalTest == dim-1) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
}
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
......
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