Commit 6d7303aa authored by Santiago Ospina's avatar Santiago Ospina

Construction of the raviart thomas working (not solving anything yet)

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent b253adf6
......@@ -142,18 +142,18 @@ int main(int argc, char** argv)
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);
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);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -168,18 +168,18 @@ int main(int argc, char** argv)
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 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!");
}
......@@ -187,21 +187,21 @@ int main(int argc, char** argv)
else{ // no adaptivity
auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile,helper);
switch(FEorder){
case 1:{
Sim<Cube<2,1>> 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;
}
// case 1:{
// Sim<Cube<2,1>> 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!");
}
......@@ -216,24 +216,24 @@ int main(int argc, char** argv)
if (gtype == "gmsh"){
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<3>>(inifile,helper);
switch(FEorder){
case 1:{
Sim<Simplex<3,1>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<Simplex<3,2>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<Simplex<3,3>> sim(helper,grid,inifile);
sim.set_policy(adapt_policy);
sim.run();
break;
}
// case 1:{
// Sim<Simplex<3,1>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 2:{
// Sim<Simplex<3,2>> sim(helper,grid,inifile);
// sim.set_policy(adapt_policy);
// sim.run();
// break;
// }
// case 3:{
// Sim<Simplex<3,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!");
}
......@@ -242,24 +242,24 @@ 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 3:{
Sim<CubeAdaptive<3,3>> 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);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......@@ -267,21 +267,21 @@ 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 3:{
Sim<Cube<3,3>> 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();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (grid.FEorder) not supported!");
}
......
#include "richards_simulation.hh"
#include "../solver/util_flux_reconstruction.hh"
#include "../solver/flux_reconstruction/rt_projection.hh"
namespace Dune{
namespace Dorie{
......@@ -130,6 +128,8 @@ void RichardsSimulation<Traits>::operator_setup()
if (helper.rank()==0)
std::cout << " Total number of DOF: " << DOFnumber << std::endl;
}
waterfrgf = std::make_shared<GFWaterFluxReconstruction>(*go0); //FIXME!
}
template<typename Traits>
......@@ -211,6 +211,7 @@ void RichardsSimulation<Traits>::update_adapters () const
condgf = std::make_shared<GFConductivity>(gv, fparam);
waterdgf = std::make_shared<GFWaterContent>(udgf, gv, fparam);
satdgf = std::make_shared<GFSaturation>(udgf, gv, fparam);
waterfrgf->update(*u);
}
template<typename Traits>
......@@ -220,66 +221,19 @@ void RichardsSimulation<Traits>::write_data () const
{
update_adapters();
constexpr int l = 0; // FIXME! this can be 0 or 1
// RT projection code
#if 1
RaviartThomasProjection<GO0,Traits::GridGeometryType,Traits::fem_order - l> proj(*go0);
#endif
// global lifting code
#if 0
constexpr bool frimp = FluxReconstructionImplemented<typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>::value;
#else
constexpr bool frimp = false;
#endif
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->template addVertexData<GFMatricHead>(udgf,"head");
vtkwriter->template addVertexData<GFWaterFlux>(fluxdgf,"flux");
vtkwriter->template addVertexData<GFConductivity>(condgf,"K_0");
vtkwriter->template addVertexData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addVertexData<GFSaturation>(satdgf,"Theta");
if constexpr (frimp)
{
using GFFluxReconstruction = Dune::Dorie::FluxReconstruction<
typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>;
std::shared_ptr<GFFluxReconstruction> frgf = std::make_shared<GFFluxReconstruction>(udgf);
vtkwriter->template addVertexData<GFFluxReconstruction>(frgf,"flux (reconstruction)");
} else {
std::cout << "WARNING: Flux reconstruction not available for this simulation!" << std::endl;
}
vtkwriter->template addVertexData<GFWaterFlux>(fluxdgf,"flux");
} else {
vtkwriter->template addCellData<GFMatricHead>(udgf,"head");
vtkwriter->template addCellData<GFWaterFlux>(fluxdgf,"flux");
vtkwriter->template addCellData<GFConductivity>(condgf,"K_0");
vtkwriter->template addCellData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addCellData<GFSaturation>(satdgf,"Theta");
if constexpr (frimp)
{
using GFFluxReconstruction = Dune::Dorie::FluxReconstruction<
typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>;
std::shared_ptr<GFFluxReconstruction> frgf = std::make_shared<GFFluxReconstruction>(udgf);
vtkwriter->template addCellData<GFFluxReconstruction>(frgf,"flux (reconstruction)");
} else {
std::cout << "WARNING: Flux reconstruction not available for this simulation!" << std::endl;
}
}
......
......@@ -30,6 +30,7 @@
#include "../solver/adapters/saturation.hh"
#include "../solver/adapters/water_flux.hh"
#include "../solver/adapters/conductivity.hh"
#include "../solver/flux_reconstruction/rt_projection.hh"
namespace Dune{
......@@ -50,6 +51,7 @@ struct RichardsSimulationTraits : public BaseTraits
{
static constexpr int dim = BaseTraits::dim;
static constexpr int fem_order = order;
static constexpr int flux_order = order - 1; //FIXME!
using RF = typename BaseTraits::RF;
using Grid = typename BaseTraits::Grid;
using GV = typename BaseTraits::GV;
......@@ -125,7 +127,8 @@ struct RichardsSimulationTraits : public BaseTraits
using GFWaterContent = Dune::Dorie::WaterContentAdapter<BaseTraits,FlowParameters,GFMatricHead>;
/// Saturation
using GFSaturation = Dune::Dorie::SaturationAdapter<BaseTraits,FlowParameters,GFMatricHead>;
/// Water Flux reconstruction
using GFWaterFluxReconstruction = Dune::Dorie::RaviartThomasFluxReconstruction<GO0,BaseTraits::GridGeometryType,flux_order>;
// -- Utility Class Definitions -- //
/// Custom VTK output writer
......@@ -229,6 +232,8 @@ private:
using GFWaterContent = typename Traits::GFWaterContent;
/// Saturation
using GFSaturation = typename Traits::GFSaturation;
/// Water Flux reconstruction
using GFWaterFluxReconstruction = typename Traits::GFWaterFluxReconstruction;
// -- Utility Class Definitions -- //
......@@ -273,12 +278,12 @@ protected:
std::shared_ptr<U> u;
mutable std::shared_ptr<GFMatricHead> udgf;
mutable std::shared_ptr<GFWaterFlux> fluxdgf;
mutable std::shared_ptr<GFConductivity> condgf;
mutable std::shared_ptr<GFWaterContent> waterdgf;
mutable std::shared_ptr<GFSaturation> satdgf;
mutable std::shared_ptr<GFMatricHead> udgf;
mutable std::shared_ptr<GFWaterFlux> fluxdgf;
mutable std::shared_ptr<GFConductivity> condgf;
mutable std::shared_ptr<GFWaterContent> waterdgf;
mutable std::shared_ptr<GFSaturation> satdgf;
mutable std::shared_ptr<GFWaterFluxReconstruction> waterfrgf;
std::unique_ptr<Writer> vtkwriter;
std::unique_ptr<AdaptivityHandler> adaptivity;
......
......@@ -9,6 +9,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/power.hh>
#include <dune/common/typetraits.hh>
#include <vector>
......@@ -22,7 +23,7 @@ namespace Dorie{
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalBasis
{
static_assert(true, "Volume Raviart Thomas basis function is not implemented!");
static_assert(Dune::AlwaysFalse<DF>::value, "Volume Raviart Thomas basis function is not implemented!");
};
template<class DF, class RF>
......
......@@ -16,7 +16,7 @@ namespace Dorie{
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalCoefficients
{
static_assert(true, "Volume Raviart Thomas coefficients is not implemented!");
static_assert(Dune::AlwaysFalse<DF>::value, "Volume Raviart Thomas coefficients is not implemented!");
};
template<class DF, class RF, int k, int dim>
......
......@@ -20,9 +20,6 @@ public:
VolumeRaviartThomasLocalCoefficients<DF,RF,k,d,gt>,
VolumeRaviartThomasLocalInterpolation>;
VolumeRaviartThomasLocalFiniteElement(const GeometryType& type)
{}
const typename Traits::LocalBasisType& localBasis() const
{
return _basis;
......
......@@ -2,12 +2,10 @@
#define DUNE_DORIE_RAVIART_THOMAS_PROJECTION_HH
#include "localengeine.hh"
#include "raviart_thomas/volumeraviartthomas.hh"
#include "skeleton_fem.hh"
#include "volume_fem.hh"
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/localfunctions/lagrange/qk.hh>
#include <dune/localfunctions/raviartthomas.hh>
#include <dune/pdelab/finiteelementmap/raviartthomasfem.hh>
#include <vector>
......@@ -15,27 +13,6 @@
namespace Dune{
namespace Dorie{
namespace {
template<class DF, class RF, int dim, int k, Dune::GeometryType::BasicType gt>
struct FaceFEHelper
{
static_assert(true,"Not implemented!");
};
template<class DF, class RF, int dim, int k>
struct FaceFEHelper<DF,RF,dim,k,Dune::GeometryType::simplex>
{
using type = Dune::PkLocalFiniteElement<DF,RF,dim-1,k>;
};
template<class DF, class RF, int dim, int k>
struct FaceFEHelper<DF,RF,dim,k,Dune::GeometryType::cube>
{
using type = Dune::QkLocalFiniteElement<DF,RF,dim-1,k>;
};
}
/*-------------------------------------------------------------------------*//**
* @brief Provide RT flux vector field from DG discretization. This code
* was adapted from Marian's code
......@@ -45,17 +22,13 @@ struct FaceFEHelper<DF,RF,dim,k,Dune::GeometryType::cube>
* @tparam order Order of the RT reconstruction.
*/
template<class GO, Dune::GeometryType::BasicType gt, int order>
class RaviartThomasProjection // : public GridFunctionBase<>
class RaviartThomasFluxReconstruction // : public GridFunctionBase<>
{
// ============
// types
// ============
using GFSU = typename GO::Traits::TrialGridFunctionSpace;
using GFSW = typename GO::Traits::TrialGridFunctionSpace;
using Domain = typename GO::Traits::Domain;
// grid view related stuff
using GV = typename GFSU::Traits::GridViewType;
using GV = typename GFSW::Traits::GridViewType;
using IndexSet = typename GV::IndexSet;
using DF = typename GO::Traits::DomainField;
using RF = typename GO::Traits::RangeField;
......@@ -63,81 +36,26 @@ class RaviartThomasProjection // : public GridFunctionBase<>
enum { dim = GV::dimension };
using FEMV = Dune::PDELab::RaviartThomasLocalFiniteElementMap<GV,DF,RF,order,gt>;
using ES = typename GFSU::Traits::EntitySet;
using GFSV = Dune::PDELab::GridFunctionSpace<ES,FEMV>;
using FEMU = Dune::PDELab::RaviartThomasLocalFiniteElementMap<GV,DF,RF,order,gt>;
using ES = typename GFSW::Traits::EntitySet;
using GFSU = Dune::PDELab::GridFunctionSpace<ES,FEMU>;
using LOP = typename GO::Traits::LocalAssembler::LocalOperator;
using MBE = typename GO::Traits::MatrixBackend;
using GOP = Dune::PDELab::GridOperator<GFSV,GFSU,LOP,MBE,DF,RF,JF>;
using GOP = Dune::PDELab::GridOperator<GFSU,GFSW,LOP,MBE,DF,RF,JF>;
using Range = typename GOP::Traits::Domain;
using LA = typename GOP::Traits::LocalAssembler;
using VolumeFE = Dune::Dorie::VolumeRaviartThomasLocalFiniteElement<DF,RF,order,dim,gt>;
using FaceFE = typename FaceFEHelper<DF,RF,dim,order,gt>::type;
using RaviartThomasEngine = Dune::Dorie::LocalRaviartThomasAssemblerEngine<LA,VolumeFE,FaceFE>;
// Basis functions to be used for interpolation
// using RTFactory = Dune::Dorie::RTBasisFactory<DF,RF,order,dim,gt>;
// using LFE = typename RTFactory::LFE; // the RT space
// using LFEFACE = typename RTFactory::TESTLFEFACE; // for face moments
// using BasisValueTypeRT = typename LFE::Traits::LocalBasisType::Traits::RangeType;
// using BasisValueTypeTestFace = typename LFEFACE::Traits::LocalBasisType::Traits::RangeType;
// using BasisValueTypeTestElem = typename RTPsi<DF,RF,order,dim,gt>::RangeType;
// enum { faces = RTFactory::faces };
// enum { facesize = RTFactory::facesize };
// enum { elementsize = RTFactory::elementsize };
// enum { totalsize = RTFactory::size };
// // Setup the grid operator projection
// using ES = typename GFSU::Traits::EntitySet;
// using GFSV = Dune::PDELab::GridFunctionSpace<ES,LFE>;
// using LOP = typename GO::Traits::LocalAssembler::LocalOperator;
// using MBE = typename GO::Traits::MatrixBackend;
// using GOP = Dune::PDELab::GridOperator<GFSU,GFSV,LOP,MBE,DF,RF,JF>;
// using Range = typename GOP::Traits::Range;
// using GFSRT = Dune::PDELab::GridFunctionSpace<ES,LFE>;
// using GFSP0 = Dune::PDELab::GridFunctionSpace<ES,LFE>;
// using GFSP1 = Dune::PDELab::GridFunctionSpace<ES,LFEFACE>;
// using VBE = Dune::PDELab::ISTL::VectorBackend<>;
// using OrderingTag = Dune::PDELab::DefaultLeafOrderingTag;
// using GFSP = Dune::PDELab::CompositeGridFunctionSpace<VBE,OrderingTag,GFSP0,GFSP1>;
// // iterators and entities
// using ElementIterator = typename GV::Traits::template Codim<0>::template Partition<Dune::Interior_Partition>::Iterator;
// using Element = typename GV::Traits::template Codim<0>::Entity;
// using JacobianInverseTransposed = typename Element::Geometry::JacobianInverseTransposed;
// using JacobianTransposed = typename Element::Geometry::JacobianTransposed;
// using IntersectionIterator = typename GV::IntersectionIterator;
// using Intersection = typename IntersectionIterator::Intersection;
// // DG local basis and discretization
// using LBTraits = typename GFSU::Traits::FiniteElementType::Traits::LocalBasisType::Traits;
// // TODO: ask Marian implications of this affine transformation
// static constexpr bool affine_transformation = (gt == Dune::GeometryType::simplex);
// // for the local system and dof storage
// typedef Dune::FieldVector<RF,totalsize> LocalDOF;
// typedef Dune::FieldMatrix<DF,totalsize,totalsize> Matrix;
// // RT basis function cache
// typedef typename LFE::Traits::LocalBasisType RTLocalBasisType;
// typedef Dune::PDELab::LocalBasisCache<RTLocalBasisType> RTCache;
using VolumeFEM = Dune::Dorie::VolumeRaviartThomasLocalFiniteElementMap<DF,RF,order,dim,gt>;
using GFSVVolume = Dune::PDELab::GridFunctionSpace<ES,VolumeFEM>;
using SkeletonFEM = typename Dune::Dorie::SkeletonFEM<DF,RF,dim,order,gt>::type;
using GFSVSkeleton = Dune::PDELab::GridFunctionSpace<ES,SkeletonFEM>;
using LocalRaviartThomasEngine = Dune::Dorie::LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton>;
// ============
// data members
// ============
const GO& _go;
GO& _go;
// grid related stuff
GV gv;
......@@ -166,7 +84,7 @@ class RaviartThomasProjection // : public GridFunctionBase<>
public:
RaviartThomasProjection(const GO& go, int intorderadd_=2)
RaviartThomasFluxReconstruction(GO& go, int intorderadd_=2)
: _go(go)
, gv(go.trialGridFunctionSpace().gridView())
// , is(gv.indexSet())
......@@ -178,17 +96,36 @@ public:
// , lastindex(-1)
{}
void update (const Domain& x)
void update (const Domain& p)
{
// // types needed for evaluation of u
// using LFSU = Dune::PDELab::LocalFunctionSpace<GFSU>;
// using LFSUCache = Dune::PDELab::LFSIndexCache<LFSU>;
// using size_type = typename LFSU::Traits::SizeType;
// using RFU = typename LFSU::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType;
// using RangeTypeU = typename LFSU::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeType;
// using JacobianTypeU = typename LFSU::Traits::FiniteElementType::Traits::LocalBasisType::Traits::JacobianType;
const auto& entity_set = _go.trialGridFunctionSpace().entitySet();
VolumeFEM volume_fem;
GFSVVolume gfsv_volume(entity_set,volume_fem);
SkeletonFEM skeleton_fem;
GFSVSkeleton gfsv_skeleton(entity_set,skeleton_fem);
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);
const auto& global_assembler_gop = gop.assembler();
const auto& local_assembler_gop = gop.localAssembler();
LocalRaviartThomasEngine local_raviart_thomas_engine(local_assembler_gop,gfsv_volume,gfsv_skeleton);
local_raviart_thomas_engine.setPrescription(p);
// local_raviart_thomas_engine.setSolution(x);
global_assembler_gop.assemble(local_raviart_thomas_engine);
}
};
......
#ifndef DUNE_DORIE_FLUX_RECONSTRUCTION_SKELETON_FEM_HH
#define DUNE_DORIE_FLUX_RECONSTRUCTION_SKELETON_FEM_HH
#include <dune/localfunctions/lagrange/pk.hh>
#include <dune/localfunctions/lagrange/qk.hh>
#include <dune/geometry/type.hh>
#include <dune/common/power.hh>
namespace Dune{
namespace Dorie{
template<class DF, class RF, int dim, int k>
class SkeletonPkFiniteElementMap
: public Dune::PDELab::SimpleLocalFiniteElementMap<Dune::PkLocalFiniteElement<DF,RF,dim-1,k>,dim-1>
{
public: