1 #ifndef DUNE_PM_ORTHONORMAL_HH
2 #define DUNE_PM_ORTHONORMAL_HH
6 #include<dune/common/fvector.hh>
7 #include<dune/common/fmatrix.hh>
8 #include<dune/common/exceptions.hh>
9 #include<dune/common/fvector.hh>
10 #include<dune/common/gmpfield.hh>
11 #include<dune/common/shared_ptr.hh>
12 #include<dune/common/array.hh>
14 #include<dune/geometry/referenceelements.hh>
15 #include<dune/geometry/quadraturerules.hh>
16 #include<dune/geometry/type.hh>
18 #include<dune/localfunctions/common/localbasis.hh>
19 #include<dune/localfunctions/common/localkey.hh>
20 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
39 template<
int k,
int d>
struct PkSize;
41 template<
int j,
int k,
int d>
49 template<
int k,
int d>
60 template<
int k,
int d>
93 template<
int k,
int d>
116 template<
int k,
int d>
167 if (d==1)
return k+1;
169 for (
int j=0; j<=k; j++)
191 alpha[
dim]=k; count++;
193 if (count==i)
return;
196 for (
int m=k-1; m>=0; m--)
201 if (count==i)
return;
212 for (
int j=0; j<d; j++) alpha[j] = 0;
214 for (
int k=1; k<25; k++)
221 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
228 for (
int i=0; i<d; ++i)
237 for (
int j = 0; j<d; ++j) {
238 alpha[j] = i % (k+1);
250 template <BasisType basisType>
256 template <
int k,
int d>
264 template <
int k,
int d>
287 template <
int k,
int d>
295 template <
int k,
int d>
330 for (
long i=k+1; i<=n; i++) nominator *= i;
332 for (
long i=2; i<=n-k; i++) denominator *= i;
333 return nominator/denominator;
338 for (
long i=n-k+1; i<=n; i++) nominator *= i;
340 for (
long i=2; i<=k; i++) denominator *= i;
341 return nominator/denominator;
351 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
358 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
364 template<
typename ComputationFieldType,
int d>
371 ComputationFieldType result(1.0);
372 for (
int i=0; i<d; i++)
374 ComputationFieldType exponent(a[i]+1);
375 result = result/exponent;
383 template<
typename ComputationFieldType>
390 ComputationFieldType one(1.0);
391 ComputationFieldType exponent0(a[0]+1);
392 return one/exponent0;
398 template<
typename ComputationFieldType>
405 ComputationFieldType sum(0.0);
406 for (
int k=0; k<=a[1]+1; k++)
410 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
411 ComputationFieldType denom(a[0]+k+1);
412 sum = sum + (nom/denom);
414 ComputationFieldType denom(a[1]+1);
421 template<
typename ComputationFieldType>
428 ComputationFieldType sumk(0.0);
429 for (
int k=0; k<=a[2]+1; k++)
433 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
434 ComputationFieldType denom(a[1]+k+1);
435 sumk = sumk + (nom/denom);
437 ComputationFieldType sumj(0.0);
438 for (
int j=0; j<=a[1]+a[2]+2; j++)
442 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
443 ComputationFieldType denom(a[0]+j+1);
444 sumj = sumj + (nom/denom);
446 ComputationFieldType denom(a[2]+1);
447 return sumk*sumj/denom;
456 template<
typename F,
int d>
459 template<
typename X,
typename A>
463 for (
int i=0; i<a[d]; i++)
472 template<
typename X,
typename A>
476 for (
int i=0; i<a[0]; i++)
511 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
524 for (
int i=0; i<d; ++i)
527 for (
int i=0; i<n; i++)
542 for (
int i=0; i<d; ++i)
545 for (
int i=0; i<n; i++)
558 return Traits::size(l,d);
562 template<
typename Po
int,
typename Result>
565 std::fill(r.begin(),r.end(),0.0);
566 for (
int j=0; j<n; ++j)
569 for (
int i=j; i<n; ++i)
570 r[i] += (*coeffs)[i][j] * monomial_value;
575 template<
typename Po
int,
typename Result>
578 std::fill(r.begin(),r.end(),0.0);
580 for (
int j=0; j<n; ++j)
583 for (
int i=j; i<n; ++i)
584 for (
int s=0;
s<d; ++
s)
585 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
590 template<
typename Po
int,
typename Result>
594 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
596 for (
int i=0; i<Traits::size(l,d); i++)
599 for (
int j=0; j<=i; j++)
606 template<
typename Po
int,
typename Result>
610 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
612 for (
int i=0; i<Traits::size(l,d); i++)
615 for (
int s=0;
s<d;
s++)
618 for (
int j=0; j<=i; j++)
621 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
627 Dune::array<Dune::shared_ptr<MultiIndex<d> >,n> alpha;
628 Dune::shared_ptr<LowprecMat> coeffs;
629 Dune::array<Dune::shared_ptr<LowprecMat>,d > gradcoeffs;
632 void orthonormalize()
647 for (
int s=0;
s<d;
s++)
648 for (
int i=0; i<n; i++)
649 for (
int j=0; j<=i; j++)
650 (*gradcoeffs[
s])[i][j] = 0;
651 for (
int i=0; i<n; i++)
652 for (
int j=0; j<=i; j++)
653 for (
int s=0; s<d; s++)
654 if ((*alpha[j])[s]>0)
657 FieldType factor = beta[
s];
659 int l = invert_index(beta);
660 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
677 int invert_index (MultiIndex<d>& a)
679 for (
int i=0; i<n; i++)
682 for (
int j=0; j<d; j++)
683 if (a[j]!=(*alpha[i])[j]) found=
false;
686 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
692 HighprecMat *
p =
new HighprecMat();
696 for (
int i=0; i<n; i++)
697 for (
int j=0; j<n; j++)
699 c[i][j] = ComputationFieldType(1.0);
701 c[i][j] = ComputationFieldType(0.0);
704 MonomialIntegrator<ComputationFieldType,bt,d> integrator;
705 for (
int i=0; i<n; i++)
708 ComputationFieldType bi[n];
711 for (
int j=0; j<i; j++)
714 bi[j] = ComputationFieldType(0.0);
715 for (
int l=0; l<=j; l++)
718 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
719 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
721 for (
int l=0; l<=j; l++)
722 c[i][l] = c[i][l] - bi[j]*c[j][l];
726 ComputationFieldType s2(0.0);
728 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
729 s2 = s2 + integrator.integrate(a);
730 for (
int j=0; j<i; j++)
731 s2 = s2 - bi[j]*bi[j];
732 ComputationFieldType
s(1.0);
734 for (
int l=0; l<=i; l++)
739 for (
int i=0; i<n; i++)
740 for (
int j=0; j<n; j++)
741 (*coeffs)[i][j] = c[i][j];
753 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=R, PB::BasisType basisType = PB::BasisType::Pk>
759 Dune::GeometryType gt;
762 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 0>
Traits;
770 unsigned int size ()
const {
return n; }
774 std::vector<typename Traits::RangeType>& out)
const {
776 opb.evaluateFunction(in,out);
782 std::vector<typename Traits::JacobianType>& out)
const {
784 opb.evaluateJacobian(in,out);
792 Dune::GeometryType
type ()
const {
return gt; }
795 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
801 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
805 std::size_t
size ()
const {
return n; }
813 std::vector<Dune::LocalKey> li;
825 template<
typename F,
typename C>
829 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
831 typedef typename LB::Traits::RangeType RangeType;
832 const int d = LB::Traits::dimDomain;
833 const Dune::QuadratureRule<RealFieldType,d>&
834 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
838 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
841 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
842 it=rule.begin(); it!=rule.end(); ++it)
845 typename LB::Traits::DomainType x;
847 for (
int i=0; i<d; i++) x[i] = it->position()[i];
851 std::vector<RangeType> phi(LB::n);
852 lb.evaluateFunction(it->position(),phi);
855 for (
int i=0; i<LB::n; i++)
856 out[i] += y*phi[i]*it->weight();
861 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=R, PB::BasisType basisType = PB::BasisType::Pk>
864 Dune::GeometryType gt;
869 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
874 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
879 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
883 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
898 return interpolation;
901 Dune::GeometryType
type ()
const {
return gt; }
Definition: l2orthonormal.hh:64
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:871
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:426
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:539
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:388
Definition: l2orthonormal.hh:94
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:805
Definition: l2orthonormal.hh:251
OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:765
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:460
Definition: l2orthonormal.hh:120
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:576
Definition: l2orthonormal.hh:247
Dune::GeometryType type() const
Definition: l2orthonormal.hh:901
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:896
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:886
Definition: l2orthonormal.hh:39
OPBLocalFiniteElement()
Definition: l2orthonormal.hh:873
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:512
OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:882
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:607
Definition: l2orthonormal.hh:45
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:310
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:891
Definition: l2orthonormal.hh:862
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:563
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:773
unsigned int size() const
Definition: l2orthonormal.hh:770
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:175
static const int dim
Definition: adaptivity.hh:82
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:822
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:903
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:473
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:768
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:788
static int size(int k, int d)
Definition: l2orthonormal.hh:304
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:517
Definition: l2orthonormal.hh:153
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:369
int size(int l)
Definition: l2orthonormal.hh:556
Dune::GeometryType type() const
Definition: l2orthonormal.hh:792
static int size(int k, int d)
Definition: l2orthonormal.hh:272
Definition: l2orthonormal.hh:817
Definition: l2orthonormal.hh:42
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:800
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:210
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:521
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:518
Definition: l2orthonormal.hh:796
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdg.hh:125
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: l2orthonormal.hh:247
const F & f
Definition: common/constraints.hh:145
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:403
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:185
MultiIndex()
Definition: l2orthonormal.hh:157
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:278
int qk_size(int k, int d)
Definition: l2orthonormal.hh:225
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:878
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, 0 > Traits
Definition: l2orthonormal.hh:762
BasisType
Definition: l2orthonormal.hh:246
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:321
Definition: l2orthonormal.hh:754
Definition: l2orthonormal.hh:97
R p(int i, D x)
Definition: qkdg.hh:62
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:356
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:235
int pk_size(int k, int d)
Definition: l2orthonormal.hh:164
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:808
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:781
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:352
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:591
Definition: l2orthonormal.hh:117
compute
Definition: l2orthonormal.hh:457
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:826
const std::string s
Definition: function.hh:1103