Commit 50a8063d authored by Jakub Both's avatar Jakub Both

Allow RT elements to be defined on arbitrary cubes.

Reuse finiteelementmap from PDELab to handle different variants
of local Raviart-Thomas elements. Due to missing dependencies
several code segments have been copied from PDELab.

Note: Tests show, that RT0 basis functions are continuous in
normal direction across edges but do not have the same sign on
edges. There are basis functions with value +1 on edges and other
with -1.
parent edf2d01d
......@@ -14,6 +14,7 @@
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
#include <dune/functions/pdelab-copy/raviartthomasfem.hh>
namespace Dune {
namespace Functions {
......@@ -88,7 +89,7 @@ public:
template<class TP>
Node<TP> node(const TP& tp) const
{
return Node<TP>{tp};
return Node<TP>{tp, gridView_};
  • This is now a O(n) operation. IMHO all informations that require O(n) setup should be stored in the factory to always keep this O(1) (where the constant can probably be large).

  • I have tried to follow your suggestion. Is it better now?

Please register or sign in to reply
}
template<class TP>
......@@ -142,12 +143,14 @@ public:
using size_type = ST;
using TreePath = TP;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = RaviartThomasCubeLocalFiniteElement<typename GV::ctype, double, dim, k>; // general RTk implementation on cube grids
using FiniteElementMap = typename PDELab::detail::RaviartThomasLocalFiniteElementMapBaseSelector<GV, dim, Dune::GeometryType::cube, typename GV::ctype, double, k>::type;
using FiniteElement =typename FiniteElementMap::Traits::FiniteElement;
RaviartThomasCubeNode(const TreePath& treePath) :
RaviartThomasCubeNode(const TreePath& treePath, const GV& gv) :
Base(treePath),
finiteElement_(nullptr),
element_(nullptr)
element_(nullptr),
finiteElementMap_(gv)
{ }
//! Return current element, throw if unbound
......@@ -169,17 +172,15 @@ public:
void bind(const Element& e)
{
element_ = &e;
finiteElement_ = new FiniteElement(5); // get local finite element
finiteElement_ = &(finiteElementMap_.find(*element_));
this->setSize(finiteElement_->size());
// TODO: The use of the constructor is hardcoded. This does not always have to work!
// TODO: See PDELAB::finiteelementmap.hh -> RTLocalFiniteElementMap, uses gridView.
}
protected:
const FiniteElement* finiteElement_;
const Element* element_;
FiniteElementMap finiteElementMap_;
};
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PDELAB_COMMON_EXCEPTIONS_HH
#define DUNE_PDELAB_COMMON_EXCEPTIONS_HH
#include <dune/common/exceptions.hh>
/**
* \file
* \brief PDELab-specific exceptions.
*/
namespace Dune {
namespace PDELab {
//! Base class for all PDELab exceptions.
class Exception
: public Dune::Exception
{};
//! GridFunctionSpace-related error.
class GridFunctionSpaceError
: public Exception
{};
//! Called a GridFunctionSpace method that requires initialization of the space.
class UninitializedGridFunctionSpaceError
: public GridFunctionSpaceError
{};
//! Called a method on a GridFunctionSpace that is not valid
//! at its current place in the function space tree.
class GridFunctionSpaceHierarchyError
: public GridFunctionSpaceError
{};
//! Ordering-related error.
class OrderingError
: public Exception
{};
//! Error related to the logical structure of an Ordering.
class OrderingStructureError
: public OrderingError
{};
//! A PermutedOrdering got a permutation vector of the wrong size.
class PermutedOrderingSizeError
: public OrderingError
{};
//! The block size of a ChunkedBlockOrdering does not divide the block count of the underlying ordering.
class ChunkedBlockOrderingSizeError
: public OrderingError
{};
} // namespace PDELab
} // namespace Dune
#endif // DUNE_PDELAB_COMMON_EXCEPTIONS_HH
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_PDELAB_FINITELEMENTMAP_HH
#define DUNE_PDELAB_FINITELEMENTMAP_HH
#include <dune/common/deprecated.hh>
#include "exceptions.hh"
#include <dune/geometry/referenceelements.hh>
namespace Dune {
namespace PDELab {
//! \ingroup FiniteElementMap
//! \{
//! \brief general FiniteElementMap exception
class FiniteElementMapError : public Exception {};
//! \brief FiniteElementMap exception concerning the computation of the FiniteElementMap size
class VariableElementSize : public FiniteElementMapError {};
//! \brief FiniteElementMap exception raised when trying to obtain a finite element for an unsupported GeometryType.
class InvalidGeometryType : public FiniteElementMapError {};
//! \brief collect types exported by a finite element map
template<class T>
struct FiniteElementMapTraits
{
//! Type of finite element from local functions
typedef T FiniteElementType;
//! Type of finite element from local functions
typedef T FiniteElement;
};
//! \brief collect types exported by a finite element map
template<class T>
struct LocalFiniteElementMapTraits : FiniteElementMapTraits<T> {};
//! \brief interface for a finite element map
template<class T, class Imp>
class LocalFiniteElementMapInterface
{
public:
//! \brief Export traits
typedef T Traits;
/** \brief Return local basis for the given entity.
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.
*/
template<class EntityType>
const typename Traits::FiniteElementType& find (const EntityType& e) const
{
return asImp().find(e);
}
/** @name Size calculation
* The FiniteElementMap provides different methods to compute
* the size of the GridFunctionSpace (if possible) without
* iterating the grid. The approach is as follows (pseudo code):
*
* \code
* computeNumberOfDofs(GridView, FEM):
* if(FEM.fixedSize()):
* sum(FEM.size(gt)*GridView.size(gt) for gt in GeometryTypes)
* else
* sum(FEM.find(E).basis().size() for E in GridView.entities<0>())
* \endcode
*/
/** @{ */
/** \brief a FiniteElementMap is fixedSize iif the size of the local
* functions space for each GeometryType is fixed.
*/
bool fixedSize() const
{
return asImp().fixedSize();
}
/** \brief if the FiniteElementMap is fixedSize, the size
* methods computes the number of DOFs for given GeometryType.
*/
std::size_t size(GeometryType gt) const
{
return asImp().size(gt);
}
/** @} */
/** \brief return if FiniteElementMap has degrees of freedom for
* given codimension
*/
bool hasDOFs(int codim) const
{
return asImp().hasDOFs(codim);
}
/** \brief compute an upper bound for the local number of DOFs.
*
* this upper bound is used to avoid reallocations in
* std::vectors used during the assembly.
*/
std::size_t maxLocalSize() const
{
return asImp().maxLocalSize();
}
private:
Imp& asImp () {return static_cast<Imp &> (*this);}
const Imp& asImp () const {return static_cast<const Imp &>(*this);}
};
//! simple implementation where all entities have the same finite element
template<class Imp>
class SimpleLocalFiniteElementMap :
public LocalFiniteElementMapInterface<LocalFiniteElementMapTraits<Imp>,
SimpleLocalFiniteElementMap<Imp> >
{
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<Imp> Traits;
//! \brief Use when Imp has a standard constructor
SimpleLocalFiniteElementMap ()
{}
//! \brief Constructor where an instance of Imp can be provided
SimpleLocalFiniteElementMap (const Imp& imp_) : imp(imp_)
{}
//! \brief get local basis functions for entity
template<class EntityType>
const typename Traits::FiniteElementType& find (const EntityType& e) const
{
return imp;
}
private:
Imp imp; // create once
};
/** \brief implementation for finite elements requiring oriented edges
*
* This is for edge elements. It works for one type of Geometry only, and
* the requirements for the local finite element are:
*
* - The default orientation of the shape functions on the edges must be
* from the vertex with the lower index within the reference element to
* the vertex with the higher index.
* - The local finite element must allow assignment.
* - The local finite element must have a default constructor.
* - the local finite element hust have a constructor of the form
* FE(unsigned int orientation), where orientation is a bitfield. If
* bit i in orientation is set it means that the orientation of the
* shape function corresponding to the edge with id i in the reference
* element is inverted.
*
* \tparam GV Type of gridview to work with
* \tparam FE Type of local finite element
* \tparam Imp Type of the final LocalFiniteElementMap implementation
*/
template<typename GV, typename FE, typename Imp>
class EdgeS0LocalFiniteElementMap
: public LocalFiniteElementMapInterface<
LocalFiniteElementMapTraits<FE>,
Imp
>
{
typedef typename GV::IndexSet IndexSet;
static const int dim = GV::dimension;
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<FE> Traits;
//! \brief construct EdgeSLocalFiniteElementMap
EdgeS0LocalFiniteElementMap (const GV& gv_)
: gv(gv_), orient(gv_.size(0))
{
typedef typename GV::Grid::ctype ct;
auto& refElem =
ReferenceElements<ct, dim>::general(FE().type());
auto &idSet = gv.grid().globalIdSet();
// create all variants
variant.resize(1 << refElem.size(dim-1));
for (unsigned int i=0; i<variant.size(); i++)
variant[i] = FE(i);
// compute orientation for all elements
auto& indexSet = gv.indexSet();
// loop once over the grid
for (const auto& element : elements(gv))
{
auto elemid = indexSet.index(element);
orient[elemid] = 0;
std::vector<typename GV::Grid::GlobalIdSet::IdType> vid(refElem.size(dim));
for(unsigned int i = 0; i < vid.size(); ++i)
vid[i] = idSet.subId(element, i, dim);
// loop over all edges of the element
for(int i = 0; i < refElem.size(dim-1); ++i) {
int v0 = refElem.subEntity(i, dim-1, 0, dim);
int v1 = refElem.subEntity(i, dim-1, 1, dim);
// if (edge orientation in refelement) != (edge orientation in indexset)
if((v0 > v1) != (vid[v0] > vid[v1]))
orient[elemid] |= 1 << i;
}
}
}
//! \brief get local basis functions for entity
template<class EntityType>
const typename Traits::FiniteElementType& find (const EntityType& e) const
{
return variant[orient[gv.indexSet().index(e)]];
}
private:
GV gv;
std::vector<FE> variant;
std::vector<unsigned char> orient;
};
//! wrap up element from local functions
//! \ingroup FiniteElementMap
template<typename GV, typename FE, typename Imp, std::size_t Variants>
class RTLocalFiniteElementMap :
public LocalFiniteElementMapInterface<
LocalFiniteElementMapTraits<FE>,
Imp >
{
typedef FE FiniteElement;
typedef typename GV::IndexSet IndexSet;
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<FE> Traits;
//! \brief Use when Imp has a standard constructor
RTLocalFiniteElementMap(const GV& gv_)
: gv(gv_), is(gv_.indexSet()), orient(gv_.size(0))
{
// create all variants
for (int i = 0; i < Variants; i++)
{
variant[i] = FiniteElement(i);
}
// compute orientation for all elements
// loop once over the grid
for(const auto& cell : elements(gv))
{
unsigned int myId = is.index(cell);
orient[myId] = 0;
for (const auto& intersection : intersections(gv,cell))
{
if (intersection.neighbor()
&& is.index(intersection.outside()) > myId)
{
orient[myId] |= 1 << intersection.indexInInside();
}
}
}
}
//! \brief get local basis functions for entity
template<class EntityType>
const typename Traits::FiniteElementType& find(const EntityType& e) const
{
return variant[orient[is.index(e)]];
}
private:
GV gv;
std::array<FiniteElement,Variants> variant;
const IndexSet& is;
std::vector<unsigned char> orient;
};
//! \} group FiniteElementMap
} // namespace PDELab
} // namespace Dune
#endif
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_RAVIARTTHOMASFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_RAVIARTTHOMASFEM_HH
#include <dune/grid/common/capabilities.hh>
#include "topologyutility.hh"
#include "rt0simplex2dfem.hh"
#include "rt1simplex2dfem.hh"
#include "rt0cube2dfem.hh"
#include "rt1cube2dfem.hh"
#include "rt2cube2dfem.hh"
#include "rt0cube3dfem.hh"
#include "rt1cube3dfem.hh"
namespace Dune {
namespace PDELab {
#ifndef DOXYGEN
namespace detail {
//! Template for determining the correct base implementation of RaviartThomasLocalFiniteElementMap.
/**
* This template must be specialized for supported variations of the Raviart-Thomas elements.
* Its member type 'type' will be used as base class of RaviartThomasLocalFiniteElementMap. That
* type must have a constructor that takes the GridView as single parameter.
*/
template<typename GV, int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
struct RaviartThomasLocalFiniteElementMapBaseSelector
{
static_assert((AlwaysFalse<GV>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
};
// ********************************************************************************
// Specializations
// ********************************************************************************
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,2,GeometryType::simplex,D,R,0>
{
typedef RT0Simplex2DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,2,GeometryType::simplex,D,R,1>
{
typedef RT1Simplex2DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,2,GeometryType::cube,D,R,0>
{
typedef RT0Cube2DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,2,GeometryType::cube,D,R,1>
{
typedef RT1Cube2DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,2,GeometryType::cube,D,R,2>
{
typedef RT2Cube2DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,3,GeometryType::cube,D,R,0>
{
typedef RT0Cube3DLocalFiniteElementMap<GV,D,R> type;
};
template<typename GV, typename D, typename R>
struct RaviartThomasLocalFiniteElementMapBaseSelector<GV,3,GeometryType::cube,D,R,1>
{
typedef RT1Cube3DLocalFiniteElementMap<GV,D,R> type;
};
} // end namespace detail
#endif // DOXYGEN
//! Raviart-Thomas elements of order k.
/**
* This LocalFiniteElementMap provides Raviart-Thomas elements of order k for the given
* GridView GV. It currently supports
*
* * RT0, RT1 and RT2 for cubes in 2D,
* * RT0 and RT1 for simplices in 2D,
* * RT0 and RT1 for cubes in 3D.
*
* We try to infer the type of the reference element (cube / simplex) from the GridView, but
* that only works for grids with a single element type that is fixed at compile time. For
* potentially mixed grids like UGGrid, you need to provide the GeometryType::BasicType of the
* cell reference element as an additional template parameter.
*
* \tparam GV The GridView on which to construct the finite element map.
* \tparam D The domain field type of the elements.
* \tparam R The range field type of the elements.
* \tparam k The order of the finite elements.
* \tparam basic_type The GeometryType::BasicType of the grid cells. You only need to provide this
* template parameter for mixed grids (if you don't provide the parameter for a
* mixed grid, you get a compiler error).
*/
template<typename GV,
typename D,
typename R,
std::size_t k,
GeometryType::BasicType basic_type = BasicTypeFromDimensionAndTopologyId<
GV::dimension,
Capabilities::hasSingleGeometryType<typename GV::Grid>::topologyId
>::value
>
class RaviartThomasLocalFiniteElementMap :
public detail::RaviartThomasLocalFiniteElementMapBaseSelector<GV,GV::dimension,basic_type,D,R,k>::type
{
public:
//! Constructs a finite element map on the GridView gv.
RaviartThomasLocalFiniteElementMap(const GV& gv)
: detail::RaviartThomasLocalFiniteElementMapBaseSelector<GV,GV::dimension,basic_type,D,R,k>::type(gv)
{}
};
#ifndef DOXYGEN
// Specialization for grids that don't provide a valid topology id for their cells.
template<typename GV, typename D, typename R, std::size_t k>
class RaviartThomasLocalFiniteElementMap<GV,D,R,k,GeometryType::none>
{
static_assert((AlwaysFalse<GV>::value),
"Your chosen grid does not export a usable topology id for its cells."
"Please provide the correct GeometryType::BasicType as an additional template parameter.");
};
#endif // DOXYGEN
} // end namespace PDELab
} // end namespace Dune
#endif // DUNE_PDELAB_FINITEELEMENTMAP_RAVIARTTHOMASFEM_HH
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE2DFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE2DFEM_HH
#include<vector>
#include<dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
#include"finiteelementmap.hh"
namespace Dune {
namespace PDELab {
//! wrap up element from local functions
//! \ingroup FiniteElementMap
template<typename GV, typename D, typename R>
class RT0Cube2DLocalFiniteElementMap :
public RTLocalFiniteElementMap<
GV,
Dune::RT0Cube2DLocalFiniteElement<D,R>,
RT0Cube2DLocalFiniteElementMap<GV,D,R>,
16>
{
typedef Dune::RT0Cube2DLocalFiniteElement<D,R> FE;
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<FE> Traits;
//! \brief Use when Imp has a standard constructor
RT0Cube2DLocalFiniteElementMap (const GV& gv)
: RTLocalFiniteElementMap<
GV,
Dune::RT0Cube2DLocalFiniteElement<D,R>,
RT0Cube2DLocalFiniteElementMap<GV,D,R>,
16>(gv)
{}
bool fixedSize() const
{
return true;
}
bool hasDOFs(int codim) const
{
return codim == 1;
}
std::size_t size(GeometryType gt) const
{
return gt.isLine() ? 1 : 0;
}
std::size_t maxLocalSize() const
{
return 4;
}
};
} // end namespace PDELab
} // end namespace Dune
#endif // DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE2DFEM_HH
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE3DFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE3DFEM_HH
#include<vector>
#include<dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
#include"finiteelementmap.hh"
namespace Dune {
namespace PDELab {
//! wrap up element from local functions
//! \ingroup FiniteElementMap
template<typename GV, typename D, typename R>
class RT0Cube3DLocalFiniteElementMap :
public RTLocalFiniteElementMap<
GV,
Dune::RT0Cube3DLocalFiniteElement<D,R>,
RT0Cube3DLocalFiniteElementMap<GV,D,R>,
64>
{
typedef Dune::RT0Cube3DLocalFiniteElement<D,R> FE;
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<FE> Traits;
//! \brief Use when Imp has a standard constructor
RT0Cube3DLocalFiniteElementMap (const GV& gv)
: RTLocalFiniteElementMap<
GV,
Dune::RT0Cube3DLocalFiniteElement<D,R>,
RT0Cube3DLocalFiniteElementMap<GV,D,R>,
64>(gv)
{}
bool fixedSize() const
{
return true;
}
bool hasDOFs(int codim) const
{
return codim == 1;
}
std::size_t size(GeometryType gt) const
{
return gt.dim() == 2 && gt.isCube() ? 1 : 0;
}
std::size_t maxLocalSize() const
{
return 6;
}
};
} // end namespace PDELab
} // end namespace Dune
#endif // DUNE_PDELAB_FINITEELEMENTMAP_RT0CUBE3DFEM_HH
// -*- tab-width: 4; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_FINITEELEMENTMAP_RT0SIMPLEX2DFEM_HH
#define DUNE_PDELAB_FINITEELEMENTMAP_RT0SIMPLEX2DFEM_HH
#include<vector>
#include<dune/localfunctions/raviartthomas/raviartthomas02d.hh>
#include"finiteelementmap.hh"
namespace Dune {
namespace PDELab {
//! wrap up element from local functions
//! \ingroup FiniteElementMap
template<typename GV, typename D, typename R>
class RT0Simplex2DLocalFiniteElementMap :
public RTLocalFiniteElementMap<
GV,
Dune::RT02DLocalFiniteElement<D,R>,
RT0Simplex2DLocalFiniteElementMap<GV,D,R>,
8>
{
typedef Dune::RT02DLocalFiniteElement<D,R> FE;
public:
//! \brief export type of the signature
typedef LocalFiniteElementMapTraits<FE> Traits;
//! \brief Use when Imp has a standard constructor
RT0Simplex2DLocalFiniteElementMap (const GV& gv)
: RTLocalFiniteElementMap<
GV,
Dune::RT02DLocalFiniteElement<D,R>,
RT0Simplex2DLocalFiniteElementMap<GV,D,R>,
8> (gv)
{}