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>
21 template<
class D,
class R,
int k>
32 for (
int order=1; order<=40; order++)
34 const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
37 matched_order = order;
38 matched_size = rule.size();
43 if (matched_order<0) DUNE_THROW(Dune::Exception,
"could not find Gauss Lobatto rule of appropriate size");
44 const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
46 for (
typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
49 size_t member=count%2;
51 if (member==1) newj=group;
else newj=k-group;
52 xi_gl[newj] = it->position()[0];
53 w_gl[newj] = it->weight();
56 for (
int j=0; j<matched_size/2; j++)
60 xi_gl[j] = xi_gl[k-j];
74 R
p (
int i, D
x)
const
77 for (
int j=0; j<=k; j++)
78 if (j!=i) result *= (
x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
83 R
dp (
int i, D
x)
const
87 for (
int j=0; j<=k; j++)
90 R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
91 for (
int l=0; l<=k; l++)
92 if (l!=i && l!=j) prod *= (
x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
123 template<
class D,
class R,
int k,
int d>
130 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
140 std::vector<typename Traits::RangeType>& out)
const
143 for (
size_t i=0; i<
size(); i++)
146 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
152 for (
int j=0; j<d; j++)
153 out[i] *= poly.
p(alpha[j],in[j]);
160 std::vector<typename Traits::JacobianType>& out)
const
165 for (
size_t i=0; i<
size(); i++)
168 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
171 for (
int j=0; j<d; j++)
175 out[i][0][j] = poly.
dp(alpha[j],in[j]);
178 for (
int l=0; l<d; l++)
180 out[i][0][j] *= poly.
p(alpha[l],in[l]);
186 void partial(
const std::array<unsigned int, Traits::dimDomain>& DUNE_UNUSED(
order),
187 const typename Traits::DomainType& DUNE_UNUSED(in),
188 std::vector<typename Traits::RangeType>& DUNE_UNUSED(out))
const {
189 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
190 if (totalOrder == 0) {
193 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
205 template<
int k,
int d,
class LB>
213 template<
typename F,
typename C>
216 typename LB::Traits::DomainType x;
223 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
226 for (
int j=0; j<d; j++)
227 x[j] = poly.
x(alpha[j]);
235 template<
int d,
class LB>
240 template<
typename F,
typename C>
243 typename LB::Traits::DomainType x(0.5);
253 template<
class D,
class R,
int k,
int d>
266 typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation>
Traits;
271 : gt(GeometryTypes::cube(d))
292 return interpolation;
316 LocalCoefficients coefficients;
317 LocalInterpolation interpolation;
327 template<
class Geometry,
class RF,
int k>
329 public ScalarLocalToGlobalFiniteElementAdaptorFactory<
330 QkDGGLLocalFiniteElement<
331 typename Geometry::ctype, RF, k, Geometry::mydimension
337 typename Geometry::ctype, RF, k, Geometry::mydimension
339 typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
341 static const LFE lfe;
348 template<
class Geometry,
class RF,
int k>
349 const typename QkDGGLFiniteElementFactory<Geometry, RF, k>::LFE
350 QkDGGLFiniteElementFactory<Geometry, RF, k>::lfe;
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Definition: qkdglagrange.hh:24
Layout map for Q1 elements.
Definition: qkdglagrange.hh:230
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:23
R w(int i) const
Definition: qkdglobatto.hh:105
GaussLobattoLagrangePolynomials()
Definition: qkdglobatto.hh:28
R dp(int i, D x) const
Definition: qkdglobatto.hh:83
R p(int i, D x) const
Definition: qkdglobatto.hh:74
R x(int i) const
Definition: qkdglobatto.hh:99
Lagrange shape functions of order k on the reference cube.
Definition: qkdglobatto.hh:125
void partial(const std::array< unsigned int, Traits::dimDomain > &DUNE_UNUSED(order), const typename Traits::DomainType &DUNE_UNUSED(in), std::vector< typename Traits::RangeType > &DUNE_UNUSED(out)) const
Evaluate partial derivative of all shape functions.
Definition: qkdglobatto.hh:186
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglobatto.hh:198
unsigned int size() const
number of shape functions
Definition: qkdglobatto.hh:133
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglobatto.hh:139
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglobatto.hh:130
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglobatto.hh:159
Definition: qkdglobatto.hh:207
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:214
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:241
Definition: qkdglobatto.hh:255
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:266
std::size_t size() const
Definition: qkdglobatto.hh:297
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:283
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:290
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:276
GeometryType type() const
Definition: qkdglobatto.hh:304
@ n
Definition: qkdglobatto.hh:262
QkDGGLLocalFiniteElement * clone() const
Definition: qkdglobatto.hh:309
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:270
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:335
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:345
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139