Skip to content
Snippets Groups Projects
Commit 4eeedf35 authored by Porrmann, Maik's avatar Porrmann, Maik
Browse files

Use GlobalValuedFE for dim == 1 too 

Fix Interpolation, assuming global derivating
parent d9c58c32
No related branches found
Tags Hermite_1d_2d_interpolationConvergence
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/functions/functionspacebases/hermitebasis.hh
+ 96
− 39
View file @ 4eeedf35
......@@ -7,10 +7,6 @@
#include <dune/common/exceptions.hh>
#include <dune/localfunctions/hermite.hh>
//needed ?? TODO remove
// #include <dune/localfunctions/lagrange/equidistantpoints.hh>
// #include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/flatmultiindex.hh>
......@@ -21,43 +17,108 @@ namespace Dune
{
namespace Functions
{
namespace Impl{
namespace Impl
{
/**
* \brief Implements a transformation from pullbacks of reference basis functions to global basisfunctions for Hermite elements.
* For more information, see globalvaluedlocalfiniteelement.hh
*
*/
struct HermiteTransformator{
template <unsigned int dim>
struct HermiteTransformator{};
//#Todo write this dimension independent
template<typename Values, typename LocalCoordinate, typename Geometry>
static auto apply(Values& values, const LocalCoordinate& xi, const Geometry& geometry)
template<>
struct HermiteTransformator<1>
{
template <typename Values, typename LocalCoordinate, typename Geometry>
static auto apply(Values &values, const LocalCoordinate &xi, const Geometry &geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
// transform the values of basisfunctions that correspond to a "derivative" node
for (std::size_t i = 0; i < 2; ++i)
{
Dune::FieldVector<double, 1> localGrad = {values[i * 2 + 1]};
Dune::FieldVector<double, 1> globalGrad(0.);
jacobianTransposed.mtv(localGrad, globalGrad);
// TODO make this nice by removing globalGrad and write directly into values
values[i * 2 + 1] = globalGrad[0];
}
}
template <typename Gradients, typename LocalCoordinate, typename Geometry>
static auto applyJacobian(Gradients &values, const LocalCoordinate &xi, const Geometry geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
// transform the values of basisfunctions that correspond to a "derivative" node
for (std::size_t i = 0; i<3; ++i){
Dune::FieldVector<double,2> localGrad = {values[i*3+1],values[i*3+2]};
Dune::FieldVector<double,2> globalGrad(0.);
jacobianTransposed.mtv(localGrad,globalGrad);
// #Todo make this nice by removing globalGrad and write directly into values
for (std::size_t i = 0; i < 2; ++i)
{
Dune::FieldVector<double, 1> globalDerivative (0.);
jacobianTransposed.mtv(Dune::FieldVector<double,1>{values[i * 2 + 1][0][0]}, globalDerivative);
// TODO make this nice by removing globalValued and write directly into values
values[i * 2 + 1][0][0] = globalDerivative[0];
}
}
//TODO think about partial interface!
template <class Function, class LocalCoordinate, class Element>
class LocalValuedFunction
{
const Function &f_;
const Element &element_;
public:
LocalValuedFunction(const Function &f, const Element &e)
: f_(f), element_(e) {}
auto operator()(const LocalCoordinate &xi) const
{
auto &&f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
return f(xi); // no transformation for functionvalues needed
}
auto derivative(const LocalCoordinate &xi) const
{
// wrapping this with makeFunctionWithCallOperator causes a loss of the geometry when used with a localFunction (e.g. from Dune::Functions::interpolate)
auto &&df = Dune::Impl::makeDerivative<LocalCoordinate>(f_);
return df(xi);;
}
};
};
template<>
struct HermiteTransformator<2>
{
template <typename Values, typename LocalCoordinate, typename Geometry>
static auto apply(Values &values, const LocalCoordinate &xi, const Geometry &geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
// transform the values of basisfunctions that correspond to a "derivative" node
for (std::size_t i = 0; i < 3; ++i)
{
Dune::FieldVector<double, 2> localGrad = {values[i * 3 + 1], values[i * 3 + 2]};
Dune::FieldVector<double, 2> globalGrad(0.);
jacobianTransposed.mtv(localGrad, globalGrad);
// TODO make this nice by removing globalGrad and write directly into values
values[i * 3 + 1] = globalGrad[0];
values[i * 3 + 2] = globalGrad[1];
}
}
template<typename Gradients, typename LocalCoordinate, typename Geometry>
static auto applyJacobian(Gradients& values, const LocalCoordinate& xi, const Geometry geometry)
template <typename Gradients, typename LocalCoordinate, typename Geometry>
static auto applyJacobian(Gradients &values, const LocalCoordinate &xi, const Geometry geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
// transform the values of basisfunctions that correspond to a "derivative" node
for (std::size_t i = 0; i < 3; ++i)
{
Dune::FieldVector<double, 2> localValued_dx = {values[i * 3 + 1][0][0], values[i * 3 + 2][0][0]}
,localValued_dy = {values[i * 3 + 1][0][1], values[i * 3 + 2][0][1]};
Dune::FieldVector<double, 2> localValued_dx = {values[i * 3 + 1][0][0], values[i * 3 + 2][0][0]}, localValued_dy = {values[i * 3 + 1][0][1], values[i * 3 + 2][0][1]};
auto globalValued_dx = localValued_dx, globalValued_dy = localValued_dy;
jacobianTransposed.mtv(localValued_dx,globalValued_dx);
jacobianTransposed.mtv(localValued_dx, globalValued_dx);
jacobianTransposed.mtv(localValued_dy, globalValued_dy);
// TODO make this nice by removing globalValued and write directly into values
......@@ -73,20 +134,18 @@ namespace Dune
template <class Function, class LocalCoordinate, class Element>
class LocalValuedFunction
{
const Function& f_;
const Element& element_;
const Function &f_;
const Element &element_;
public:
LocalValuedFunction(const Function& f, const Element& e)
: f_(f),element_(e){}
LocalValuedFunction(const Function &f, const Element &e)
: f_(f), element_(e) {}
auto operator()(const LocalCoordinate& xi) const
auto operator()(const LocalCoordinate &xi) const
{
auto&& f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
return f(xi); // no transformation for functionvalues needed
auto &&f = Dune::Impl::makeFunctionWithCallOperator<LocalCoordinate>(f_);
return f(xi); // no transformation for functionvalues needed
}
auto derivative(const LocalCoordinate &xi) const
{
......@@ -102,7 +161,6 @@ namespace Dune
return df(xi);
;
}
};
};
......@@ -111,14 +169,15 @@ namespace Dune
template <typename GV, typename R>
class HermiteNode : public LeafBasisNode
{
static constexpr int dim = GV::dimension;
using LocalValuedFE = HermiteLocalFiniteElement<typename GV::ctype, R, GV::dimension>;
static constexpr unsigned int dim = GV::dimension;
using LocalValuedFE = HermiteLocalFiniteElement<typename GV::ctype, R, dim>;
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = typename std::conditional<(dim == 1), LocalValuedFE, Impl::GlobalValuedLocalFiniteElement<Impl::HermiteTransformator, LocalValuedFE, Element> >::type;
using FiniteElement = //typename std::conditional<(dim == 1), LocalValuedFE,
Impl::GlobalValuedLocalFiniteElement<Impl::HermiteTransformator<dim>, LocalValuedFE, Element> ;//>::type;
HermiteNode() : localValuedFiniteElement_(std::make_shared<LocalValuedFE>()),
finiteElement_(std::make_shared<FiniteElement>()),
......@@ -152,7 +211,7 @@ namespace Dune
if (not e.type().isSimplex())
DUNE_THROW(Dune::NotImplemented, "HermiteBasis can only be bound to simplex elements");
element_ = &e;
if constexpr (dim != 1)
// if constexpr (dim != 1)
finiteElement_->bind(*localValuedFiniteElement_, *element_);
}
......@@ -196,7 +255,7 @@ namespace Dune
//! Constructor for a given grid view object
HermitePreBasis(const GV &gv) : gridView_(gv)
{
if (dim >2)
if (dim > 2)
DUNE_THROW(Dune::NotImplemented, "HermitePreBasis only implemented for dim <= 2");
}
......@@ -276,7 +335,7 @@ namespace Dune
*it = {{2 * ((size_type)(gridIndexSet.subIndex(element, localKey.subEntity(), 1))) + localKey.index()}};
continue;
}
else if(dim == 2)
else if (dim == 2)
{
// TODO think again if you can order this elementwise to keep blockstructure in matrices!
*it = {{(size_type)((localKey.codim() == 2) ? (3 * gridIndexSet.subIndex(element, localKey.subEntity(), 2) + localKey.index()) : 3 * gridIndexSet.size(2) + gridIndexSet.subIndex(element, localKey.subEntity(), 0))}};
......@@ -292,7 +351,6 @@ namespace Dune
}; // class HermitePreBasis
namespace BasisFactory
{
......@@ -300,15 +358,14 @@ namespace Dune
auto hermite()
{
return [](const auto &gridView)
{
return HermitePreBasis<std::decay_t<decltype(gridView)>, R>(gridView);
};
}
} //namespace BasisFactory
} //Functions
} // Dune
#endif
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