Document the flux reconstruction code.

Signed-off-by: Santiago Ospina De Los Ríos's avatarSantiago Ospina <santiago.ospina@iup.uni-heidelberg.de>
parent d803ff0d
......@@ -18,5 +18,8 @@ USE_MATHJAX = YES
# subdirectory from a directory tree whose root is specified with the INPUT tag.
# exclude versions of the operators that are not used
# EXCLUDE += @top_srcdir@/src/operator_DG_old.hh \
# @top_srcdir@/src/operator_DG_no_upwind.hh
# EXCLUDE += @top_srcdir@/...
# Include documents with other extesions
FILE_PATTERNS += *.dox
add_subdirectory("flux_reconstruction")
#install headers
install(FILES
boundary_condition.hh
......
// These groups should gather information of code that belong to fairly similar
// categories. In them we aim to explain the code: why is done in that specific
// way, how carful someone has to be to modify it, specific name conventions, etc.
// These descriptions are not for regular users but for future developers (maybe
// the same developer in some months/years later).
/**
@defgroup Common Common
@{
@todo document models
@}
@defgroup FluxReconstruction Flux reconstruction
@{
@ingroup Common
@todo document the flux reconstruction
@}
**/
add_subdirectory("raviart_thomas")
#install headers
install(FILES
assembler.hh
minimal_local_function_space.hh
raviart_thomas.hh
volume_fem.hh
local_engeine.hh
skeleton_fem.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/common/flux_reconstruction)
\ No newline at end of file
This diff is collapsed.
#ifndef DUNE_DORIE_MINIMAL_LOCAL_FUNCTION_SPACE_HH
#define DUNE_DORIE_MINIMAL_LOCAL_FUNCTION_SPACE_HH
#include <dune/localfunctions/common/interfaceswitch.hh>
#include <dune/pdelab/gridfunctionspace/localfunctionspace.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Class for minimal local function space.
* @details This is a special local function space for the raviart thomas
* engine. In particular, it is thought to work with one grid
* function space, that is, without composite or vector spaces so
* that we avoid anything of the type tree and ordering nightmare.
* It is simply not needed here. Indeed we don't even need the dof
* indices since the engine does a local solution and there is no
* need to map to entities. In any case, the reason of this object
* is because the user is expecting a local function space as an
* argument of the local operator (even though no many really know
* how a local function space looks like ¬_¬).
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
*
* @tparam GFS The GridFunctionSpace.
* @tparam TAG The dune-pdelab space tag.
*/
template <class GFS, class TAG=Dune::PDELab::AnySpaceTag>
class MinimalLocalFunctionSpace
{
static_assert(not GFS::Traits::isComposite,
"This local function space does not work with "
"composited or power grid function spaces");
//! The dof index type to store indices of each dof.
using DOFIndex = typename Dune::PDELab::gfs_to_lfs<GFS>::DOFIndex;
public:
//! Traits following the Dune::PDELab::LocalFunctionSpace structure.
using Traits = typename Dune::PDELab::LocalFunctionSpace<GFS>::Traits;
//! Finite element switch following the Dune::FiniteElementInterfaceSwitch
//! structure defined by dune-localfunctions.
using FESwitch = Dune::FiniteElementInterfaceSwitch<
typename Traits::FiniteElementType>;
/*-----------------------------------------------------------------------*//**
* @brief Constructs the local function space.
*
* @param[in] gfs The pointer to a grid function space.
*/
MinimalLocalFunctionSpace(std::shared_ptr<const GFS> gfs)
: pgfs(gfs)
{}
/*-----------------------------------------------------------------------*//**
* @brief Constructs the local function space.
*
* @param[in] gfs The grid function space.
*/
MinimalLocalFunctionSpace(const GFS & gfs)
: pgfs(stackobject_to_shared_ptr(gfs))
{}
/**
* @brief Copy constructs the local function space.
*
* @param[in] lfs The local function space.
*/
MinimalLocalFunctionSpace(const MinimalLocalFunctionSpace & lfs) {}
/*-----------------------------------------------------------------------*//**
* @brief Maps an index of the shape function to the index into the local
* function space.
* @details In composite and power grid function spaces is here where the
* magic of the typetree happens.
*
* @param[in] index The index of the shape function.
*
* @return Index representing a degree of freedom of the local function
* space.
*/
typename Traits::IndexContainer::size_type
localIndex(typename Traits::IndexContainer::size_type index) const
{
return index;
}
/*-----------------------------------------------------------------------*//**
* @brief Total size of the local function space.
*
* @return The size.
*/
typename Traits::IndexContainer::size_type size () const
{
return n;
}
/**
* @brief Maps given index in this local function space to its
* corresponding global MultiIndex.
*
* @param[in] index The local index value from the range 0,...,size()-1
*
* @return A const reference to the associated, globally unique
* MultiIndex. Note that the returned object may (and must) be
* copied if it needs to be stored beyond the time of the next
* call to bind() on this LocalFunctionSpace (e.g. when the
* MultiIndex is used as a DOF identifier in a constraints
* container).
*/
const typename Traits::DOFIndex&
dofIndex(typename Traits::IndexContainer::size_type index) const
{
DUNE_THROW(Dune::NotImplemented,
"This local function space never constructs dof indices");
return Traits::DOFIndex();
}
/*-----------------------------------------------------------------------*//**
* @brief Binds a given entity or intersection to the respective finite
* element corresponding to the finite element map.
* @details This is indeed instantiated as an intersection here. Since we
* don't care about the dof intex, a template is is ok to forward
* it to the finite element map. In other cases one whould have to
* hard code the type.
*
* @param e The entity or intersection.
*
* @tparam Entity The entity or intersection type.
*/
template<class Entity>
void bind(Entity& e)
{
FESwitch::setStore(pfe, pgfs->finiteElementMap().find(e));
n = FESwitch::basis(*pfe).size();
}
/*-----------------------------------------------------------------------*//**
* @brief Get finite element
*
* @return The finite element.
*/
const typename Traits::FiniteElementType& finiteElement () const
{
assert(pfe);
return *pfe;
}
private:
std::shared_ptr<const GFS> pgfs;
typename Traits::IndexContainer::size_type n;
typename FESwitch::Store pfe;
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_MINIMAL_LOCAL_FUNCTION_SPACE_HH
#ifndef DUNE_DORIE_RAVIART_THOMAS_PROJECTION_HH
#define DUNE_DORIE_RAVIART_THOMAS_PROJECTION_HH
#include "localengeine.hh"
#include "skeleton_fem.hh"
#include "volume_fem.hh"
#include "raviart_thomas_operator/assembler.hh"
#include <dune/pdelab/finiteelementmap/raviartthomasfem.hh>
#include <vector>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/skeleton_fem.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/volume_fem.hh>
#include <dune/dorie/common/flux_reconstruction/assembler.hh>
#include <dune/dorie/common/flux_reconstruction/local_engeine.hh>
namespace Dune{
namespace Dorie{
......
add_subdirectory("volume")
#install headers
install(FILES
volume_raviart_thomas.hh
volume_fem.hh
skeleton_fem.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/common/flux_reconstruction/raviart_thomas)
\ No newline at end of file
#install headers
install(FILES
volumertlocalbasis.hh
volumertlocalcoefficients.hh
volumertlocalinterpolation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/dorie/common/flux_reconstruction/raviartthomas/volume)
\ No newline at end of file
......@@ -16,16 +16,112 @@
namespace Dune{
namespace Dorie{
//======================================================================
// The Volume Raviart Thomas basis functions
//======================================================================
/*-------------------------------------------------------------------------*//**
* @brief Class for local basis of the volume raviart thomas finite
* element.
* @details This class is used to reconstruct fluxes from local operators
* that solve fluxes in the residual vector. In particular, this
* volume raviart thomas is equivalent to the usual raviart thomas
* element but excluding the degrees of freedom located on the faces
* (codim = 1).
*
* @author Santiago Ospina De Los Ríos implemented it in DORiE out of source
* files sent by Marian Piatkowski.
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
* @see Dune::Dorie::VolumeRaviartThomasLocalFiniteElement.
*
* @tparam DF The domain field.
* @tparam RF The range field.
* @tparam k The order.
* @tparam d The dimension.
* @tparam gt The geometry type.
*/
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalBasis
{
static_assert(Dune::AlwaysFalse<DF>::value, "Volume Raviart Thomas basis function is not implemented!");
//! Identify when there is no implementation of this class.
static_assert(Dune::AlwaysFalse<DF>::value,
"Volume Raviart Thomas basis function is not implemented!");
#ifdef DOXYGEN /* Documentation of the code */
//! Order of the polynomials.
static constexpr int k;
//! World dimension.
static constexpr int dim;
//! Number of shape functions
static constexpr int volsize;
/*-----------------------------------------------------------------------*//**
* @brief Traits of local basis.
* @details It follows the traits structure for local basis defined by
* dune-localfunctions (Dune::LocalBasisTraits).
*/
struct Traits {};
/*-----------------------------------------------------------------------*//**
* @brief Number of shape functions (i.e. degrees of freedom).
*
* @return The number of shape functions.
*/
constexpr unsigned int size() const;
/*-----------------------------------------------------------------------*//**
* @brief Evaluates the each local base function (in the vector out) in
* the position in.
*
* @param[in] in Position to evaluate the base function in local
* coordinates.
* @param out Vector of values. Each value represents one local base
* function evaluated in the position in.
*/
inline void evaluateFunction(const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const;
/*-----------------------------------------------------------------------*//**
* @brief Evaluates the jacobian of each local base function (in the
* vector out) in the position in.
*
* @param[in] in Position to evaluate the base function in local
* coordinates.
* @param out Vector of jacobians. Each entry i in the vector
* represents the jacobian of the local base function i
* evaluated in the position in.
*/
inline void
evaluateJacobian(const typename Traits::DomainType& in,
std::vector<typename Traits::JacobianType>& out) const;
/*-----------------------------------------------------------------------*//**
* @brief Evaluate partial derivatives of any order of all shape
* functions.
*
* @param[in] order Order of the partial derivatives, in the classic
* multi-index notation.
* @param[in] in Position to evaluate the base function in local
* coordinates.
* @param out The desired partial derivatives.
*/
void partial(const std::array<unsigned int,dim>& order,
const typename Traits::DomainType& in,
std::vector<typename Traits::RangeType>& out) const;
/*-------------------------------------------------------------------------*//**
* @brief Polynomial order of the shape functions for quadrature
*
* @return The order.
*/
constexpr unsigned int order() const;
#endif /* Documentation of the code */
};
#ifndef DOXYGEN /* Actual code */
// This classes only implement the parts needed for the flux reconstruction.
// Never use them to generate a grid function space in dune-pdelab!
template<class DF, class RF>
class VolumeRaviartThomasLocalBasis<DF,RF,0,2,Dune::GeometryType::simplex>
{
......@@ -731,6 +827,8 @@ private:
Dune::QkStuff::GaussLobattoLagrangePolynomials<DF,RF,order-1> poly_km1;
};
#endif /* Actual code */
} // namespace Dorie
} // namespace Dune
......
......@@ -13,46 +13,44 @@
namespace Dune{
namespace Dorie{
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
/*-------------------------------------------------------------------------*//**
* @brief Class for local coefficients of the volume raviart thomas finite
* element.
* @details This class is used to reconstruct fluxes from local operators
* that solve fluxes in the residual vector. In particular, since
* the volume part is always associated to the codimension 0, all
* the local keys are the same.
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
* @see Dune::Dorie::VolumeRaviartThomasLocalFiniteElement.
*
* @tparam DF The domain field.
* @tparam RF The range field.
* @tparam k The order.
* @tparam dim The dimension.
* @tparam gt The geometry type.
*/
template<class DF, class RF, int k, int dim, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalCoefficients
{
static_assert(Dune::AlwaysFalse<DF>::value, "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:
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
for (unsigned int i = 0; i < size(); ++i)
li[i] = LocalKey(0,0,i);
}
constexpr std::size_t size() const
{
return volsize;
}
const LocalKey& localKey (std::size_t i) const
{
return li[i];
}
private:
std::vector<LocalKey> li;
};
template<class DF, class RF, int k, int dim>
class VolumeRaviartThomasLocalCoefficients<DF,RF,k,dim,Dune::GeometryType::cube>
{
static constexpr int volsize = dim*Dune::StaticPower<k+1,dim-1>::power*k;
//! Identify when there is no implementation of this class.
static_assert( (gt == Dune::GeometryType::simplex) ||
(gt == Dune::GeometryType::cube),
"Volume Raviart Thomas coefficients is not implemented!");
//! Number of shape functions
static constexpr int volsize = (gt == Dune::GeometryType::simplex)
? ((dim == 2) ? dim*k*(k+1)/2 : dim*k*(k+1)*(k+2)/6) // dim for simplex
: (dim*Dune::StaticPower<k+1,dim-1>::power*k); // dim for cubes
public:
/*-----------------------------------------------------------------------*//**
* @brief Constructs the coeffcients.
* @details All locall keys are associated to codimension 0.
*/
VolumeRaviartThomasLocalCoefficients()
: li(size())
{
......@@ -60,17 +58,32 @@ public:
li[i] = LocalKey(0,0,i);
}
/*-----------------------------------------------------------------------*//**
* @brief Number of local coefficients
*
* @return The size of the local coefficients
*/
constexpr std::size_t size() const
{
return volsize;
}
/**
* @brief Local key for a given index.
*
* @param[in] i The index that identifies a shape function (see
* Dune::Dorie::VolumeRaviartThomasLocalBasis).
*
* @return The local key associated to the index i. (see Dune:.LocalKey in
* dune-localfunctions)
*/
const LocalKey& localKey (std::size_t i) const
{
return li[i];
}
private:
//! The vector of local keys.
std::vector<LocalKey> li;
};
......
......@@ -8,6 +8,20 @@
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Class for local interpolation of the volume raviart thomas finite
* element.
* @details This class is used to reconstruct fluxes from local operators
* that solve fluxes in the residual vector. In particular, the
* content of the class is never needed for the reconstruction, and
* so, the usage of this class always throws an error.
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
* @see Dune::Dorie::VolumeRaviartThomasLocalFiniteElement.
*/
class VolumeRaviartThomasLocalInterpolation
{
public:
......
#ifndef DUNE_DORIE_FLUX_RECONSTRUCTION_VOLUME_FEM_HH
#define DUNE_DORIE_FLUX_RECONSTRUCTION_VOLUME_FEM_HH
#include <dune/geometry/type.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/volume_raviart_thomas.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Finite element map for volume part of the raviart thomas.
* @details This class basically provides a volume raviart thomas finite
* element for any entity following the finite element map defined
* by dune-pdelab (see
* Dune::PDELab::LocalFiniteElementMapInterface).
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
* @see Dune::Dorie::VolumeRaviartThomasLocalFiniteElement.
*
* @tparam DF The domain field.
* @tparam RF The range field.
* @tparam k The order.
* @tparam d The dimension.
* @tparam gt The geometry type.
*/
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalFiniteElementMap
: public Dune::PDELab::SimpleLocalFiniteElementMap<
VolumeRaviartThomasLocalFiniteElement<DF,RF,k,d,gt>,d>
{
public:
#ifdef DOXYGEN /* Documentation of the code */
/*-----------------------------------------------------------------------*//**
* @brief Return finite element for the given entity.
*
* @details The return value is a reference to Traits::LocalBasisType. If
* there is a different local basis 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()
{
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)
{
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
{
return (geometry_type.dim() == d) ? _fe.size(): 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
{
return _fe.size();
}
private:
//! The raviart thomas finite element.
const VolumeRaviartThomasLocalFiniteElement<DF,RF,k,d,gt> _fe;
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_FLUX_RECONSTRUCTION_VOLUME_FEM_HH
\ No newline at end of file
#ifndef DUNE_DORIE_VOLUME_RAVIART_THOMAS_HH
#define DUNE_DORIE_VOLUME_RAVIART_THOMAS_HH
#include "volume/volumertlocalbasis.hh"
#include "volume/volumertlocalcoefficients.hh"
#include "volume/volumertlocalinterpolation.hh"
#include <dune/geometry/type.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/volume/volumertlocalbasis.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/volume/volumertlocalcoefficients.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/volume/volumertlocalinterpolation.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Class for volume part of the raviart thomas finite element.
* @details This class is used to reconstruct fluxes from local operators
* that solve fluxes in the residual vector. In particular, this
* volume raviart thomas is equivalent to the usual raviart thomas
* element but excluding the degrees of freedom located on the faces
* (i.e. codim = 1).
*
* @author Santiago Ospina De Los Ríos
* @date 2018
* @copyright MIT License.
* @ingroup FluxReconstruction
* @see Dune::Dorie::VolumeRaviartThomasLocalBasis.
* @see Dune::Dorie::VolumeRaviartThomasLocalCoefficients.
* @see Dune::Dorie::VolumeRaviartThomasLocalInterpolation.
*
* @tparam DF The domain field.
* @tparam RF The range field.
* @tparam k The order.
* @tparam d The dimension.
* @tparam gt The geometry type.
*/
template<class DF, class RF, int k, int d, Dune::GeometryType::BasicType gt>
class VolumeRaviartThomasLocalFiniteElement
{
public:
//! Traits of local basis. It follows the traits structure for finite
//! elements defined by dune-localfunctions (Dune::LocalFiniteElementTraits).
using Traits = Dune::LocalFiniteElementTraits<
VolumeRaviartThomasLocalBasis<DF,RF,k,d,gt>,
VolumeRaviartThomasLocalCoefficients<DF,RF,k,d,gt>,
VolumeRaviartThomasLocalInterpolation>;
/*-----------------------------------------------------------------------*//**
* @brief The local basis for this finite element.
* @see Dune::Dorie::VolumeRaviartThomasLocalBasis.
*
* @return The local basis.