Commit 2ce71e92 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

[Flux Reconstruction] Fix constributions for each test function in Richards operator

parent 9a8273c2
......@@ -106,7 +106,7 @@ class RaviartThomasFluxReconstruction
SkeletonFiniteElementMap<DF,RF,dim,order,gt>;
using GFSVSkeleton = Dune::PDELab::GridFunctionSpace<ES,SkeletonFEM>;
using LocalRaviartThomasEngine = Dune::Dorie::
LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton,true>;
LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton,false>;
// using Assember = typename GOP::Assembler;
using Assember = Dune::Dorie::DefaultAssembler<GFSU,GFSW,
Dune::PDELab::EmptyTransformation,
......
......@@ -179,6 +179,8 @@ protected:
std::vector<Cache> cache; //!< Local basis cache
enum class LOPCase {DG,RTVolume,RTSkeleton}; //!< Local operator case
public:
// pattern assembly flags
......@@ -263,6 +265,7 @@ public:
* the ansatz space.*/
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
constexpr LOPCase lopcase = (dimRangeTest==dim) ? LOPCase::RTVolume : LOPCase::DG;
const int order = lfsu.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * order;
......@@ -303,7 +306,7 @@ public:
// discrete gradient sign.
double dg_sign = 1.;
if constexpr (dimRangeTest == dim) {
if constexpr (lopcase == LOPCase::RTVolume) {
// (we assume lfsv = RTk)
lfsv.finiteElement().localBasis().evaluateFunction(p_local,gradphiv);
// Piola transformation
......@@ -315,9 +318,12 @@ public:
BTransposed.mtv(y,gradphiv[i]);
}
dg_sign = -1.;
} else {
} else if constexpr (lopcase == LOPCase::DG) {
// (we assume Galerkin method lfsu = lfsv)
gradphiv = gradphiu;
} else {
static_assert(Dune::AlwaysTrue<EG>::value,
"alpha_volume() does not support RTSkeleton case!");
}
// compute gradient of u
......@@ -334,10 +340,9 @@ public:
const RF factor = it.weight() * gt.integrationElement(p_local);
// update residual
for (unsigned int i = 0; i<lfsv.size(); i++)
{
r.accumulate(lfsv,i, dg_sign * cond * (gradu*gradphiv[i]) * factor);
}
if constexpr (lopcase != LOPCase::RTSkeleton)
for (unsigned int i = 0; i<lfsv.size(); i++)
r.accumulate(lfsv,i, dg_sign * cond * (gradu*gradphiv[i]) * factor);
}
}
......@@ -378,6 +383,8 @@ public:
constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
constexpr LOPCase lopcase = (dimDomainLocalTest==dim-1) ? LOPCase::RTSkeleton :
(dimRangeTest==dim ? LOPCase::RTVolume : LOPCase::DG);
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
......@@ -444,7 +451,7 @@ public:
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
if constexpr (lopcase == LOPCase::RTSkeleton) {
// Evaluate position in local coordinates of the inside entity of codim 1
using LocalCoordinate = typename IG::LocalGeometry::LocalCoordinate;
LocalCoordinate sub_p_local_s, sub_p_local_n;
......@@ -467,8 +474,8 @@ public:
lfsv_n.finiteElement().localBasis().evaluateFunction(sub_p_local_n,phiv_n);
// discrete gradient sign.
dg_sign = 0.;
} else if constexpr (dimRangeTest == dim) {
dg_sign = -1.;
} else if constexpr (lopcase == LOPCase::RTVolume) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(p_local_n,gradphiv_n);
......@@ -581,7 +588,7 @@ public:
relCond_upwind = cond_upwind / satCond_upwind;
}
if constexpr (dimRangeTest != dim)
if constexpr (lopcase != LOPCase::RTVolume)
{
// update residual (flux term)
// diffusion term
......@@ -590,52 +597,52 @@ public:
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * numFlux * phiv_s[i] * factor);
r_s.accumulate(lfsv_s, i, numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - dg_sign * numFlux * phiv_n[i] * factor);
r_n.accumulate(lfsv_n, i, - numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * relCond_upwind * numFlux * phiv_s[i] * factor);
r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, - dg_sign * relCond_upwind * numFlux * phiv_n[i] * factor);
r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * relCond_upwind * numFlux_upwind * phiv_s[i] * factor);
r_s.accumulate(lfsv_s, i, relCond_upwind * numFlux_upwind * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - dg_sign * relCond_upwind * numFlux_upwind * phiv_n[i] * factor);
r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phiv_n[i] * factor);
}
}
if constexpr (dimDomainLocalTest == dim-1)
return;
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_s * satCond_s * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_n * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
if constexpr (lopcase != LOPCase::RTSkeleton)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * (gradphiv_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_s * satCond_s * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_n * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * (gradphiv_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
}
}
}
}
......@@ -660,6 +667,8 @@ public:
constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
constexpr LOPCase lopcase = (dimDomainLocalTest==dim-1) ? LOPCase::RTSkeleton :
(dimRangeTest==dim ? LOPCase::RTVolume : LOPCase::DG);
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
......@@ -691,25 +700,25 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto js = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> gradphiu_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
js.mv(jsu_s[i][0],gradphiu_s[i]);
// transform gradients to real element
const auto js = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> gradphiu_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
js.mv(jsu_s[i][0],gradphiu_s[i]);
std::vector<Vector> gradphiv_s(lfsv_s.size());
// discrete gradient sign.
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
if constexpr (lopcase == LOPCase::RTSkeleton) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
dg_sign = 0.;
} else if constexpr (dimRangeTest == dim) {
dg_sign = -1.;
} else if constexpr (lopcase == LOPCase::RTVolume) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
// Piola transformation
......@@ -721,7 +730,6 @@ public:
BTransposed_s.mtv(y_s,gradphiv_s[i]);
}
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
......@@ -765,8 +773,9 @@ public:
normal_flux *= std::abs( ig.centerUnitOuterNormal() * param.gravity() );
// update residual (flux term)
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * normal_flux * phiv_s[i] * factor);
if constexpr (lopcase != LOPCase::RTVolume)
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, normal_flux * phiv_s[i] * factor);
continue;
}
......@@ -816,19 +825,25 @@ public:
relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
}
// update residual (flux term)
// diffusion term
// consistency term
// + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * relCond * numFlux * phiv_s[i] * factor);
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond * satCond_s * jump * factor);
if constexpr (lopcase != LOPCase::RTVolume)
{
// update residual (flux term)
// diffusion term
// consistency term
// + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, relCond * numFlux * phiv_s[i] * factor);
}
if constexpr (lopcase != LOPCase::RTSkeleton)
{
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond * satCond_s * jump * factor);
}
continue;
}
......
Supports Markdown
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