Commit b253adf6 authored by Santiago Ospina's avatar Santiago Ospina

Added volume raviart thomas finite elements and an sketch of the local engeine

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 2f453f62
......@@ -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,8 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::cube>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<2>,Geo::simplex>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
......@@ -8,8 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,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<3>,Geo::cube>,2>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,2>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,2>>;
} // namespace Dorie
} // namespace Dune
......@@ -8,8 +8,8 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3>>;
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::cube>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<UGGrid<3>,Geo::simplex>,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<2>,Geo::cube>,1>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<2>,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<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
......@@ -8,7 +8,7 @@
namespace Dune{
namespace Dorie{
template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>;
// template class RichardsSimulation<RichardsSimulationTraits<BaseTraits<YaspGrid<3>,Geo::cube>,3>>;
} // namespace Dorie
} // namespace Dune
\ No newline at end of file
#include "richards_simulation.hh"
#include "../solver/util_flux_reconstruction.hh"
#include "../solver/flux_reconstruction/rt_projection.hh"
namespace Dune{
namespace Dorie{
......@@ -215,17 +216,29 @@ void RichardsSimulation<Traits>::update_adapters () const
template<typename Traits>
void RichardsSimulation<Traits>::write_data () const
{
constexpr int l = 0; // FIXME! this can be 0 or 1
if (vtkwriter)
{
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");
......
#ifndef DUNE_DORIE_RAVIART_THOMAS_RESIDUALENGINE_HH
#define DUNE_DORIE_RAVIART_THOMAS_RESIDUALENGINE_HH
#include <dune/pdelab/gridfunctionspace/localvector.hh>
#include <dune/pdelab/gridoperator/common/assemblerutilities.hh>
#include <dune/pdelab/gridoperator/common/localassemblerenginebase.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/localoperator/callswitch.hh>
namespace Dune{
namespace Dorie{
/**
\brief The local assembler engine for DUNE grids which
assembles the residual vector
\tparam LA The local assembler
*/
template<typename LA, typename VolumeFE, typename SkeletonFE>
class LocalRaviartThomasAssemblerEngine
: public Dune::PDELab::LocalAssemblerEngineBase
{
public:
template<typename TrialConstraintsContainer, typename TestConstraintsContainer>
bool needsConstraintsCaching(const TrialConstraintsContainer& cu, const TestConstraintsContainer& cv) const
{
return false;
}
//! The type of the wrapping local assembler
typedef LA LocalAssembler;
//! The type of the local operator
typedef typename LA::LocalOperator LOP;
//! The type of the residual vector
typedef typename LA::Traits::Residual Residual;
typedef typename Residual::ElementType ResidualElement;
//! The type of the solution vector
typedef typename LA::Traits::Solution Solution;
typedef typename Solution::ElementType SolutionElement;
//! The local function spaces
typedef typename LA::LFSU LFSU;
typedef typename LA::LFSUCache LFSUCache;
typedef typename LFSU::Traits::GridFunctionSpace GFSU;
typedef typename LA::LFSV LFSV;
typedef typename LA::LFSVCache LFSVCache;
typedef typename LFSV::Traits::GridFunctionSpace GFSV;
typedef typename Solution::template ConstLocalView<LFSUCache> SolutionView;
typedef typename Residual::template LocalView<LFSVCache> ResidualView;
/**
\brief Constructor
\param [in] local_assembler_ The local assembler object which
creates this engine
*/
LocalRaviartThomasAssemblerEngine(
const LocalAssembler & local_assembler_,
const VolumeFE & volume_fe_,
const SkeletonFE & skeleton_fe_
)
: local_assembler(local_assembler_),
pvolume_fe(stackobject_to_shared_ptr(volume_fe_)),
pskeleton_fe(stackobject_to_shared_ptr(skeleton_fe_)),
lop(local_assembler_.localOperator()),
rl_view(rl,1.0),
rn_view(rn,1.0)
{}
LocalRaviartThomasAssemblerEngine(
const LocalAssembler & local_assembler_,
const std::shared_ptr<const VolumeFE> & volume_fe_,
const std::shared_ptr<const SkeletonFE> & skeleton_fe_
)
: local_assembler(local_assembler_),
pvolume_fe(volume_fe_),
pskeleton_fe(skeleton_fe_),
lop(local_assembler_.localOperator()),
rl_view(rl,1.0),
rn_view(rn,1.0)
{}
//! Query methods for the global grid assembler
//! @{
bool requireSkeleton() const
{ return ( local_assembler.doAlphaSkeleton() || local_assembler.doLambdaSkeleton() ); }
bool requireSkeletonTwoSided() const
{ return local_assembler.doSkeletonTwoSided(); }
bool requireUVVolume() const
{ return local_assembler.doAlphaVolume(); }
bool requireVVolume() const
{ return local_assembler.doLambdaVolume(); }
bool requireUVSkeleton() const
{ return local_assembler.doAlphaSkeleton(); }
bool requireVSkeleton() const
{ return local_assembler.doLambdaSkeleton(); }
bool requireUVBoundary() const
{ return local_assembler.doAlphaBoundary(); }
bool requireVBoundary() const
{ return local_assembler.doLambdaBoundary(); }
bool requireUVVolumePostSkeleton() const
{ return local_assembler.doAlphaVolumePostSkeleton(); }
bool requireVVolumePostSkeleton() const
{ return local_assembler.doLambdaVolumePostSkeleton(); }
//! @}
//! Public access to the wrapping local assembler
const LocalAssembler & localAssembler() const
{
return local_assembler;
}
//! Trial space constraints
const typename LocalAssembler::Traits::TrialGridFunctionSpaceConstraints& trialConstraints() const
{
return localAssembler().trialConstraints();
}
//! Test space constraints
const typename LocalAssembler::Traits::TestGridFunctionSpaceConstraints& testConstraints() const
{
return localAssembler().testConstraints();
}
//! Set current residual vector. Should be called prior to
//! assembling.
void setResidual(Residual & residual_)
{
global_rl_view.attach(residual_);
global_rn_view.attach(residual_);
}
//! Set current solution vector. Should be called prior to
//! assembling.
void setSolution(const Solution & solution_)
{
global_sl_view.attach(solution_);
global_sn_view.attach(solution_);
}
//! 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 & lfsv_cache)
{
global_sl_view.bind(lfsu_cache);
xl.resize(lfsu_cache.size());
}
template<typename EG, typename LFSVC>
void onBindLFSV(const EG & eg, const LFSVC & lfsv_cache)
{
global_rl_view.bind(lfsv_cache);
rl.assign(lfsv_cache.size(),0.0);
}
template<typename IG, typename LFSUC, typename LFSVC>
void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSVC & lfsv_cache)
{
global_sl_view.bind(lfsu_cache);
xl.resize(lfsu_cache.size());
}
template<typename IG, typename LFSUC, typename LFSVC>
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)
{
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 & lfsv_cache)
{
global_rl_view.bind(lfsv_cache);
rl.assign(lfsv_cache.size(),0.0);
}
template<typename IG, typename LFSVC>
void onBindLFSVOutside(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache)
{
global_rn_view.bind(lfsv_n_cache);
rn.assign(lfsv_n_cache.size(),0.0);
}
//! @}
//! 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 & lfsv_cache)
{
global_rl_view.add(rl);
global_rl_view.commit();
}
template<typename IG, typename LFSVC>
void onUnbindLFSVInside(const IG & ig, const LFSVC & lfsv_cache)
{
global_rl_view.add(rl);
global_rl_view.commit();
}
template<typename IG, typename LFSVC>
void onUnbindLFSVOutside(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache)
{
global_rn_view.add(rn);
global_rn_view.commit();
}
//! @}
//! Methods for loading of the local function's coefficients
//! @{
template<typename LFSUC>
void loadCoefficientsLFSUInside(const LFSUC & lfsu_s_cache)
{
global_sl_view.read(xl);
}
template<typename LFSUC>
void loadCoefficientsLFSUOutside(const LFSUC & lfsu_n_cache)
{
global_sn_view.read(xn);
}
template<typename LFSUC>
void loadCoefficientsLFSUCoupling(const LFSUC & lfsu_c_cache)
{
DUNE_THROW(Dune::NotImplemented,"No coupling lfsu available for ");
}
//! @}
//! Notifier functions, called immediately before and after assembling
//! @{
void postAssembly(const GFSU& gfsu, const GFSV& gfsv)
{
if(local_assembler.doPostProcessing())
Dune::PDELab::constrain_residual(local_assembler.testConstraints(),global_rl_view.container());
}
//! @}
//! Assembling methods
//! @{
/** Assemble on a given cell without function spaces.
\return If true, the assembling for this cell is assumed to
be complete and the assembler continues with the next grid
cell.
*/
template<typename EG>
bool assembleCell(const EG & eg)
{
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 & lfsv_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolume>::
// alpha_volume(lop,eg,lfsu_cache.localFunctionSpace(),xl,lfsv_cache.localFunctionSpace(),rl_view);
}
template<typename EG, typename LFSVC>
void assembleVVolume(const EG & eg, const LFSVC & lfsv_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaVolume>::
// lambda_volume(lop,eg,lfsv_cache.localFunctionSpace(),rl_view);
}
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)
{
// 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);
}
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 LFSUC, typename LFSVC>
void assembleUVBoundary(const IG & ig, const LFSUC & lfsu_s_cache, const LFSVC & lfsv_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);
}
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 LFSUC, typename LFSVC>
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)
{
DUNE_THROW(Dune::NotImplemented,"Assembling of coupling spaces is not implemented for ");
}
template<typename IG, typename LFSVC>
static void assembleVEnrichedCoupling(const IG & ig,
const LFSVC & lfsv_s_cache,
const LFSVC & lfsv_n_cache,
const LFSVC & lfsv_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 & lfsv_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
// alpha_volume_post_skeleton(lop,eg,lfsu_cache.localFunctionSpace(),xl,lfsv_cache.localFunctionSpace(),rl_view);
}
template<typename EG, typename LFSVC>
void assembleVVolumePostSkeleton(const EG & eg, const LFSVC & lfsv_cache)
{
// rl_view.setWeight(local_assembler.weight());
// Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doLambdaVolumePostSkeleton>::
// lambda_volume_post_skeleton(lop,eg,lfsv_cache.localFunctionSpace(),rl_view);
}
//! @}
private:
//! Reference to the wrapping local assembler object which
//! constructed this engine
const LocalAssembler & local_assembler;
//! Reference to the local operator
const LOP & lop;
//! Pointer to the current residual vector in which to assemble
ResidualView global_rl_view;
ResidualView global_rn_view;
//! Pointer to the current residual vector in which to assemble
SolutionView global_sl_view;
SolutionView global_sn_view;
//! The local vectors and matrices as required for assembling
//! @{
typedef Dune::PDELab::TrialSpaceTag LocalTrialSpaceTag;
typedef Dune::PDELab::TestSpaceTag LocalTestSpaceTag;
typedef Dune::PDELab::LocalVector<SolutionElement, LocalTrialSpaceTag> SolutionVector;
typedef Dune::PDELab::LocalVector<ResidualElement, LocalTestSpaceTag> ResidualVector;
//! Inside local coefficients
SolutionVector xl;
//! Outside local coefficients
SolutionVector xn;
//! Inside local residual
ResidualVector rl;
//! Outside local residual
ResidualVector rn;
//! Inside local residual weighted view
typename ResidualVector::WeightedAccumulationView rl_view;
//! Outside local residual weighted view
typename ResidualVector::WeightedAccumulationView rn_view;
//! @}
const std::shared_ptr<const VolumeFE> pvolume_fe;
const std::shared_ptr<const SkeletonFE> pskeleton_fe;
}; // End of class LocalRaviartThomasAssemblerEngine
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_RAVIART_THOMAS_RESIDUALENGINE_HH
#ifndef DUNE_DORIE_VOLUME_RAVIART_THOMAS_COEFFICIENTS_HH
#define DUNE_DORIE_VOLUME_RAVIART_THOMAS_COEFFICIENTS_HH
#include <dune/localfunctions/common/localkey.hh>
#include <dune/geometry/type.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/power.hh>
#include <vector>
namespace Dune{
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!");
};
template<class DF, class RF, int k, int dim>
class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::simplex>
{
static constexpr int volsize = (dim == 2) ? dim*k*(k+1)/2
: dim*k*(k+1)*(k+2)/6;
public: