More documentation on raviart thomas flux reconstruction

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent 7ea21f40
...@@ -27,10 +27,10 @@ namespace Dorie{ ...@@ -27,10 +27,10 @@ namespace Dorie{
* @ingroup FluxReconstruction * @ingroup FluxReconstruction
* @ingroup FiniteElement * @ingroup FiniteElement
* *
* @tparam GV { description } * @tparam GV Grid view type.
* @tparam DF { description } * @tparam DF Domain field type.
* @tparam RF { description } * @tparam RF Range field type.
* @tparam k { description } * @tparam k Order.
*/ */
template<class GV, class DF, class RF, std::size_t k> template<class GV, class DF, class RF, std::size_t k>
class RaviartThomasSimplexDGLocalFiniteElementMap class RaviartThomasSimplexDGLocalFiniteElementMap
...@@ -65,29 +65,75 @@ public: ...@@ -65,29 +65,75 @@ public:
, _fe(Dune::GeometryTypes::simplex(GV::dimension),k) , _fe(Dune::GeometryTypes::simplex(GV::dimension),k)
{} {}
#ifdef DOXYGEN /* Documentation of the code */
/*-----------------------------------------------------------------------*//**
* @brief Return finite element for the given entity.
*
* @details The return value is a reference to Traits::FiniteElementType.
* If there is a different finite element for two elements then
* this type must be polymorphic. In this specific case, any
* entity maps to the same finite element.
*
* @param[in] e The entity that maps to the finite element.
*
* @tparam EntityType The type of the entity. Usually an entity as
* defined by dune-grid.
*
* @return The finite element associated to the entity e.
*/
template<class EntityType>
const typename Traits::FiniteElementType&
find (const EntityType& e);
#endif /* Documentation of the code */
/*-----------------------------------------------------------------------*//**
* @brief A FiniteElementMap is fixedSize iif the size of the local
* functions space for each GeometryType is fixed.
*
* @return A fixed size bool.
*/
static constexpr bool fixedSize() static constexpr bool fixedSize()
{ {
return true; return true;
} }
/*-----------------------------------------------------------------------*//**
* @brief Return if FiniteElementMap has degrees of freedom for given
* codimension.
*
* @param[in] codim The codimension.
*
* @return True if has degreed of freedom for codimension codim, False
* otherwise.
*/
static constexpr bool hasDOFs(int codim) static constexpr bool hasDOFs(int codim)
{ {
return (codim == 0); return (codim == 0);
} }
/*-----------------------------------------------------------------------*//**
* @brief If the FiniteElementMap is fixedSize, the size methods computes
* the number of DOFs for given GeometryType.
*
* @param[in] geometry_type The geometry type.
*
* @return The size of the fem for fixed size map.
*/
std::size_t size(GeometryType geometry_type) const std::size_t size(GeometryType geometry_type) const
{ {
switch (geometry_type.dim()) return (geometry_type.dim() == d) ? _fe.size(): 0;
{
case dim:
return size_rt;
default:
return 0;
}
} }
/*-----------------------------------------------------------------------*//**
* @brief Computes an upper bound for the local number of DOFs.
*
* @details This upper bound is used to avoid reallocations in std::vectors
* used during the assembly.
*
* @return The maximum local size of all finite elements used by the map.
*/
std::size_t maxLocalSize() const std::size_t maxLocalSize() const
{ {
return _fe.size(); return _fe.size();
......
...@@ -15,14 +15,15 @@ namespace Dorie{ ...@@ -15,14 +15,15 @@ namespace Dorie{
* *
* @author Santiago Ospina De Los Ríos * @author Santiago Ospina De Los Ríos
* @date 2018 * @date 2018
* *
* @tparam GV Grid view. * @tparam GV Grid view.
* @tparam DF Domain field. * @tparam DF Domain field.
* @tparam RF Range field. * @tparam RF Range field.
* @tparam k Order. * @tparam k Order.
* @tparam gt Geometry type. * @tparam gt Geometry type.
* @tparam <unnamed> Never instantiate this argument! Implementation * @tparam <unnamed> Never instantiate this argument! Implementation
* detail: Curiously recurring template pattern. * detail: Curiously recurring template pattern is used
* to select the correct map.
*/ */
template<class GV, class DF, class RF, std::size_t k, template<class GV, class DF, class RF, std::size_t k,
GeometryType::BasicType gt, class = void> GeometryType::BasicType gt, class = void>
...@@ -56,14 +57,14 @@ struct RaviartThomasLocalFiniteElementMapBaseSelector< ...@@ -56,14 +57,14 @@ struct RaviartThomasLocalFiniteElementMapBaseSelector<
* @f$\mathcal{RT}_k^{(dg)}\f$. * @f$\mathcal{RT}_k^{(dg)}\f$.
* @details This class is equivalent to * @details This class is equivalent to
* Dune::PDELab::RaviartThomasLocalFiniteElementMap plus elements * Dune::PDELab::RaviartThomasLocalFiniteElementMap plus elements
* for simplices in 3D (see * for raviart thomas in dg in simplices (see
* Dune::Dorie::RaviartThomasSimplexDGLocalFiniteElementMap). A big * Dune::Dorie::RaviartThomasSimplexDGLocalFiniteElementMap). A big
* difference with the dune-pdelab map is that these extra elements * difference with the dune-pdelab map is that these extra elements
* are not oriented, which make them to be not conforming in normal * are not oriented, which make them to be not conforming in normal
* direction. In other words, they cannot be used for usual finite * direction. In other words, they cannot be used for conforming
* element approximations. Nevertheless they are still useful for * finite element approximations. Nevertheless they are still useful
* flux reconstruction since we prescribe the continuity by using * for flux reconstruction since we prescribe the continuity by
* the numerical fluxes solved by the dG/FV methods. * using the numerical fluxes solved by the dG/FV methods.
* *
* Available elements: <table> <tr> <th colspan="2">Order</th> * Available elements: <table> <tr> <th colspan="2">Order</th>
* <th>0</th> <th>1</th> <th>2</th> <th>3</th> </tr> <tr> <td * <th>0</th> <th>1</th> <th>2</th> <th>3</th> </tr> <tr> <td
......
...@@ -49,10 +49,10 @@ public: ...@@ -49,10 +49,10 @@ public:
/*-----------------------------------------------------------------------*//** /*-----------------------------------------------------------------------*//**
* @brief Return finite element for the given entity. * @brief Return finite element for the given entity.
* *
* @details The return value is a reference to Traits::LocalBasisType. If * @details The return value is a reference to Traits::FiniteElementType.
* there is a different local basis for two elements then this * If there is a different finite element for two elements then
* type must be polymorphic. In this specific case, any entity * this type must be polymorphic. In this specific case, any
* maps to the same finite element. * entity maps to the same finite element.
* *
* @param[in] e The entity that maps to the finite element. * @param[in] e The entity that maps to the finite element.
* *
......
...@@ -60,15 +60,14 @@ class SkeletonPkFiniteElementMap ...@@ -60,15 +60,14 @@ class SkeletonPkFiniteElementMap
public: public:
#ifdef DOXYGEN /* Documentation of the code */ #ifdef DOXYGEN /* Documentation of the code */
/*-----------------------------------------------------------------------*//** /*-----------------------------------------------------------------------*//**
* @brief Return finite element \f$\mathcal{P}_k\f$ for the given * @brief Return finite element \f$\mathcal{P}_k\f$ for the given
* intersection. * intersection.
* *
* @details The return value is a reference to Traits::LocalBasisType. If * @details The return value is a reference to Traits::FiniteElementType.
* there is a different local basis for two elements then this * If there is a different finite element for two elements then
* type must be polymorphic. In this specific case, any entity * this type must be polymorphic. In this specific case, any
* maps to the same finite element. * entity maps to the same finite element.
* *
* @param[in] e The entity or intersection that maps to the finite * @param[in] e The entity or intersection that maps to the finite
* element. In case of intersection, it maps to the * element. In case of intersection, it maps to the
...@@ -192,10 +191,10 @@ public: ...@@ -192,10 +191,10 @@ public:
* @brief Return finite element \f$\mathcal{Q}_k\f$ for the given * @brief Return finite element \f$\mathcal{Q}_k\f$ for the given
* intersection. * intersection.
* *
* @details The return value is a reference to Traits::LocalBasisType. If * @details The return value is a reference to Traits::FiniteElementType.
* there is a different local basis for two elements then this * If there is a different finite element for two elements then
* type must be polymorphic. In this specific case, any entity * this type must be polymorphic. In this specific case, any
* maps to the same finite element. * entity maps to the same finite element.
* *
* @param[in] e The entity or intersection that maps to the finite * @param[in] e The entity or intersection that maps to the finite
* element. In case of intersection, it maps to the * element. In case of intersection, it maps to the
......
...@@ -10,7 +10,51 @@ ...@@ -10,7 +10,51 @@
@defgroup FluxReconstruction Flux reconstruction @defgroup FluxReconstruction Flux reconstruction
@{ @{
@ingroup Common @ingroup Common
@todo document flux reconstruction
### Flux reconstruction with Raviart Thomas elements.
Following the notation of Di Petro 2012 the discrete problem for discontinuous
galerkin problems can be written as:
For all @f$T\in\mathcal{T}_h@f$ and for all @f$\xi\in\mathbb{P}^k_d(T)@f$,
@f[
\int_T G_h^l(u_h)\cdot\nabla\xi + \sum_{F\in\mathcal{F}_T}\epsilon_{T,F}
\int_F\phi_F(u_h)\xi=\int_Tf\xi \qquad (1)
@f]
with @f$l\in\{k-1,k\}@f$,@f$\epsilon_{T,F}=n_F\cdot n_F@f$, and the numerical
flux
@f[
\phi(u_h):=-\{\{\nabla_hu_h\}\}\cdot n_F+\frac{\eta}{h_F}[[u_h]]
@f]
Following that, one can solve the degrees of freedom of raviart thomas element
@f$ \sigma_h^*\in\mathbb{RTN}^l_d(\mathcal{T}),\,l\in\{k-1,k\}@f$ problem
solving
(i) For all @f$F\in\mathcal{F_h}@f$ and all @f$q\in\mathbb{P}^l_{d-1}(F)@f$,
@f[
\int_F(\sigma_h^*\cdot n_F)q:=\int_F\phi_F(u_h)q \qquad (2a)
@f]
(ii) For all @f$T\in\mathcal{T_h}@f$ and all
@f$r\in[\mathbb{P}^{l-1}_d(T)]^d@f$,
@f[
\int_T\sigma_h^*r:=-\int_TG_h^l(u_h)\cdot r \qquad (2b)
@f]
In the context of the flux reconstruction code
@f$r\in[\mathbb{P}^{l-1}_d(T)]^d@f$ is what we call the volume raviart thomas
elements (see Dune::Dorie::RaviartThomasVolumeLocalFiniteElementMap),
@f$q\in\mathbb{P}^l_{d-1}(F)@f$ the skeleton elements (see
Dune::Dorie::SkeletonFiniteElementMap), @f$u_h@f$ the prescribed function and
@f$ \sigma_h^*@f$ the ansatz function (see
Dune::Dorie::RaviartThomasLocalFiniteElementMap).
Since local operators usually solve the equation (1), it is not to hard to
modify them in order to accept non-conforming local functions spaces that
solve (2). Since integrals for skeleton and for volume integrals need
different test functions, it is implemented the
Dune::Dorie::LocalRaviartThomasAssemblerEngine which can forward the right
local function spaces the two different integrals to the local operator.
@} @}
**/ **/
...@@ -33,26 +33,34 @@ public: ...@@ -33,26 +33,34 @@ public:
}; };
/*-------------------------------------------------------------------------*//** /*-------------------------------------------------------------------------*//**
* @brief Provide RT flux vector field from DG discretization. This code * @brief Raviart Thomas flux reconstruction for discontinuous Galerkin
* was adapted from Marian's code * grid operators.
* * @details The local operator used in the grid operator must be able to
* handle non-conforming finite elements see @ref FluxReconstruction
* for details.
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @ingroup FluxReconstruction
*
*
* @tparam GO The GridOperator type. * @tparam GO The GridOperator type.
* @tparam gt Geometry type of the grid. * @tparam gt Geometry type of the grid.
* @tparam order Order of the RT reconstruction. * @tparam order Order of the Raviart Thomas reconstruction.
*/ */
template<class GO, Dune::GeometryType::BasicType gt, unsigned int order> template<class GO, Dune::GeometryType::BasicType gt, unsigned int order>
class RaviartThomasFluxReconstruction class RaviartThomasFluxReconstruction
: public Dune::PDELab::GridFunctionBase< : public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits< Dune::PDELab::GridFunctionTraits<
typename GO::Traits::TrialGridFunctionSpace::Traits::GridViewType, typename GO::Traits::TrialGridFunctionSpace::Traits::GridViewType,
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType, typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension, GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension,
Dune::FieldVector< Dune::FieldVector<
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType, typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension
> >, > >,
RaviartThomasFluxReconstruction<GO,gt,order> RaviartThomasFluxReconstruction<GO,gt,order>
> >
{ {
using GFSW = typename GO::Traits::TrialGridFunctionSpace; using GFSW = typename GO::Traits::TrialGridFunctionSpace;
using Domain = typename GO::Traits::Domain; using Domain = typename GO::Traits::Domain;
...@@ -94,13 +102,23 @@ public: ...@@ -94,13 +102,23 @@ public:
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension
> >; > >;
RaviartThomasFluxReconstruction(GO& go, int intorderadd_=2) /*---------------------------------------------------------------------*//**
* @brief Constructor of the class
*
* @warning This class keeps a reference to the grid operator.
*
* @param go Grid operator.
* @param[in] int_order_add Integration integer to add to the quadrature
* rule.
*/
RaviartThomasFluxReconstruction(GO& go, int int_order_add=2)
: _go(go) : _go(go)
, _gv(go.trialGridFunctionSpace().gridView()) , _gv(go.trialGridFunctionSpace().gridView())
, _femu(_gv) , _femu(_gv)
, _gfsu(_go.trialGridFunctionSpace().entitySet(),_femu) , _gfsu(_go.trialGridFunctionSpace().entitySet(),_femu)
, _x(_gfsu,0.0) , _x(_gfsu,0.0)
, _dgfp(_gfsu,_x) , _dgfp(_gfsu,_x)
, _int_order_add(int_order_add)
{} {}
void update (const Domain& p) void update (const Domain& p)
...@@ -129,7 +147,7 @@ public: ...@@ -129,7 +147,7 @@ public:
const auto& local_assembler_gop = gop.localAssembler(); const auto& local_assembler_gop = gop.localAssembler();
LocalRaviartThomasEngine local_raviart_thomas_engine(local_assembler_gop,gfsv_volume,gfsv_skeleton); LocalRaviartThomasEngine local_raviart_thomas_engine(local_assembler_gop,gfsv_volume,gfsv_skeleton,_int_order_add);
local_raviart_thomas_engine.setPrescription(p); local_raviart_thomas_engine.setPrescription(p);
local_raviart_thomas_engine.setSolution(_x); local_raviart_thomas_engine.setSolution(_x);
...@@ -161,6 +179,7 @@ private: ...@@ -161,6 +179,7 @@ private:
GFSU _gfsu; GFSU _gfsu;
Range _x; Range _x;
Dune::PDELab::DiscreteGridFunctionPiola<GFSU,Range> _dgfp; Dune::PDELab::DiscreteGridFunctionPiola<GFSU,Range> _dgfp;
const int _int_order_add;
}; };
} // namespace Dorie } // namespace Dorie
......
...@@ -14,16 +14,20 @@ namespace Dorie{ ...@@ -14,16 +14,20 @@ namespace Dorie{
/*-------------------------------------------------------------------------*//** /*-------------------------------------------------------------------------*//**
* @brief The assembler for standard DUNE grid. * @brief The assembler for standard DUNE grid.
* @details This is a modification of the Dune::PDELab::DefaultAssembler in * @details This is a modification of the Dune::PDELab::DefaultAssembler in
* PDELab. It binds and unbinds inside local function spaces with * PDELab. It binds and unbinds inside and outside local function
* the intersection. Notice that to merge this into PDELab, it is * spaces with the intersection. Notice that to merge this into
* necessary to remove the inside bind/unbind methods in all the * PDELab, it is necessary to remove the inside/outside bind/unbind
* local assamblers! They are not used currently, but if we add this * methods in all the local assamblers! They are not used currently,
* changes, local function spaces will be binded several times for * but if we add this changes, local function spaces will be binded
* the same entity. * several times for the same entity.
* *
* @author <a href="https://github.com/dune-project/dune-pdelab/blob/master/LICENSE.md">dune-pdelab team</a> * @author <a
* href="https://github.com/dune-project/dune-pdelab/blob/master/LICENSE.md">dune-pdelab
* team</a>
* @date 2011-2018 * @date 2011-2018
* @copyright <a href="https://github.com/dune-project/dune-pdelab/blob/master/LICENSE.md">dune-pdelab licence</a> * @copyright <a
* href="https://github.com/dune-project/dune-pdelab/blob/master/LICENSE.md">dune-pdelab
* licence</a>
* @ingroup FluxReconstruction * @ingroup FluxReconstruction
* *
* @tparam GFSU GridFunctionSpace for ansatz functions. * @tparam GFSU GridFunctionSpace for ansatz functions.
...@@ -35,14 +39,14 @@ template<typename GFSU, typename GFSV, typename CU, typename CV> ...@@ -35,14 +39,14 @@ template<typename GFSU, typename GFSV, typename CU, typename CV>
class DefaultAssembler { class DefaultAssembler {
public: public:
//! Types related to current grid view //! @name Types related to current grid view
//! @{ //! @{
using EntitySet = typename GFSU::Traits::EntitySet; using EntitySet = typename GFSU::Traits::EntitySet;
using Element = typename EntitySet::Element; using Element = typename EntitySet::Element;
using Intersection = typename EntitySet::Intersection; using Intersection = typename EntitySet::Intersection;
//! @} //! @}
//! Grid function spaces //! @name Grid function spaces
//! @{ //! @{
typedef GFSU TrialGridFunctionSpace; typedef GFSU TrialGridFunctionSpace;
typedef GFSV TestGridFunctionSpace; typedef GFSV TestGridFunctionSpace;
...@@ -54,6 +58,14 @@ public: ...@@ -54,6 +58,14 @@ public:
//! Static check on whether this is a Galerkin method //! Static check on whether this is a Galerkin method
static const bool isGalerkinMethod = std::is_same<GFSU,GFSV>::value; static const bool isGalerkinMethod = std::is_same<GFSU,GFSV>::value;
/*-----------------------------------------------------------------------*//**
* @brief Constructs the object.
*
* @param[in] gfsu_ The ansatz grid function space
* @param[in] gfsv_ The test grid function space
* @param[in] cu_ Constraints for ansatz functions.
* @param[in] cv_ Constraints for test functions.
*/
DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_, const CU& cu_, const CV& cv_) DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_, const CU& cu_, const CV& cv_)
: gfsu(gfsu_) : gfsu(gfsu_)
, gfsv(gfsv_) , gfsv(gfsv_)
...@@ -65,6 +77,12 @@ public: ...@@ -65,6 +77,12 @@ public:
, lfsvn(gfsv_) , lfsvn(gfsv_)
{ } { }
/*-----------------------------------------------------------------------*//**
* @brief Constructs the object.
*
* @param[in] gfsu_ The ansatz grid function space
* @param[in] gfsv_ The test grid function space
*/
DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_) DefaultAssembler (const GFSU& gfsu_, const GFSV& gfsv_)
: gfsu(gfsu_) : gfsu(gfsu_)
, gfsv(gfsv_) , gfsv(gfsv_)
...@@ -88,6 +106,14 @@ public: ...@@ -88,6 +106,14 @@ public:
return gfsv; return gfsv;
} }
/*-----------------------------------------------------------------------*//**
* @brief Assemble a local engeine that follows the structure of the
* Dune::PDELab::LocalAssemblerEngineBase.
*
* @param assembler_engine The local assembler engine
*
* @tparam LocalAssemblerEngine The local assembler engeine type
*/
template<class LocalAssemblerEngine> template<class LocalAssemblerEngine>
void assemble(LocalAssemblerEngine & assembler_engine) const void assemble(LocalAssemblerEngine & assembler_engine) const
{ {
......
...@@ -9,7 +9,7 @@ except NameError: ...@@ -9,7 +9,7 @@ except NameError:
pass pass
# paths set by cmake # paths set by cmake
DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Release/dorie" DORIEDIR = "/home/saospina/Codes/DUNE_INSTALL/Debug/dorie"
MPIEXEC = "/usr/bin/mpiexec" MPIEXEC = "/usr/bin/mpiexec"
MPIEXEC_NUMPROC_FLAG = "-n" MPIEXEC_NUMPROC_FLAG = "-n"
MPIEXEC_PREFLAG = "" MPIEXEC_PREFLAG = ""
......
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