Commit c7b3252a authored by Santiago Ospina's avatar Santiago Ospina

Update flux reconstruction

* Added available higher orders
* Fixed sign for lifting
* Enabled reconstruction for adaptivity cases
* Fixed bug in local keys of volume raviart thomas
Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 5cccbd55
......@@ -136,18 +136,18 @@ int main(int argc, char** argv)
if (gtype == "gmsh"){
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<2>>(inifile,helper);
switch(FEorder){
// case 1:{
// Sim<Simplex<2,1>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{
// Sim<Simplex<2,2>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 1:{
Sim<Simplex<2,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<Simplex<2,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 3:{ *
// Sim<Simplex<2,3>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
......@@ -162,24 +162,24 @@ int main(int argc, char** argv)
if(adaptivity){
auto grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<2>>(inifile,helper);
switch(FEorder){
// case 1:{
// Sim<CubeAdaptive<2,1>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{
// Sim<CubeAdaptive<2,2>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 3:{
// Sim<CubeAdaptive<2,3>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 1:{
Sim<CubeAdaptive<2,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<CubeAdaptive<2,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<CubeAdaptive<2,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -192,16 +192,16 @@ int main(int argc, char** argv)
sim.run();
break;
}
// case 2:{
// Sim<Cube<2,2>> sim(helper,grid,inifile);
// sim.run();
// break;
// }
// case 3:{
// Sim<Cube<2,3>> sim(helper,grid,inifile);
// sim.run();
// break;
// }
case 2:{
Sim<Cube<2,2>> sim(helper,grid,inifile);
sim.run();
break;
}
case 3:{
Sim<Cube<2,3>> sim(helper,grid,inifile);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -242,18 +242,18 @@ int main(int argc, char** argv)
if(adaptivity){
auto grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<3>>(inifile,helper);
switch(FEorder){
// case 1:{
// Sim<CubeAdaptive<3,1>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{
// Sim<CubeAdaptive<3,2>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
case 1:{
Sim<CubeAdaptive<3,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<CubeAdaptive<3,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 3:{ *
// Sim<CubeAdaptive<3,3>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
......@@ -267,16 +267,16 @@ int main(int argc, char** argv)
else{ // no adaptivity
auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<3>>(inifile,helper);
switch(FEorder){
// case 1:{
// Sim<Cube<3,1>> sim(helper,grid,inifile);
// sim.run();
// break;
// }
// case 2:{
// Sim<Cube<3,2>> sim(helper,grid,inifile);
// sim.run();
// break;
// }
case 1:{
Sim<Cube<3,1>> sim(helper,grid,inifile);
sim.run();
break;
}
case 2:{
Sim<Cube<3,2>> sim(helper,grid,inifile);
sim.run();
break;
}
// case 3:{ *
// Sim<Cube<3,3>> sim(helper,grid,inifile);
// sim.run();
......
......@@ -8,8 +8,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,8 +8,8 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>; *
} // namespace Dorie
......
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1>>; *
} // namespace Dorie
......
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>; *
} // namespace Dorie
......
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,Geo::cube>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,1>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,2>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -209,9 +209,6 @@ public:
std::cout << "lfsv_skeleton.debug():" << std::endl; lfsv_skeleton.debug();
std::cout << "ig.entity().subEntities(1): " << ig.inside().subEntities(1) << std::endl;
// Check that local matrix will be conforming
assert(lfsu_cache.size() == lfsw_cache.size());
}
template<typename IG, typename LFSUC, typename LFSWC>
......@@ -404,6 +401,8 @@ 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)
{
// if (not ig.intersection().conforming())
// std::cout << WARNING: Non-conforming entities produce non-continuos fluxes in normal direction of the entity << std::endl;
std::cout << "assembleUVSkeleton" << std::endl;
......@@ -415,12 +414,19 @@ public:
lfsw_n_cache.localFunctionSpace(),pn,lfsv_skeleton_n,
rl_view,rn_view);
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersection().indexInInside();
double weight = 1.;
if (ig.intersection().conforming()){
auto sub_entity_s_id = ig.intersection().indexInInside();
auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
weight = sub_entity_s.geometry().volume()/ig.geometry().volume();
}
// Load residual view into the local residual vector
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));
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << " , " << rn(lfsv_skeleton_n,i) << std::endl;
r_vec[i+offset] += weight*rl(lfsv_skeleton,i);
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
......@@ -449,6 +455,18 @@ public:
// position of quadrature point in local coordinates of elements
const auto p_local_s = ig.geometryInInside().global(it.position());
using LocalCoordinate = typename IG::LocalGeometry::LocalCoordinate;
LocalCoordinate sub_p_local_s, sub_p_local_n;
if (ig.intersection().conforming()) {
sub_p_local_s = it.position();
sub_p_local_n = it.position();
} else {
auto sub_entity_s_id = ig.intersection().indexInInside();
auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
auto sub_p_local_s = sub_entity_s.geometry().local(gtface.global(it.position()));
}
Dune::FieldVector<SolutionElement,IG::Entity::mydimension> y;
auto n_F = ig.unitOuterNormal(it.position());
auto BTransposed = ig.inside().geometry().jacobianTransposed(p_local_s);
......@@ -460,7 +478,7 @@ public:
// evaluate basis functions
lfsu_s.finiteElement().localBasis().evaluateFunction(p_local_s,phiu_s);
lfsv_skeleton.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
lfsv_skeleton.finiteElement().localBasis().evaluateFunction(sub_p_local_s,phiv_s);
// integration factor
const auto factor = it.weight() * ig.geometry().integrationElement(it.position());
......@@ -471,12 +489,10 @@ public:
}
}
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>
......@@ -493,7 +509,7 @@ public:
std::cout << "assembleUVBoundary" << std::endl;
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersectionIndex();
unsigned int offset = lfsv_volume.size() + lfsv_skeleton.size()*ig.intersection().indexInInside();
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
std::cout << i+offset << ": " << rl(lfsv_skeleton,i) << std::endl;
......
......@@ -27,6 +27,7 @@ class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::simpl
public:
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
for (unsigned int i = 0; i < size(); ++i)
li[i] = LocalKey(0,0,i);
......@@ -53,6 +54,7 @@ class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::cube>
public:
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
for (unsigned int i = 0; i < size(); ++i)
li[i] = LocalKey(0,0,i);
......
......@@ -298,7 +298,7 @@ public:
y /= detB;
BTransposed.mtv(y,gradphiv[i]);
}
dg_sign = 0.;
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
gradphiv = gradphiu;
......@@ -409,11 +409,29 @@ public:
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
// 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;
if (ig.intersection().conforming()) {
sub_p_local_s = it.position();
sub_p_local_n = it.position();
} else {
auto sub_entity_s_id = ig.intersection().indexInInside();
auto sub_entity_s = ig.inside().template subEntity<1>(sub_entity_s_id);
sub_p_local_s = sub_entity_s.geometry().local(gtface.global(it.position()));
auto sub_entity_n_id = ig.intersection().indexInOutside();
auto sub_entity_n = ig.outside().template subEntity<1>(sub_entity_n_id);
sub_p_local_n = sub_entity_n.geometry().local(gtface.global(it.position()));
}
// (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);
lfsv_s.finiteElement().localBasis().evaluateFunction(sub_p_local_s,phiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(sub_p_local_n,phiv_n);
// discrete gradient sign.
dg_sign = 0.;
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
......@@ -633,7 +651,7 @@ public:
if constexpr (dimDomainLocalTest == dim-1) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
dg_sign = 0.;
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
......
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