Go to the documentation of this file.
4 #ifndef DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
5 #define DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
12 #include<dune/common/fvector.hh>
13 #include<dune/common/fmatrix.hh>
14 #include<dune/common/gmpfield.hh>
15 #include<dune/common/exceptions.hh>
16 #include<dune/common/fvector.hh>
18 #include<dune/geometry/referenceelements.hh>
19 #include<dune/geometry/quadraturerules.hh>
20 #include<dune/geometry/type.hh>
22 #include<dune/localfunctions/common/localbasis.hh>
23 #include<dune/localfunctions/common/localkey.hh>
24 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
66 for (
int j=0; j<=k; j++)
82 template<
int k,
int d>
98 alpha[
dim]=k; count++;
100 if (count==i)
return;
103 for (
int m=k-1; m>=0; m--)
108 if (count==i)
return;
119 for (
int j=0; j<d; j++) alpha[j] = 0;
121 for (
int k=1; k<25; k++)
128 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
135 for (
int i=0; i<d; ++i)
144 for (
int j = 0; j<d; ++j) {
145 alpha[j] = i % (k+1);
157 template <BasisType basisType>
163 template <
int k,
int d>
171 template <
int k,
int d>
194 template <
int k,
int d>
202 template <
int k,
int d>
237 for (
long i=k+1; i<=n; i++) nominator *= i;
239 for (
long i=2; i<=n-k; i++) denominator *= i;
240 return nominator/denominator;
245 for (
long i=n-k+1; i<=n; i++) nominator *= i;
247 for (
long i=2; i<=k; i++) denominator *= i;
248 return nominator/denominator;
258 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
265 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
271 template<
typename ComputationFieldType,
int d>
278 ComputationFieldType result(1.0);
279 for (
int i=0; i<d; i++)
281 ComputationFieldType exponent(a[i]+1);
282 result = result/exponent;
290 template<
typename ComputationFieldType>
297 ComputationFieldType one(1.0);
298 ComputationFieldType exponent0(a[0]+1);
299 return one/exponent0;
305 template<
typename ComputationFieldType>
312 ComputationFieldType sum(0.0);
313 for (
int k=0; k<=a[1]+1; k++)
317 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
318 ComputationFieldType denom(a[0]+k+1);
319 sum = sum + (nom/denom);
321 ComputationFieldType denom(a[1]+1);
328 template<
typename ComputationFieldType>
335 ComputationFieldType sumk(0.0);
336 for (
int k=0; k<=a[2]+1; k++)
340 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
341 ComputationFieldType denom(a[1]+k+1);
342 sumk = sumk + (nom/denom);
344 ComputationFieldType sumj(0.0);
345 for (
int j=0; j<=a[1]+a[2]+2; j++)
349 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
350 ComputationFieldType denom(a[0]+j+1);
351 sumj = sumj + (nom/denom);
353 ComputationFieldType denom(a[2]+1);
354 return sumk*sumj/denom;
363 template<
typename F,
int d>
366 template<
typename X,
typename A>
370 for (
int i=0; i<a[d]; i++)
379 template<
typename X,
typename A>
383 for (
int i=0; i<a[0]; i++)
418 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
431 for (
int i=0; i<d; ++i)
434 for (
int i=0; i<
n; i++)
449 for (
int i=0; i<d; ++i)
452 for (
int i=0; i<
n; i++)
465 return Traits::size(l,d);
469 template<
typename Po
int,
typename Result>
472 std::fill(r.begin(),r.end(),0.0);
473 for (
int j=0; j<
n; ++j)
476 for (
int i=j; i<
n; ++i)
477 r[i] += (*coeffs)[i][j] * monomial_value;
482 template<
typename Po
int,
typename Result>
485 std::fill(r.begin(),r.end(),0.0);
487 for (
int j=0; j<
n; ++j)
490 for (
int i=j; i<
n; ++i)
491 for (
int s=0;
s<d; ++
s)
492 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
497 template<
typename Po
int,
typename Result>
501 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
503 for (
int i=0; i<Traits::size(l,d); i++)
506 for (
int j=0; j<=i; j++)
513 template<
typename Po
int,
typename Result>
517 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
519 for (
int i=0; i<Traits::size(l,d); i++)
522 for (
int s=0;
s<d;
s++)
525 for (
int j=0; j<=i; j++)
528 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
534 std::array<std::shared_ptr<MultiIndex<d> >,
n> alpha;
535 std::shared_ptr<LowprecMat> coeffs;
536 std::array<std::shared_ptr<LowprecMat>,d > gradcoeffs;
539 void orthonormalize()
554 for (
int s=0;
s<d;
s++)
555 for (
int i=0; i<
n; i++)
556 for (
int j=0; j<=i; j++)
557 (*gradcoeffs[
s])[i][j] = 0;
558 for (
int i=0; i<
n; i++)
559 for (
int j=0; j<=i; j++)
560 for (
int s=0;
s<d;
s++)
561 if ((*alpha[j])[
s]>0)
564 FieldType factor = beta[
s];
566 int l = invert_index(beta);
567 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
584 int invert_index (MultiIndex<d>& a)
586 for (
int i=0; i<
n; i++)
589 for (
int j=0; j<d; j++)
590 if (a[j]!=(*alpha[i])[j]) found=
false;
593 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
603 for (
int i=0; i<
n; i++)
604 for (
int j=0; j<
n; j++)
606 c[i][j] = ComputationFieldType(1.0);
608 c[i][j] = ComputationFieldType(0.0);
611 MonomialIntegrator<ComputationFieldType,bt,d> integrator;
612 for (
int i=0; i<
n; i++)
615 ComputationFieldType bi[
n];
618 for (
int j=0; j<i; j++)
621 bi[j] = ComputationFieldType(0.0);
622 for (
int l=0; l<=j; l++)
625 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
626 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
628 for (
int l=0; l<=j; l++)
629 c[i][l] = c[i][l] - bi[j]*c[j][l];
633 ComputationFieldType s2(0.0);
635 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
636 s2 = s2 + integrator.integrate(a);
637 for (
int j=0; j<i; j++)
638 s2 = s2 - bi[j]*bi[j];
639 ComputationFieldType
s(1.0);
642 for (
int l=0; l<=i; l++)
647 for (
int i=0; i<
n; i++)
648 for (
int j=0; j<
n; j++)
649 (*coeffs)[i][j] = c[i][j];
661 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
667 Dune::GeometryType gt;
670 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
673 DUNE_NO_DEPRECATED_BEGIN
680 DUNE_NO_DEPRECATED_END
682 unsigned int size ()
const {
return n; }
686 std::vector<typename Traits::RangeType>& out)
const {
694 std::vector<typename Traits::JacobianType>& out)
const {
704 Dune::GeometryType
type ()
const {
return gt; }
707 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
713 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
717 std::size_t
size ()
const {
return n; }
725 std::vector<Dune::LocalKey> li;
737 template<
typename F,
typename C>
741 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
743 typedef typename LB::Traits::RangeType RangeType;
744 const int d = LB::Traits::dimDomain;
745 const Dune::QuadratureRule<RealFieldType,d>&
746 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
750 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
753 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
754 it=rule.begin(); it!=rule.end(); ++it)
757 typename LB::Traits::DomainType x;
759 for (
int i=0; i<d; i++) x[i] = it->position()[i];
763 std::vector<RangeType> phi(LB::n);
764 lb.evaluateFunction(it->position(),phi);
767 for (
int i=0; i<LB::n; i++)
768 out[i] += y*phi[i]*it->weight();
773 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
776 Dune::GeometryType gt;
781 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
785 DUNE_NO_DEPRECATED_BEGIN
788 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
793 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
796 DUNE_NO_DEPRECATED_END
799 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
814 return interpolation;
817 Dune::GeometryType
type ()
const {
return gt; }
828 #endif // DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:720
BasisType
Definition: l2orthonormal.hh:153
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:367
Definition: l2orthonormal.hh:662
compute
Definition: l2orthonormal.hh:364
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:792
@ Pk
Definition: l2orthonormal.hh:154
DUNE_NO_DEPRECATED_END OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:798
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:685
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:428
static const int dim
Definition: adaptivity.hh:84
Dune::GeometryType type() const
Definition: l2orthonormal.hh:817
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:425
double DefaultComputationalFieldType
Definition: l2orthonormal.hh:42
Dune::GeometryType type() const
Definition: l2orthonormal.hh:704
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:819
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:295
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:419
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:263
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:446
Dune::LocalFiniteElementTraits< OPBLocalBasis< D, R, k, d, bt, ComputationFieldType, basisType >, OPBLocalCoefficients< k, d, basisType >, OPBLocalInterpolation< OPBLocalBasis< D, R, k, d, bt, ComputationFieldType, basisType > > > Traits
Definition: l2orthonormal.hh:783
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:693
int size(int l)
Definition: l2orthonormal.hh:463
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:276
DUNE_NO_DEPRECATED_END unsigned int size() const
Definition: l2orthonormal.hh:682
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:738
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:812
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:514
DUNE_NO_DEPRECATED_BEGIN OPBLocalFiniteElement()
Definition: l2orthonormal.hh:787
Definition: l2orthonormal.hh:50
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:424
constexpr int pk_size(int k, int d)
Definition: l2orthonormal.hh:61
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:333
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:228
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:117
MultiIndex()
Definition: l2orthonormal.hh:54
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:259
Definition: l2orthonormal.hh:708
constexpr int qk_size(int k, int d)
Definition: l2orthonormal.hh:132
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:72
static int size(int k, int d)
Definition: l2orthonormal.hh:211
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:117
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:678
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:185
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:807
Definition: l2orthonormal.hh:774
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:717
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:470
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:734
@ Qk
Definition: l2orthonormal.hh:154
Definition: l2orthonormal.hh:158
@ n
Definition: l2orthonormal.hh:671
@ value
Definition: l2orthonormal.hh:86
DUNE_NO_DEPRECATED_BEGIN OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:675
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:498
const std::string s
Definition: function.hh:830
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
@ n
Definition: l2orthonormal.hh:423
Definition: l2orthonormal.hh:83
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:483
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:92
static int size(int k, int d)
Definition: l2orthonormal.hh:179
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:217
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:700
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:712
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:380
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:802
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:142
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:310
Definition: l2orthonormal.hh:729
const P & p
Definition: constraints.hh:147
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: l2orthonormal.hh:670