1 #ifndef DUNE_PM_ORTHONORMAL_HH
2 #define DUNE_PM_ORTHONORMAL_HH
8 #include<dune/common/fvector.hh>
9 #include<dune/common/fmatrix.hh>
10 #include<dune/common/gmpfield.hh>
11 #include<dune/common/exceptions.hh>
12 #include<dune/common/fvector.hh>
13 #include<dune/common/array.hh>
15 #include<dune/geometry/referenceelements.hh>
16 #include<dune/geometry/quadraturerules.hh>
17 #include<dune/geometry/type.hh>
19 #include<dune/localfunctions/common/localbasis.hh>
20 #include<dune/localfunctions/common/localkey.hh>
21 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
46 template<
int k,
int d>
struct PkSize;
48 template<
int j,
int k,
int d>
56 template<
int k,
int d>
67 template<
int k,
int d>
100 template<
int k,
int d>
123 template<
int k,
int d>
174 if (d==1)
return k+1;
176 for (
int j=0; j<=k; j++)
198 alpha[
dim]=k; count++;
200 if (count==i)
return;
203 for (
int m=k-1; m>=0; m--)
208 if (count==i)
return;
219 for (
int j=0; j<d; j++) alpha[j] = 0;
221 for (
int k=1; k<25; k++)
228 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
235 for (
int i=0; i<d; ++i)
244 for (
int j = 0; j<d; ++j) {
245 alpha[j] = i % (k+1);
257 template <BasisType basisType>
263 template <
int k,
int d>
271 template <
int k,
int d>
294 template <
int k,
int d>
302 template <
int k,
int d>
337 for (
long i=k+1; i<=n; i++) nominator *= i;
339 for (
long i=2; i<=n-k; i++) denominator *= i;
340 return nominator/denominator;
345 for (
long i=n-k+1; i<=n; i++) nominator *= i;
347 for (
long i=2; i<=k; i++) denominator *= i;
348 return nominator/denominator;
358 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
365 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
371 template<
typename ComputationFieldType,
int d>
378 ComputationFieldType result(1.0);
379 for (
int i=0; i<d; i++)
381 ComputationFieldType exponent(a[i]+1);
382 result = result/exponent;
390 template<
typename ComputationFieldType>
397 ComputationFieldType one(1.0);
398 ComputationFieldType exponent0(a[0]+1);
399 return one/exponent0;
405 template<
typename ComputationFieldType>
412 ComputationFieldType sum(0.0);
413 for (
int k=0; k<=a[1]+1; k++)
417 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
418 ComputationFieldType denom(a[0]+k+1);
419 sum = sum + (nom/denom);
421 ComputationFieldType denom(a[1]+1);
428 template<
typename ComputationFieldType>
435 ComputationFieldType sumk(0.0);
436 for (
int k=0; k<=a[2]+1; k++)
440 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
441 ComputationFieldType denom(a[1]+k+1);
442 sumk = sumk + (nom/denom);
444 ComputationFieldType sumj(0.0);
445 for (
int j=0; j<=a[1]+a[2]+2; j++)
449 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
450 ComputationFieldType denom(a[0]+j+1);
451 sumj = sumj + (nom/denom);
453 ComputationFieldType denom(a[2]+1);
454 return sumk*sumj/denom;
463 template<
typename F,
int d>
466 template<
typename X,
typename A>
470 for (
int i=0; i<a[d]; i++)
479 template<
typename X,
typename A>
483 for (
int i=0; i<a[0]; i++)
518 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
529 : coeffs(new LowprecMat)
531 for (
int i=0; i<d; ++i)
532 gradcoeffs[i].reset(
new LowprecMat());
534 for (
int i=0; i<n; i++)
547 : coeffs(new LowprecMat)
549 for (
int i=0; i<d; ++i)
550 gradcoeffs[i].reset(
new LowprecMat());
552 for (
int i=0; i<n; i++)
565 return Traits::size(l,d);
569 template<
typename Po
int,
typename Result>
572 std::fill(r.begin(),r.end(),0.0);
573 for (
int j=0; j<n; ++j)
576 for (
int i=j; i<n; ++i)
577 r[i] += (*coeffs)[i][j] * monomial_value;
582 template<
typename Po
int,
typename Result>
585 std::fill(r.begin(),r.end(),0.0);
587 for (
int j=0; j<n; ++j)
590 for (
int i=j; i<n; ++i)
591 for (
int s=0;
s<d; ++
s)
592 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
597 template<
typename Po
int,
typename Result>
601 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
603 for (
int i=0; i<Traits::size(l,d); i++)
606 for (
int j=0; j<=i; j++)
613 template<
typename Po
int,
typename Result>
617 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
619 for (
int i=0; i<Traits::size(l,d); i++)
622 for (
int s=0;
s<d;
s++)
625 for (
int j=0; j<=i; j++)
628 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
634 Dune::array<std::shared_ptr<MultiIndex<d> >,n> alpha;
635 std::shared_ptr<LowprecMat> coeffs;
636 Dune::array<std::shared_ptr<LowprecMat>,d > gradcoeffs;
639 void orthonormalize()
654 for (
int s=0;
s<d;
s++)
655 for (
int i=0; i<n; i++)
656 for (
int j=0; j<=i; j++)
657 (*gradcoeffs[
s])[i][j] = 0;
658 for (
int i=0; i<n; i++)
659 for (
int j=0; j<=i; j++)
660 for (
int s=0; s<d; s++)
661 if ((*alpha[j])[s]>0)
664 FieldType factor = beta[
s];
666 int l = invert_index(beta);
667 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
684 int invert_index (MultiIndex<d>& a)
686 for (
int i=0; i<n; i++)
689 for (
int j=0; j<d; j++)
690 if (a[j]!=(*alpha[i])[j]) found=
false;
693 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
699 HighprecMat *
p =
new HighprecMat();
703 for (
int i=0; i<n; i++)
704 for (
int j=0; j<n; j++)
706 c[i][j] = ComputationFieldType(1.0);
708 c[i][j] = ComputationFieldType(0.0);
711 MonomialIntegrator<ComputationFieldType,bt,d> integrator;
712 for (
int i=0; i<n; i++)
715 ComputationFieldType bi[n];
718 for (
int j=0; j<i; j++)
721 bi[j] = ComputationFieldType(0.0);
722 for (
int l=0; l<=j; l++)
725 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
726 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
728 for (
int l=0; l<=j; l++)
729 c[i][l] = c[i][l] - bi[j]*c[j][l];
733 ComputationFieldType s2(0.0);
735 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
736 s2 = s2 + integrator.integrate(a);
737 for (
int j=0; j<i; j++)
738 s2 = s2 - bi[j]*bi[j];
739 ComputationFieldType
s(1.0);
742 for (
int l=0; l<=i; l++)
747 for (
int i=0; i<n; i++)
748 for (
int j=0; j<n; j++)
749 (*coeffs)[i][j] = c[i][j];
761 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
767 Dune::GeometryType gt;
770 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 0>
Traits;
778 unsigned int size ()
const {
return n; }
782 std::vector<typename Traits::RangeType>& out)
const {
790 std::vector<typename Traits::JacobianType>& out)
const {
800 Dune::GeometryType
type ()
const {
return gt; }
803 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
809 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
813 std::size_t
size ()
const {
return n; }
821 std::vector<Dune::LocalKey> li;
833 template<
typename F,
typename C>
837 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
839 typedef typename LB::Traits::RangeType RangeType;
840 const int d = LB::Traits::dimDomain;
841 const Dune::QuadratureRule<RealFieldType,d>&
842 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
846 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
849 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
850 it=rule.begin(); it!=rule.end(); ++it)
853 typename LB::Traits::DomainType x;
855 for (
int i=0; i<d; i++) x[i] = it->position()[i];
859 std::vector<RangeType> phi(LB::n);
860 lb.evaluateFunction(it->position(),phi);
863 for (
int i=0; i<LB::n; i++)
864 out[i] += y*phi[i]*it->weight();
869 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
872 Dune::GeometryType gt;
877 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
882 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
887 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
891 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
906 return interpolation;
909 Dune::GeometryType
type ()
const {
return gt; }
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:776
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:904
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:525
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:376
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:813
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:217
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:410
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:359
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:192
Definition: l2orthonormal.hh:258
Definition: l2orthonormal.hh:52
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:363
Definition: l2orthonormal.hh:49
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:524
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:528
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:879
Definition: l2orthonormal.hh:71
Definition: l2orthonormal.hh:254
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:546
int pk_size(int k, int d)
Definition: l2orthonormal.hh:171
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:911
OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:890
Definition: l2orthonormal.hh:762
compute
Definition: l2orthonormal.hh:464
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:328
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:834
Dune::GeometryType type() const
Definition: l2orthonormal.hh:800
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:796
int qk_size(int k, int d)
Definition: l2orthonormal.hh:232
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:317
MultiIndex()
Definition: l2orthonormal.hh:164
unsigned int size() const
Definition: l2orthonormal.hh:778
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:583
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:395
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:598
Dune::GeometryType type() const
Definition: l2orthonormal.hh:909
Definition: l2orthonormal.hh:160
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:614
Definition: l2orthonormal.hh:254
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:894
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:480
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:886
double DefaultComputationalFieldType
Definition: l2orthonormal.hh:39
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:519
int size(int l)
Definition: l2orthonormal.hh:563
static int size(int k, int d)
Definition: l2orthonormal.hh:311
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:433
Definition: l2orthonormal.hh:104
Definition: adaptivity.hh:27
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdg.hh:125
OPBLocalFiniteElement()
Definition: l2orthonormal.hh:881
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:816
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:182
static const int dim
Definition: adaptivity.hh:83
Definition: l2orthonormal.hh:101
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:242
Definition: l2orthonormal.hh:124
const std::string s
Definition: function.hh:1102
Definition: l2orthonormal.hh:825
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:899
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:789
Definition: l2orthonormal.hh:870
const P & p
Definition: constraints.hh:146
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:830
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:467
BasisType
Definition: l2orthonormal.hh:253
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:781
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:808
Definition: l2orthonormal.hh:127
OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:773
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:285
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, 0 > Traits
Definition: l2orthonormal.hh:770
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:570
Definition: l2orthonormal.hh:46
static int size(int k, int d)
Definition: l2orthonormal.hh:279
Definition: l2orthonormal.hh:804