Commit 2bdf3a0d authored by Santiago Ospina's avatar Santiago Ospina

Working at least on cubes and low order schemes

parent 12d40e0d
......@@ -181,37 +181,29 @@ public:
template<typename EG, typename LFSWC>
void onBindLFSV(const EG & eg, const LFSWC & lfsw_cache)
{
lfsv_volume.bind( eg.entity() );
std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
global_pl_view.bind(lfsw_cache);
pl.resize(lfsw_cache.size());
lfsv_volume.bind( eg.entity() );
rl.assign(lfsv_volume.size(),0.0);
}
template<typename IG, typename LFSUC, typename LFSWC>
void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSWC & lfsw_cache)
{
// This is never called
}
template<typename IG, typename LFSUC, typename LFSWC>
void onBindLFSUVOutside(const IG & ig,
const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
{
lfsv_skeleton.bind( ig.intersection() );
lfsv_skeleton_n.bind( ig.intersection() );
// inside part
std::cout << "lfsu_s_cache.size(): " << lfsu_s_cache.size() << std::endl;
std::cout << "lfsu_s.debug():" << std::endl; lfsu_s_cache.localFunctionSpace().debug();
std::cout << "lfsw_s_cache.size(): " << lfsw_s_cache.size() << std::endl;
std::cout << "lfsw_s.debug():" << std::endl; lfsw_s_cache.localFunctionSpace().debug();
std::cout << "lfsu_cache.size(): " << lfsu_cache.size() << std::endl;
std::cout << "lfsu.debug():" << std::endl; lfsu_cache.localFunctionSpace().debug();
std::cout << "lfsv_volume.size(): " << lfsv_volume.size() << std::endl;
std::cout << "lfsv_volume.debug():" << std::endl; lfsv_volume.debug();
std::cout << "lfsw_cache.size(): " << lfsw_cache.size() << std::endl;
std::cout << "lfsw.debug():" << std::endl; lfsw_cache.localFunctionSpace().debug();
std::cout << "lfsv_skeleton.size(): " << lfsv_skeleton.size() << std::endl;
std::cout << "lfsv_skeleton.debug():" << std::endl; lfsv_skeleton.debug();
......@@ -219,16 +211,18 @@ public:
std::cout << "ig.entity().subEntities(1): " << ig.inside().subEntities(1) << std::endl;
// Check that local matrix will be conforming
assert(lfsu_s_cache.size() == lfsw_s_cache.size());
assert(lfsu_s_cache.size() == (lfsv_volume.size() + lfsv_skeleton.size() * ig.inside().subEntities(1)));
global_sl_view.bind(lfsu_s_cache);
xl.assign(lfsu_s_cache.size(),0.0);
assert(lfsu_cache.size() == lfsw_cache.size());
}
global_pl_view.bind(lfsw_s_cache);
pl.resize(lfsw_s_cache.size());
template<typename IG, typename LFSUC, typename LFSWC>
void onBindLFSUVOutside(const IG & ig,
const LFSUC & lfsu_s_cache, const LFSWC & lfsw_s_cache,
const LFSUC & lfsu_n_cache, const LFSWC & lfsw_n_cache)
{
lfsv_skeleton_n.bind( ig.intersection() ); // TODO switch inside/outside
// outside part
// Check that local matrix will be conforming
assert(lfsu_s_cache.size() == (lfsv_volume.size() + lfsv_skeleton.size() * ig.inside().subEntities(1)));
assert(lfsv_skeleton.size() == lfsv_skeleton_n.size());
global_sn_view.bind(lfsu_n_cache);
......@@ -241,7 +235,8 @@ public:
template<typename IG, typename LFSWC>
void onBindLFSVInside(const IG & ig, const LFSWC & lfsw_cache)
{
// This is never called
// global_rl_view.bind(lfsw_s_cache);
rl.assign(lfsw_cache.size(),0.0);
}
template<typename IG, typename LFSWC>
......@@ -249,8 +244,6 @@ public:
const LFSWC & lfsw_s_cache,
const LFSWC & lfsw_n_cache)
{
// global_rl_view.bind(lfsw_s_cache);
rl.assign(lfsw_s_cache.size(),0.0);
// global_rn_view.bind(lfsw_n_cache);
rn.assign(lfsw_n_cache.size(),0.0);
}
......@@ -262,11 +255,19 @@ public:
const LFSUC & lfsu_cache,
const LFSVC & lfsv_cache)
{
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
x_vec.resize(lfsu_cache.size(),0.0);
m_matrix.solve(x_vec,r_vec);
auto& lfsu = lfsu_cache.localFunctionSpace();
// for (unsigned int i = 0; i < lfsu.size(); ++i)
// xl(lfsu,i) = x_vec[i];
// global_sl_view.add(xl);
for (unsigned int i = 0; i < lfsu.size(); ++i)
xl(lfsu,i) = x_vec[i];
global_sl_view.write(xl);
}
//! @}
......@@ -381,8 +382,15 @@ public:
for (unsigned int i = 0; i<lfsv_volume.size(); i++)
for (unsigned int j = 0; j<lfsu.size(); j++)
m_matrix[i][j] = (phiu[j]*phiv[i])/detB * factor/detB;
m_matrix[i][j] += (phiu[j]*phiv[i])/detB * factor/detB;
}
std::cout << "assembleUVVolume" << std::endl;
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
}
template<typename EG, typename LFSWC>
......@@ -393,6 +401,9 @@ public:
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)
{
std::cout << "assembleUVSkeleton" << std::endl;
rl_view.setWeight(local_assembler.weight());
rn_view.setWeight(local_assembler.weight());
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
......@@ -404,8 +415,10 @@ public:
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
r_vec[i+offset] += rl(lfsv_skeleton,i) + rn(lfsv_skeleton_n,i);
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << " + " << rn(lfsv_skeleton_n,i) << std::endl;
r_vec[i+offset] += rl(lfsv_skeleton,i);// + rn(lfsv_skeleton_n,i);
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
......@@ -451,10 +464,16 @@ public:
for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
for (unsigned int j = 0; j<lfsu_s.size(); j++){
auto xx = y*phiu_s[j]/detB * phiv_s[i+offset] * factor;
m_matrix[i+offset][j] += xx;
m_matrix[i+offset][j] += y*phiu_s[j]/detB * phiv_s[i] * factor;
}
}
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
}
template<typename IG, typename LFSWC>
......@@ -468,10 +487,15 @@ public:
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_skeleton,rl_view);
std::cout << "assembleUVBoundary" << std::endl;
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i)
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << std::endl;
r_vec[i+offset] += rl(lfsv_skeleton,i);
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
......@@ -517,8 +541,13 @@ public:
for (unsigned int i = 0; i<lfsv_skeleton.size(); i++)
for (unsigned int j = 0; j<lfsu_s.size(); j++)
m_matrix[i+offset][j] += (y*phiu_s[j])/detB * phiv_s[i+offset] * factor;
m_matrix[i+offset][j] += (y*phiu_s[j])/detB * phiv_s[i] * factor;
}
std::cout << "r_vec" << std::endl;
std::cout << r_vec << std::endl;
std::cout << "m_matrix" << std::endl;
std::cout << m_matrix << std::endl;
}
template<typename IG, typename LFSWC>
......@@ -549,9 +578,6 @@ public:
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);
x_vec.resize(lfsu_cache.size(),0.0);
m_matrix.solve(x_vec,r_vec);
}
template<typename EG, typename LFSWC>
......
......@@ -284,6 +284,9 @@ public:
std::vector<Vector> gradphiv;
// discrete gradient sign.
double dg_sign = 1.;
if constexpr (dimRangeTest == dim) {
// (we assume lfsv = RTk)
lfsv.finiteElement().localBasis().evaluateFunction(p_local,gradphiv);
......@@ -295,6 +298,7 @@ public:
y /= detB;
BTransposed.mtv(y,gradphiv[i]);
}
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
gradphiv = gradphiu;
......@@ -316,7 +320,7 @@ public:
// update residual
for (unsigned int i = 0; i<lfsv.size(); i++)
{
r.accumulate(lfsv,i, cond * (gradu*gradphiv[i]) * factor);
r.accumulate(lfsv,i, dg_sign * cond * (gradu*gradphiv[i]) * factor);
}
}
}
......@@ -401,11 +405,15 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
std::vector<Scalar> phiv_n(lfsv_n.size());
// discrete gradient sign.
double dg_sign = 1.;
// evaluate basis functions of test function
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);
// discrete gradient sign.
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
......@@ -543,23 +551,23 @@ public:
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_s * satCond_s * jump * factor);
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (tgradphi_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, theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (tgradphi_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, theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (tgradphi_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, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (tgradphi_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, theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
r_s.accumulate(lfsv_s, i, dg_sign * theta * (tgradphi_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, theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
r_n.accumulate(lfsv_n, i, dg_sign * theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
}
}
}
......@@ -619,10 +627,13 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
// discrete gradient sign.
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
......@@ -819,7 +830,7 @@ public:
// (non-)symmetric IP term
// symmetry term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
r_s.accumulate(lfsv_s,i, dg_sign * theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
continue;
}
......
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