4 #ifndef DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
5 #define DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
9 #include <dune/localfunctions/common/localbasis.hh>
10 #include <dune/localfunctions/common/localkey.hh>
11 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
12 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
13 #include <dune/localfunctions/common/localinterpolation.hh>
22 template<
class D,
class R,
int k>
33 for (
int order=1; order<=40; order++)
35 const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
38 matched_order = order;
39 matched_size = rule.size();
44 if (matched_order<0) DUNE_THROW(Dune::Exception,
"could not find Gauss Lobatto rule of appropriate size");
45 const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
47 for (
typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
50 size_t member=count%2;
52 if (member==1) newj=group;
else newj=k-group;
53 xi_gl[newj] = it->position()[0];
54 w_gl[newj] = it->weight();
57 for (
int j=0; j<matched_size/2; j++)
61 xi_gl[j] = xi_gl[k-j];
75 R
p (
int i, D
x)
const
78 for (
int j=0; j<=k; j++)
79 if (j!=i) result *= (
x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
84 R
dp (
int i, D
x)
const
88 for (
int j=0; j<=k; j++)
91 R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
92 for (
int l=0; l<=k; l++)
93 if (l!=i && l!=j) prod *= (
x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
124 template<
class D,
class R,
int k,
int d>
131 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
141 std::vector<typename Traits::RangeType>& out)
const
144 for (
size_t i=0; i<
size(); i++)
147 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
153 for (
int j=0; j<d; j++)
154 out[i] *= poly.
p(alpha[j],in[j]);
161 std::vector<typename Traits::JacobianType>& out)
const
166 for (
size_t i=0; i<
size(); i++)
169 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
172 for (
int j=0; j<d; j++)
176 out[i][0][j] = poly.
dp(alpha[j],in[j]);
179 for (
int l=0; l<d; l++)
181 out[i][0][j] *= poly.
p(alpha[l],in[l]);
188 const typename Traits::DomainType& in,
189 std::vector<typename Traits::RangeType>& out)
const {
190 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
191 if (totalOrder == 0) {
194 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
206 template<
int k,
int d,
class LB>
214 template<
typename F,
typename C>
217 typename LB::Traits::DomainType x;
218 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
225 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
228 for (
int j=0; j<d; j++)
229 x[j] = poly.
x(alpha[j]);
237 template<
int d,
class LB>
242 template<
typename F,
typename C>
245 typename LB::Traits::DomainType x(0.5);
246 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
256 template<
class D,
class R,
int k,
int d>
269 typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation>
Traits;
274 : gt(GeometryTypes::cube(d))
295 return interpolation;
319 LocalCoefficients coefficients;
320 LocalInterpolation interpolation;
330 template<
class Geometry,
class RF,
int k>
332 public ScalarLocalToGlobalFiniteElementAdaptorFactory<
333 QkDGGLLocalFiniteElement<
334 typename Geometry::ctype, RF, k, Geometry::mydimension
340 typename Geometry::ctype, RF, k, Geometry::mydimension
342 typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
344 static const LFE lfe;
351 template<
class Geometry,
class RF,
int k>
352 const typename QkDGGLFiniteElementFactory<Geometry, RF, k>::LFE
353 QkDGGLFiniteElementFactory<Geometry, RF, k>::lfe;
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Definition: qkdglagrange.hh:26
Layout map for Q1 elements.
Definition: qkdglagrange.hh:232
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:24
R w(int i) const
Definition: qkdglobatto.hh:106
GaussLobattoLagrangePolynomials()
Definition: qkdglobatto.hh:29
R dp(int i, D x) const
Definition: qkdglobatto.hh:84
R p(int i, D x) const
Definition: qkdglobatto.hh:75
R x(int i) const
Definition: qkdglobatto.hh:100
Lagrange shape functions of order k on the reference cube.
Definition: qkdglobatto.hh:126
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglobatto.hh:199
unsigned int size() const
number of shape functions
Definition: qkdglobatto.hh:134
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglobatto.hh:140
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglobatto.hh:131
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglobatto.hh:160
void partial(const std::array< unsigned int, Traits::dimDomain > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivative of all shape functions.
Definition: qkdglobatto.hh:187
Definition: qkdglobatto.hh:208
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:215
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:243
Definition: qkdglobatto.hh:258
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:269
std::size_t size() const
Definition: qkdglobatto.hh:300
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:286
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:293
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:279
GeometryType type() const
Definition: qkdglobatto.hh:307
@ n
Definition: qkdglobatto.hh:265
QkDGGLLocalFiniteElement * clone() const
Definition: qkdglobatto.hh:312
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:273
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:338
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:348
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139