4#ifndef DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
5#define DUNE_LOCALFUNCTIONS_QKLOCALBASIS_HH
9#include <dune/common/fvector.hh>
10#include <dune/common/fmatrix.hh>
11#include <dune/common/power.hh>
13#include <dune/geometry/type.hh>
33 template<
class D,
class R,
int k,
int d>
36 enum { n = StaticPower<k+1,d>::power };
39 static R p (
int i, D x)
42 for (
int j=0; j<=k; j++)
43 if (j!=i) result *= (k*x-j)/(i-j);
48 static R dp (
int i, D x)
52 for (
int j=0; j<=k; j++)
55 R prod( (k*1.0)/(i-j) );
56 for (
int l=0; l<=k; l++)
58 prod *= (k*x-l)/(i-l);
65 static Dune::FieldVector<int,d> multiindex (
int i)
67 Dune::FieldVector<int,d> alpha;
68 for (
int j=0; j<d; j++)
77 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 1>
Traits;
82 return StaticPower<k+1,d>::power;
87 std::vector<typename Traits::RangeType>& out)
const
90 for (
size_t i=0; i<
size(); i++)
93 Dune::FieldVector<int,d> alpha(multiindex(i));
99 for (
int j=0; j<d; j++)
100 out[i] *= p(alpha[j],in[j]);
110 std::vector<typename Traits::JacobianType>& out)
const
115 for (
size_t i=0; i<
size(); i++)
118 Dune::FieldVector<int,d> alpha(multiindex(i));
121 for (
int j=0; j<d; j++)
125 out[i][0][j] = dp(alpha[j],in[j]);
128 for (
int l=0; l<d; l++)
130 out[i][0][j] *= p(alpha[l],in[l]);
142 std::vector<typename Traits::RangeType>& out)
const
144 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
156 for (
size_t i=0; i<
size(); i++)
159 Dune::FieldVector<int,d> alpha(multiindex(i));
165 for (std::size_t l=0; l<d; l++)
166 out[i][0] *= (
order[l]) ? dp(alpha[l],in[l]) : p(alpha[l],in[l]);
171 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
180 template<
int diffOrder>
182 const std::array<int,1>& direction,
184 std::vector<typename Traits::RangeType>& out)
const
186 static_assert(diffOrder == 1,
"We only can compute first derivatives");
190 for (
size_t i=0; i<
size(); i++)
193 Dune::FieldVector<int,d> alpha(multiindex(i));
196 std::size_t j = direction[0];
200 out[i][0] = dp(alpha[j],in[j]);
203 for (std::size_t l=0; l<d; l++)
205 out[i][0] *= p(alpha[l],in[l]);
Definition brezzidouglasmarini1cube2d.hh:14
Type traits for LocalBasisVirtualInterface.
Definition localbasis.hh:38
D DomainType
domain type
Definition localbasis.hh:49
Lagrange shape functions of order k on the reference cube.
Definition qklocalbasis.hh:35
void partial(const std::array< unsigned int, d > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition qklocalbasis.hh:140
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition qklocalbasis.hh:86
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, 1 > Traits
Definition qklocalbasis.hh:77
unsigned int size() const
number of shape functions
Definition qklocalbasis.hh:80
void evaluate(const std::array< int, 1 > &direction, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate derivative in a given direction.
Definition qklocalbasis.hh:181
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition qklocalbasis.hh:109
unsigned int order() const
Polynomial order of the shape functions.
Definition qklocalbasis.hh:210