...
 
Commits (33)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINICUBEBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINICUBEBASIS_HH
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube2d.hh>
#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1cube3d.hh>
#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini1simplex2d.hh>
#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2cube2d.hh>
#include <dune/localfunctions/brezzidouglasmarini/brezzidouglasmarini2simplex2d.hh>
#include <dune/typetree/leafnode.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
namespace Dune {
namespace Functions {
namespace Impl {
template<int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
struct BDMLocalInfo
{
static_assert((AlwaysFalse<D>::value),"The requested type of BDM element is not implemented, sorry!");
};
template<typename D, typename R>
struct BDMLocalInfo<2,GeometryType::simplex,D,R,1>
{
using FiniteElement = BDM1Simplex2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct BDMLocalInfo<2,GeometryType::simplex,D,R,2>
{
using FiniteElement = BDM2Simplex2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct BDMLocalInfo<2,GeometryType::cube,D,R,1>
{
using FiniteElement = BDM1Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct BDMLocalInfo<2,GeometryType::cube,D,R,2>
{
using FiniteElement = BDM2Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct BDMLocalInfo<3,GeometryType::cube,D,R,1>
{
using FiniteElement = BDM1Cube3DLocalFiniteElement<D,R>;
static const std::size_t Variants = 64;
};
template<typename GV, int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
class BDMLocalFiniteElementMap
{
static const std::size_t Variants = BDMLocalInfo<dim, basic_type, D, R, k>::Variants;
public:
using FiniteElement = typename BDMLocalInfo<dim, basic_type, D, R, k>::FiniteElement;
BDMLocalFiniteElementMap(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 FiniteElement& find(const EntityType& e) const
{
return variant_[orient_[is_->index(e)]];
}
private:
GV gv_;
std::array<FiniteElement,Variants> variant_;
const typename GV::IndexSet* is_;
std::vector<unsigned char> orient_;
};
} // namespace Impl
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
// BrezziDouglasMariniNodeFactory
// BrezziDouglasMariniNodeIndexSet
// BrezziDouglasMariniNode
//
// The factory allows to create the others and is the owner of possible shared
// state. These three components do _not_ depend on the global basis or index
// set and can be used without a global basis.
// *****************************************************************************
template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNode;
template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNodeIndexSet;
template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNodeFactory;
template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNodeFactory
{
static const int dim = GV::dimension;
using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
public:
/** \brief The grid view that the FE space is defined on */
using GridView = GV;
using size_type = ST;
// Precompute the number of dofs per entity type depending on the entity's codimension
std::vector<int> dofsPerCodim {dim*(k-1)*3, dim+(k-1)}; // hardcoded
template<class TP>
using Node = BrezziDouglasMariniNode<GV, k, size_type, TP, basic_type>;
template<class TP>
using IndexSet = BrezziDouglasMariniNodeIndexSet<GV, k, MI, TP, ST, basic_type>;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using SizePrefix = Dune::ReservedVector<size_type, 1>;
/** \brief Constructor for a given grid view object */
BrezziDouglasMariniNodeFactory(const GridView& gv) :
gridView_(gv),
finiteElementMap_(gv)
{ }
void initializeIndices()
{
CodimOffset_.resize(3);
CodimOffset_[0] = 0;
CodimOffset_[1] = CodimOffset_[0] + dofsPerCodim[0] * gridView_.size(0);
//if (dim==3) CodimOffset_[2] = CodimOffset_[1] + dofsPerCodim[1] * gridView_.size(1);
}
/** \brief Obtain the grid view that the basis is defined on
*/
const GridView& gridView() const
{
return gridView_;
}
template<class TP>
Node<TP> node(const TP& tp) const
{
return Node<TP>{tp, &finiteElementMap_};
}
template<class TP>
IndexSet<TP> indexSet() const
{
return IndexSet<TP>{*this};
}
size_type size() const
{
return dofsPerCodim[0] * gridView_.size(0) + dofsPerCodim[1] * gridView_.size(1); // only 2d
}
//! Return number possible values for next position in multi index
size_type size(const SizePrefix prefix) const
{
assert(prefix.size() == 0 || prefix.size() == 1);
return (prefix.size() == 0) ? size() : 0;
}
/** \todo This method has been added to the interface without prior discussion. */
size_type dimension() const
{
return size();
}
size_type maxNodeSize() const
{
return StaticPower<(k+1),GV::dimension>::power;
}
//protected:
const GridView gridView_;
std::vector<size_t> CodimOffset_;
FiniteElementMap finiteElementMap_;
};
template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNode :
public LeafBasisNode<ST, TP>
{
static const int dim = GV::dimension;
static const int maxSize = StaticPower<(k+1),GV::dimension>::power;
using Base = LeafBasisNode<ST,TP>;
public:
using size_type = ST;
using TreePath = TP;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElementMap = typename Impl::BDMLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
using FiniteElement = typename FiniteElementMap::FiniteElement;
BrezziDouglasMariniNode(const TreePath& treePath, const FiniteElementMap* finiteElementMap) :
Base(treePath),
finiteElement_(nullptr),
element_(nullptr),
finiteElementMap_(finiteElementMap)
{ }
//! Return current element, throw if unbound
const Element& element() const
{
return *element_;
}
/** \brief Return the LocalFiniteElement for the element we are bound to
*
* The LocalFiniteElement implements the corresponding interfaces of the dune-localfunctions module
*/
const FiniteElement& finiteElement() const
{
return *finiteElement_;
}
//! Bind to element.
void bind(const Element& e)
{
element_ = &e;
finiteElement_ = &(finiteElementMap_->find(*element_));
this->setSize(finiteElement_->size());
}
protected:
const FiniteElement* finiteElement_;
const Element* element_;
const FiniteElementMap* finiteElementMap_;
};
template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
class BrezziDouglasMariniNodeIndexSet
{
enum {dim = GV::dimension};
public:
using size_type = ST;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using NodeFactory = BrezziDouglasMariniNodeFactory<GV, k, MI, ST, basic_type>;
using Node = typename NodeFactory::template Node<TP>;
BrezziDouglasMariniNodeIndexSet(const NodeFactory& nodeFactory) :
nodeFactory_(&nodeFactory)
{}
/** \brief Bind the view to a grid element
*
* Having to bind the view to an element before being able to actually access any of its data members
* offers to centralize some expensive setup code in the 'bind' method, which can save a lot of run-time.
*/
void bind(const Node& node)
{
node_ = &node;
}
/** \brief Unbind the view
*/
void unbind()
{
node_ = nullptr;
}
/** \brief Size of subtree rooted in this node (element-local)
*/
size_type size() const
{
return node_->finiteElement().size();
}
//! Maps from subtree index set [0..size-1] to a globally unique multi index in global basis
//! Assume dim \in \lbrace 2, 3 \rbrace.
MultiIndex index(size_type i) const
{
Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
const auto& gridIndexSet = nodeFactory_->gridView().indexSet();
const auto& element = node_->element();
// The dimension of the entity that the current dof is related to
size_t subentity = localKey.subEntity();
size_t codim = localKey.codim();
// Throw if Element is no cube or simplex
if (not(basic_type==GeometryType::BasicType::cube and element.type().isCube()) and
not(basic_type==GeometryType::BasicType::simplex and element.type().isSimplex())) DUNE_THROW(Dune::NotImplemented, "BDMNodalBasis only implemented for cube and simplex elements.");
if (not(codim==0 or codim==1)) DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the BDMThomasBasis");
// if (codim==0) { // element dof
return { nodeFactory_->CodimOffset_[codim] +
nodeFactory_->dofsPerCodim[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
// }
// // treat edges individually due to orientation
// if (codim==1) { // edge/face dof
// const Dune::ReferenceElement<double,dim>& refElement
// = Dune::ReferenceElements<double,dim>::general(element.type());
//
// // we have to reverse the numbering if the local triangle edge is
// // not aligned with the global edge
// size_t v0 = gridIndexSet.subIndex(element,refElement.subEntity(subentity,codim,0,dim),dim);
// size_t v1 = gridIndexSet.subIndex(element,refElement.subEntity(subentity,codim,1,dim),dim);
// bool flip = (v0 > v1);
// //if (k==0) flip=true;
// return { (flip)
// ? nodeFactory_->CodimOffset_[codim]
// + nodeFactory_->dofsPerCodim[codim] * gridIndexSet.subIndex(element,subentity,codim) + (nodeFactory_->dofsPerCodim[codim]-1) - localKey.index()
// : nodeFactory_->CodimOffset_[codim]
// + nodeFactory_->dofsPerCodim[codim] * gridIndexSet.subIndex(element,subentity,codim) + localKey.index() };
// }
// DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the BrezziDouglasMariniCubeBasis");
}
protected:
const NodeFactory* nodeFactory_;
const Node* node_;
};
namespace BasisBuilder {
namespace Imp {
template<std::size_t k, GeometryType::BasicType basic_type>
struct BrezziDouglasMariniNodeFactoryBuilder
{
static const std::size_t requiredMultiIndexSize=1;
template<class MultiIndex, class GridView, class size_type=std::size_t>
auto build(const GridView& gridView)
-> BrezziDouglasMariniNodeFactory<GridView, k, MultiIndex, size_type, basic_type>
{
return {gridView};
}
};
} // end namespace BasisBuilder::Imp
template<std::size_t k, GeometryType::BasicType basic_type>
Imp::BrezziDouglasMariniNodeFactoryBuilder<k, basic_type> bdm()
{
return{};
}
} // end namespace BasisBuilder
// *****************************************************************************
// This is the actual global basis implementation based on the reusable parts.
// *****************************************************************************
/** \brief Nodal basis of a scalar k-th-order BDM finite element space on simplicial elements
*
* TODO
*
* \tparam GV The GridView that the space is defined on
* \tparam k The order of the basis
*/
template<typename GV, int k, GeometryType::BasicType basic_type, class ST = std::size_t>
using BrezziDouglasMariniCubeNodalBasis = DefaultGlobalBasis<BrezziDouglasMariniNodeFactory<GV, k, FlatMultiIndex<ST>, ST, basic_type> >;
} // end namespace Functions
} // end namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_BREZZIDOUGLASMARINICUBEBASIS_HH
......@@ -75,6 +75,32 @@ struct FlatVectorBackend<typename Dune::FieldMatrix<K, n, m> >
};
template<class K, int n, int m>
struct FlatVectorBackend<typename Dune::FieldVector< typename Dune::FieldVector<K, n>, m> >
{
// Note: The ordering of first and second index is reverse compared to the FieldMatrix
// specialization. This is due to the 'wrong' definition of the JacobianRange in
// discreteglobalbasisfunction.hh. When fixing that issue, the order here has to be changed
// as well.
template<class VV, class Index>
static auto getEntry(VV&& v, const Index& i) -> decltype(v[i%m][i/m])
{
return v[i%m][i/m];
}
template<class VV>
static std::size_t size(VV&& v)
{
auto size = 0;
for (std::size_t i=0; i<v.size(); ++i)
size += v[i].size();
return size;
}
};
} // namespace Dune::Functions
} // namespace Dune
......
......@@ -110,6 +110,7 @@ public:
{
using FiniteElement = typename Node::FiniteElement;
using FiniteElementRange = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<FiniteElement>::type;
// Note that we capture j by reference. Hence we can switch
......@@ -126,7 +127,7 @@ public:
using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<FiniteElementRange(LocalDomain), decltype(localFj), FunctionBaseClass>;
auto interpolationValues = std::vector<FiniteElementRange>();
auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
auto&& fe = node.finiteElement();
......@@ -139,7 +140,7 @@ public:
{
// We cannot use localFj directly because interpolate requires a Dune::VirtualFunction like interface
fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationValues);
fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationCoefficients);
for (size_t i=0; i<fe.localBasis().size(); ++i)
{
auto multiIndex = localIndexSet_.index(node.localIndex(i));
......@@ -149,7 +150,7 @@ public:
if (interpolateHere)
{
auto&& vectorBlock = vector_[multiIndex];
FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationValues[i];
FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationCoefficients[i];
}
}
}
......
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/raviartthomas.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas0cube2d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas0cube3d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas02d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas1cube2d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas1cube3d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas12d.hh>
#include <dune/localfunctions/raviartthomas/raviartthomas2cube2d.hh>
#include <dune/typetree/leafnode.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
namespace Dune {
namespace Functions {
namespace Impl {
template<int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
struct RaviartThomasLocalInfo
{
static_assert((AlwaysFalse<D>::value),"The requested type of Raviart-Thomas element is not implemented, sorry!");
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<2,GeometryType::simplex,D,R,0>
{
using FiniteElement = RT02DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<2,GeometryType::simplex,D,R,1>
{
using FiniteElement = RT12DLocalFiniteElement<D,R>;
static const std::size_t Variants = 8;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,0>
{
using FiniteElement = RT0Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,1>
{
using FiniteElement = RT1Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<2,GeometryType::cube,D,R,2>
{
using FiniteElement = RT2Cube2DLocalFiniteElement<D,R>;
static const std::size_t Variants = 16;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<3,GeometryType::cube,D,R,0>
{
using FiniteElement = RT0Cube3DLocalFiniteElement<D,R>;
static const std::size_t Variants = 64;
};
template<typename D, typename R>
struct RaviartThomasLocalInfo<3,GeometryType::cube,D,R,1>
{
using FiniteElement = RT1Cube3DLocalFiniteElement<D,R>;
static const std::size_t Variants = 64;
};
template<typename GV, int dim, GeometryType::BasicType basic_type, typename D, typename R, std::size_t k>
class RaviartThomasLocalFiniteElementMap
{
static const std::size_t Variants = RaviartThomasLocalInfo<dim, basic_type, D, R, k>::Variants;
public:
using FiniteElement = typename RaviartThomasLocalInfo<dim, basic_type, D, R, k>::FiniteElement;
RaviartThomasLocalFiniteElementMap(const GV& gv)
: gv_(gv), is_(&(gv_.indexSet())), orient_(gv.size(0))
{
// create all variants
for (size_t 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 FiniteElement& find(const EntityType& e) const
{
return variant_[orient_[is_->index(e)]];
}
private:
GV gv_;
std::array<FiniteElement,Variants> variant_;
const typename GV::IndexSet* is_;
std::vector<unsigned char> orient_;
};
} // namespace Impl
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
// RaviartThomasNodeFactory
// RaviartThomasNodeIndexSet
// RaviartThomasNode
//
// The factory allows to create the others and is the owner of possible shared
// state. These three components do _not_ depend on the global basis or index
// set and can be used without a global basis.
// *****************************************************************************
template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
class RaviartThomasNode;
template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
class RaviartThomasNodeIndexSet;
template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
class RaviartThomasNodeFactory;
template<typename GV, int k, class MI, class ST, GeometryType::BasicType basic_type>
class RaviartThomasNodeFactory
{
static const int dim = GV::dimension;
using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
private:
template<typename, int, class, class, class, GeometryType::BasicType>
friend class RaviartThomasNodeIndexSet;
public:
/** \brief The grid view that the FE space is defined on */
using GridView = GV;
using size_type = ST;
// Precompute the number of dofs per entity type depending on the entity's codimension and the grid's type (only valid for cube and simplex grids)
// for 3D only for cubes k=0,1
// Note: dofsPerElement = dofsPerFace * dim for cubes, k=0,1
const static int dofsPerFace = dim == 2 ? k+1 : 3*k+1;
const static int dofsPerElement = dim == 2 ? (basic_type == GeometryType::cube ? k*(k+1)*dim : k*dim) : k*(k+1)*(k+1)*dim;
const std::vector<int> dofsPerCodim {dofsPerElement, dofsPerFace};
template<class TP>
using Node = RaviartThomasNode<GV, k, size_type, TP, basic_type>;
template<class TP>
using IndexSet = RaviartThomasNodeIndexSet<GV, k, MI, TP, ST, basic_type>;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using SizePrefix = Dune::ReservedVector<size_type, 1>;
/** \brief Constructor for a given grid view object */
RaviartThomasNodeFactory(const GridView& gv) :
gridView_(gv),
finiteElementMap_(gv)
{ }
void initializeIndices()
{
codimOffset_.resize(2);
codimOffset_[0] = 0;
codimOffset_[1] = codimOffset_[0] + dofsPerCodim[0] * gridView_.size(0);
}
/** \brief Obtain the grid view that the basis is defined on
*/
const GridView& gridView() const
{
return gridView_;
}
/* \brief Update the stored grid view, to be called if the grid has changed */
void update (const GridView& gv)
{
gridView_ = gv;
}
template<class TP>
Node<TP> node(const TP& tp) const
{
return Node<TP>{tp, &finiteElementMap_};
}
template<class TP>
IndexSet<TP> indexSet() const
{
return IndexSet<TP>{*this};
}
size_type size() const
{
return dofsPerCodim[0] * gridView_.size(0) + dofsPerCodim[1] * gridView_.size(1);
}
//! Return number possible values for next position in multi index
size_type size(const SizePrefix prefix) const
{
assert(prefix.size() == 0 || prefix.size() == 1);
return (prefix.size() == 0) ? size() : 0;
}
/** \todo This method has been added to the interface without prior discussion. */
size_type dimension() const
{
return size();
}
size_type maxNodeSize() const
{
return StaticPower<(k+1),GV::dimension>::power;
}
protected:
const GridView gridView_;
std::vector<size_t> codimOffset_;
FiniteElementMap finiteElementMap_;
};
template<typename GV, int k, typename ST, typename TP, GeometryType::BasicType basic_type>
class RaviartThomasNode :
public LeafBasisNode<ST, TP>
{
static const int dim = GV::dimension;
static const int maxSize = StaticPower<(k+1),GV::dimension>::power;
using Base = LeafBasisNode<ST,TP>;
public:
using size_type = ST;
using TreePath = TP;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElementMap = typename Impl::RaviartThomasLocalFiniteElementMap<GV, dim, basic_type, typename GV::ctype, double, k>;
using FiniteElement = typename FiniteElementMap::FiniteElement;
RaviartThomasNode(const TreePath& treePath, const FiniteElementMap* finiteElementMap) :
Base(treePath),
finiteElement_(nullptr),
element_(nullptr),
finiteElementMap_(finiteElementMap)
{ }
//! Return current element, throw if unbound
const Element& element() const
{
return *element_;
}
/** \brief Return the LocalFiniteElement for the element we are bound to
*
* The LocalFiniteElement implements the corresponding interfaces of the dune-localfunctions module
*/
const FiniteElement& finiteElement() const
{
return *finiteElement_;
}
//! Bind to element.
void bind(const Element& e)
{
element_ = &e;
finiteElement_ = &(finiteElementMap_->find(*element_));
this->setSize(finiteElement_->size());
}
protected:
const FiniteElement* finiteElement_;
const Element* element_;
const FiniteElementMap* finiteElementMap_;
};
template<typename GV, int k, class MI, class TP, class ST, GeometryType::BasicType basic_type>
class RaviartThomasNodeIndexSet
{
enum {dim = GV::dimension};
public:
using size_type = ST;
/** \brief Type used for global numbering of the basis vectors */
using MultiIndex = MI;
using NodeFactory = RaviartThomasNodeFactory<GV, k, MI, ST, basic_type>;
using Node = typename NodeFactory::template Node<TP>;
RaviartThomasNodeIndexSet(const NodeFactory& nodeFactory) :
nodeFactory_(&nodeFactory)
{}
/** \brief Bind the view to a grid element
*
* Having to bind the view to an element before being able to actually access any of its data members
* offers to centralize some expensive setup code in the 'bind' method, which can save a lot of run-time.
*/
void bind(const Node& node)
{
node_ = &node;
}
/** \brief Unbind the view
*/
void unbind()
{
node_ = nullptr;
}
/** \brief Size of subtree rooted in this node (element-local)
*/
size_type size() const
{
return node_->finiteElement().size();
}
//! Maps from subtree index set [0..size-1] to a globally unique multi index in global basis
//! Assume dim \in \lbrace 2, 3 \rbrace.
MultiIndex index(size_type i) const
{
Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i);
const auto& gridIndexSet = nodeFactory_->gridView().indexSet();
const auto& element = node_->element();
// The dimension of the entity that the current dof is related to
size_t subentity = localKey.subEntity();
size_t codim = localKey.codim();
// throw if Element is not of predefined type
if (not(basic_type==GeometryType::BasicType::cube and element.type().isCube()) and
not(basic_type==GeometryType::BasicType::simplex and element.type().isSimplex())) DUNE_THROW(Dune::NotImplemented, "RaviartThomasNodalBasis only implemented for cube and simplex elements.");
if (not(codim==0 or codim==1)) DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the RaviartThomasBasis");
return { nodeFactory_->codimOffset_[codim] +
nodeFactory_->dofsPerCodim[codim] * gridIndexSet.subIndex(element, subentity, codim) + localKey.index() };
}
protected:
const NodeFactory* nodeFactory_;
const Node* node_;
};
namespace BasisBuilder {
namespace Imp {
template<std::size_t k, GeometryType::BasicType basic_type>
struct RaviartThomasNodeFactoryBuilder
{
static const std::size_t requiredMultiIndexSize=1;
template<class MultiIndex, class GridView, class size_type=std::size_t>
auto build(const GridView& gridView)
-> RaviartThomasNodeFactory<GridView, k, MultiIndex, size_type, basic_type>
{
return {gridView};
}
};
} // end namespace BasisBuilder::Imp
template<std::size_t k, GeometryType::BasicType basic_type>
Imp::RaviartThomasNodeFactoryBuilder<k, basic_type> rt()
{
return{};
}
} // end namespace BasisBuilder
// *****************************************************************************
// This is the actual global basis implementation based on the reusable parts.
// *****************************************************************************
/** \brief Nodal basis of a scalar k-th-order Raviart Thomas finite element space
*
* TODO
*
* \tparam GV The GridView that the space is defined on
* \tparam k The order of the basis
*/
template<typename GV, int k, GeometryType::BasicType basic_type, class ST = std::size_t>
using RaviartThomasNodalBasis = DefaultGlobalBasis<RaviartThomasNodeFactory<GV, k, FlatMultiIndex<ST>, ST, basic_type> >;
} // end namespace Functions
} // end namespace Dune
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_RAVIARTTHOMASBASIS_HH
......@@ -7,12 +7,16 @@
#include <dune/common/shared_ptr.hh>
#include <dune/localfunctions/finiteelementtypes/functionspacetypes.hh>
#include <dune/functions/functionspacebases/defaultnodetorangemap.hh>
#include <dune/functions/functionspacebases/flatvectorbackend.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/functions/common/treedata.hh>
namespace Dune {
namespace Functions {
......@@ -81,6 +85,13 @@ public:
using Domain = typename EntitySet::GlobalCoordinate;
using Range = R;
// How it should be (not working for Range = double):
//using RangeField = typename Range::value_type;
//using JacobianRange = FieldVector< FieldVector<RangeField, Domain::dimension>, Range::dimension>;
using RangeField = double;
// For vectorial ranges, this works only if Range::dimension == Domain::dimension, otherwise the assert
// in the evaluation method for the derivative throws an error.
using JacobianRange = FieldVector< Range, Domain::dimension>;
using LocalDomain = typename EntitySet::LocalCoordinate;
using Element = typename EntitySet::Element;
......@@ -95,11 +106,16 @@ public:
template<class Node>
using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
template<class Node>
using LocalBasisJacobian = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
template<class Node>
using NodeData = typename std::vector<LocalBasisRange<Node>>;
template<class Node>
using NodeDataJacobian = typename std::vector<LocalBasisJacobian<Node>>;
using ShapeFunctionValueContainer = TreeData<SubTree, NodeData, true>;
using ShapeFunctionJacobianContainer = TreeData<SubTree, NodeDataJacobian, true>;
struct LocalEvaluateVisitor
: public TypeTree::TreeVisitor
......@@ -128,10 +144,41 @@ public:
auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
localBasis.evaluateFunction(x_, shapeFunctionValues);
// Apply function space type dependent continuity preserving transformation.
// Distinguish between H, Hdiv and Hcurl types. Use transformations:
// H: φ → φ
// Hdiv: φ → 1/|det J| J^T φ
// Hcurl: φ → J^(-T) φ
auto&& type = fe.functionSpaceType();
// if type==FunctionSpace::Type::H) no further transformation needed.
if (type==FunctionSpace::Type::Hdiv)
{
auto element = node.element();
auto geometry = element.geometry();
auto jacobianTransposed = geometry.jacobianTransposed(x_);
auto integrationElement = geometry.integrationElement(x_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto tmp = shapeFunctionValues[i];
jacobianTransposed.mtv(tmp, shapeFunctionValues[i]);
shapeFunctionValues[i] /= integrationElement;
}
}
else if (type==FunctionSpace::Type::Hcurl)
{
auto element = node.element();
auto geometry = element.geometry();
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto tmp = shapeFunctionValues[i];
jacobianInverseTransposed.mtv(tmp, shapeFunctionValues[i]);
}
}
// Get range entry associated to this node
auto&& re = nodeToRangeEntry_(node, y_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
......@@ -169,12 +216,158 @@ public:
ShapeFunctionValueContainer& shapeFunctionValueContainer_;
};
// TODO: To be tested properly!
struct LocalJacobianEvaluateVisitor
: public TypeTree::TreeVisitor
, public TypeTree::DynamicTraversal
{
LocalJacobianEvaluateVisitor(const LocalDomain& x, JacobianRange& y, const LocalIndexSet& localIndexSet, const Vector& coefficients, const NodeToRangeEntry& nodeToRangeEntry, ShapeFunctionJacobianContainer& shapeFunctionJacobianContainer):
x_(x),
y_(y),
localIndexSet_(localIndexSet),
coefficients_(coefficients),
nodeToRangeEntry_(nodeToRangeEntry),
shapeFunctionJacobianContainer_(shapeFunctionJacobianContainer)
{}
template<typename Node, typename TreePath>
void leaf(Node& node, TreePath treePath)
{
using LocalBasisJacobianRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
using MultiIndex = typename LocalIndexSet::MultiIndex;
using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
using RangeBlock = typename std::decay<decltype(nodeToRangeEntry_(node, y_))>::type;
auto&& fe = node.finiteElement();
auto&& localBasis = fe.localBasis();
auto&& shapeFunctionJacobians = shapeFunctionJacobianContainer_[node];
localBasis.evaluateJacobian(x_, shapeFunctionJacobians);
// Apply function space type dependent continuity preserving transformation.
// Distinguish between H, Hdiv and Hcurl types. Use transformations:
// H: ∇φ → J^(-T) ∇φ
// Hdiv: ∇φ → 1/|det J| J ∇φ J^(-1)
// Hcurl: ∇φ → ... not implemented
auto&& type = fe.functionSpaceType();
auto element = node.element();
auto geometry = element.geometry();
if (type==FunctionSpace::Type::H)
{
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto referenceShapeFunctionJacobian = shapeFunctionJacobians[i];
jacobianInverseTransposed.mv(referenceShapeFunctionJacobian[0], shapeFunctionJacobians[i][0]);
}
}
else if (type==FunctionSpace::Type::Hdiv)
{
auto integrationElement = geometry.integrationElement(x_);
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x_);
auto jacobianTransposed = geometry.jacobianTransposed(x_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto referenceShapeFunctionJacobian = shapeFunctionJacobians[i];
jacobianTransposed.mtv(referenceShapeFunctionJacobian[0], shapeFunctionJacobians[i][0]);
referenceShapeFunctionJacobian = shapeFunctionJacobians[i];
jacobianInverseTransposed.mtv(referenceShapeFunctionJacobian[0], shapeFunctionJacobians[i][0]);
shapeFunctionJacobians[i] /= integrationElement;
}
}
else if (type==FunctionSpace::Type::Hcurl)
{
DUNE_THROW(Dune::NotImplemented, "Derivatives for curl function not supported by dune-functions.");
}
// Get range entry associated to this node
auto&& re = nodeToRangeEntry_(node, y_);
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
// Get coefficient associated to i-th shape function
auto&& c = coefficients_[multiIndex];
// Get value of i-th shape function
auto&& v = shapeFunctionJacobians[i];
// Notice that the range entry re, the coefficient c, and the shape functions
// value v may all be scalar, vector, matrix, or general container valued.
// The matching of their entries is done via the multistage procedure described
// in the class documentation of DiscreteGlobalBasisFunction.
auto dimC = FlatVectorBackend<CoefficientBlock>::size(c);
auto dimV = FlatVectorBackend<LocalBasisJacobianRange>::size(v);
assert(dimC*dimV == FlatVectorBackend<RangeBlock>::size(re));
for(size_type j=0; j<dimC; ++j)
{
auto&& c_j = FlatVectorBackend<CoefficientBlock>::getEntry(c, j);
for(size_type k=0; k<dimV; ++k)
{
auto&& v_k = FlatVectorBackend<LocalBasisJacobianRange>::getEntry(v, k);
FlatVectorBackend<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
}
}
}
}
const LocalDomain& x_;
JacobianRange& y_;
const LocalIndexSet& localIndexSet_;
const Vector& coefficients_;
const NodeToRangeEntry& nodeToRangeEntry_;
ShapeFunctionJacobianContainer& shapeFunctionJacobianContainer_;
};
struct LocalDofVisitor
: public TypeTree::TreeVisitor
, public TypeTree::DynamicTraversal
{
LocalDofVisitor (std::vector<RangeField>& localDofs, const LocalIndexSet& localIndexSet, const Vector& coefficients) :
localDofs_(localDofs),
localIndexSet_(localIndexSet),
coefficients_(coefficients),
counter(0)
{}
template<typename Node, typename TreePath>
void leaf(Node& node, TreePath treePath)
{
using MultiIndex = typename LocalIndexSet::MultiIndex;
using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
auto&& fe = node.finiteElement();
auto&& localBasis = fe.localBasis();
for (size_type i = 0; i < localBasis.size(); ++i)
{
auto&& multiIndex = localIndexSet_.index(node.localIndex(i));
// Get coefficient associated to i-th shape function
auto&& c = coefficients_[multiIndex];
// Write modus depends on index type (interleafed, ...)
localDofs_[counter++] = c;
auto dimC = FlatVectorBackend<CoefficientBlock>::size(c);
if (dimC != 1 )std::cout << "[Warning] getLocalDofs() does only work for scalar local coefficients. (dune-functions)" << std::endl;
}
}
std::vector<RangeField>& localDofs_;
const LocalIndexSet& localIndexSet_;
const Vector& coefficients_;
size_t counter=0;
};
public:
using GlobalFunction = DiscreteGlobalBasisFunction;
using Domain = LocalDomain;
using Range = GlobalFunction::Range;
using JacobianRange = GlobalFunction::JacobianRange;
using Element = GlobalFunction::Element;
LocalFunction(const DiscreteGlobalBasisFunction& globalFunction)
......@@ -186,6 +379,7 @@ public:
// and queried for tree indices even in unbound state.
subTree_ = &TypeTree::child(localBasisView_.tree(), globalFunction_->treePath());
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
// localDoFs_.reserve(localBasisView_.maxSize());
}
......@@ -198,6 +392,7 @@ public:
// and queried for tree indices even in unbound state.
subTree_ = &TypeTree::child(localBasisView_.tree(), globalFunction_->treePath());
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
}
LocalFunction operator=(const LocalFunction& other)
......@@ -210,6 +405,7 @@ public:
// Here we assume that the tree can be accessed, traversed,
// and queried for tree indices even in unbound state.
shapeFunctionValueContainer_.init(*subTree_);
shapeFunctionJacobianContainer_.init(*subTree_);
}
/**
......@@ -224,8 +420,8 @@ public:
localIndexSet_.bind(localBasisView_);
// Read dofs associated to bound element
// localDoFs_.resize(subTree.size());
// for (size_type i = 0; i < subTree.size(); ++i)
// localDoFs_.resize(subTree_->size());
// for (size_type i = 0; i < subTree_->size(); ++i)
// localDoFs_[i] = globalFunction_->dofs()[localIndexSet_.index(i)];
}
......@@ -254,6 +450,31 @@ public:
return y;
}
/**
* \brief Evaluate derivative of LocalFunction at bound element.
*/
JacobianRange derivative(const Domain& x) const
{
JacobianRange y;
for (size_t j=0; j<y.size(); ++j)
y[j] = 0.0;
LocalJacobianEvaluateVisitor localJacobianEvaluateVisitor(x, y, localIndexSet_, globalFunction_->dofs(), globalFunction_->nodeToRangeEntry(), shapeFunctionJacobianContainer_);
TypeTree::applyToTree(*subTree_, localJacobianEvaluateVisitor);
return y;
}
std::vector<RangeField> localDofs() const
{
std::vector<RangeField> localDofs(subTree_->size());
LocalDofVisitor localDofVisitor(localDofs, localIndexSet_, globalFunction_->dofs());
TypeTree::applyToTree(*subTree_, localDofVisitor);
return localDofs;
}
const Element& localContext() const
{
return localBasisView_.element();
......@@ -271,6 +492,7 @@ public:
LocalIndexSet localIndexSet_;
mutable ShapeFunctionValueContainer shapeFunctionValueContainer_;
mutable ShapeFunctionJacobianContainer shapeFunctionJacobianContainer_;
// std::vector<typename V::value_type> localDoFs_;
const SubTree* subTree_;
};
......
......@@ -6,3 +6,12 @@ target_link_dune_default_libraries("poisson-pq2")
add_executable("stokes-taylorhood" stokes-taylorhood.cc)
target_link_dune_default_libraries("stokes-taylorhood")
add_executable("poisson-mfem" poisson-mfem.cc)
target_link_dune_default_libraries("poisson-mfem")
add_executable("rt0-basistest" rt0-basistest.cc)
target_link_dune_default_libraries("rt0-basistest")
add_executable("derivative-test" derivative-test.cc)
target_link_dune_default_libraries("derivative-test")
#include <config.h>
#include <dune/common/parallel/mpihelper.hh>
#include<dune/geometry/type.hh>
#include<dune/geometry/quadraturerules.hh>
#include <dune/istl/bvector.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/grid/yaspgrid.hh>
#include "dune/functions/functionspacebases/interpolate.hh"
#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <math.h>
using namespace Dune;
template <class GridView, class FunctionA, class FunctionB>
double l2difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionA(quadPos);
diff -= functionB(quadPos);
auto norm = diff.two_norm2();
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
template <class GridView, class FunctionA, class FunctionB>
double h1difference (const GridView& gridView, FunctionA& functionA, FunctionB& functionB)
{
double difference = 0.0;
for (const auto& element: elements(gridView))
{
auto geometry = element.geometry();
functionA.bind(element);
functionB.bind(element);
const auto& quad = QuadratureRules<double, GridView::dimension>::rule(element.type(), 2);
for (const auto& quadPoint : quad)
{
auto quadPos = quadPoint.position();
auto diff = functionB(quadPos);
auto grad = functionA.derivative(quadPos);
for (size_t i=0; i<diff.N(); ++i)
for (size_t j=0; j<diff.M(); ++j)
diff[i][j] -= grad[i][j];
auto norm = diff.frobenius_norm2();
auto integrationElement = geometry.integrationElement(quadPos);
difference += norm * quadPoint.weight() * integrationElement;
}
}
difference = sqrt(difference);
return difference;
}
int main (int argc, char *argv[]) try
{
// Set up MPI, if available
MPIHelper::instance(argc, argv);
//////////////////////////////////////////////////////////////////////////////////////////////
// Generate the grid
//////////////////////////////////////////////////////////////////////////////////////////////
const int dim = 3;
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> ur(1.0);
std::array<int,dim> numElements = {10,10,10};
GridType grid(ur,numElements);
typedef GridType::LeafGridView GridView;
GridView gridView = grid.leafGridView();
for (size_t refineLevel = 0; refineLevel < 3; ++refineLevel)
{
//////////////////////////////////////////////////////////////////////////////////////////////
// Choose a finite element space
//////////////////////////////////////////////////////////////////////////////////////////////
using namespace Functions::BasisBuilder;
auto basis = makeBasis(gridView, power<2>(pq<1>(),flatInterleaved()));
//////////////////////////////////////////////////////////////////////////////////////////////
// Define ranges
//////////////////////////////////////////////////////////////////////////////////////////////
using FunctionRange = FieldVector<double,2>;
using DerivativeRange = FieldMatrix<double,2,dim>;
//////////////////////////////////////////////////////////////////////////////////////////////
// Array type formats
//////////////////////////////////////////////////////////////////////////////////////////////
typedef BlockVector<FieldVector<double,1> > VectorType;
typedef Dune::Functions::HierarchicVectorWrapper<VectorType, double> HierarchicVectorView;
//////////////////////////////////////////////////////////////////////////////////////////////
// Define previous time step and current iterate (pressure, flux), exact solution
//////////////////////////////////////////////////////////////////////////////////////////////
VectorType x;
HierarchicVectorView(x).resize(basis);
x = 0;
//////////////////////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
//////////////////////////////////////////////////////////////////////////////////////////////
auto discreteFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FunctionRange>(basis, TypeTree::hybridTreePath(), HierarchicVectorView(x));
//////////////////////////////////////////////////////////////////////////////////////////////
// Transform discrete functions to local functions and store in a tuple
//////////////////////////////////////////////////////////////////////////////////////////////
auto localDiscreteFunction = localFunction(discreteFunction);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define continuous function and its derivative
//////////////////////////////////////////////////////////////////////////////////////////////
auto pi = std::acos(-1.0);
auto function = [&pi] (const auto& x)
{
FieldVector<double, 2> result(0.0);
result[0] = sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[1] = pow(x[0],2) * exp(x[1]) - x[2];
return result;
};
auto derivative = [&pi] (const auto& x)
{
FieldMatrix<double, 2, 3> result(0.0);
result[0][0] = 2.*pi * cos(2.*pi*x[0]) * sin(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][1] = 2.*pi * sin(2.*pi*x[0]) * cos(2.*pi*x[1]) * sin(2.*pi*x[2]);
result[0][2] = 2.*pi * sin(2.*pi*x[0]) * sin(2.*pi*x[1]) * cos(2.*pi*x[2]);
result[1][0] = 2. * x[0] * exp(x[1]);
result[1][1] = pow(x[0],2) * exp(x[1]);
result[1][2] = - 1.;
return result;
};
auto gridFunction = localFunction(Functions::makeGridViewFunction(function, gridView));
auto gridDerivative = localFunction(Functions::makeGridViewFunction(derivative, gridView));
//////////////////////////////////////////////////////////////////////////////////////////////
// Interpolate function
//////////////////////////////////////////////////////////////////////////////////////////////
interpolate(basis, Dune::TypeTree::hybridTreePath(), HierarchicVectorView(x), function);
//////////////////////////////////////////////////////////////////////////////////////////////
// Define VTK writer and output
//////////////////////////////////////////////////////////////////////////////////////////////
// SubsamplingVTKWriter<GridView> vtkWriter(gridView,1);
// vtkWriter.addVertexData(discreteFunction, VTK::FieldInfo("function", VTK::FieldInfo::Type::vector, 2));
// vtkWriter.write("derivative-test");
//////////////////////////////////////////////////////////////////////////////////////////////
// Check l2 difference for function and its derivative
//////////////////////////////////////////////////////////////////////////////////////////////
auto l2diff = l2difference(gridView, localDiscreteFunction, gridFunction);
auto h1diff = h1difference(gridView, localDiscreteFunction, gridDerivative);
std::cout << "Refine level: " << refineLevel << std::endl;
std::cout << "L2 interpolation error: " << l2diff << std::endl;
std::cout << "H1-seminorm interpolation error: " << h1diff << std::endl << std::endl;
grid.globalRefine(1);
}
}
// Error handling
catch (Exception e) {
std::cout << e << std::endl;
}
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <vector>
#include <cmath>