4 #ifndef DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH 5 #define DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH 11 #include<dune/common/fvector.hh> 12 #include<dune/common/fmatrix.hh> 13 #include<dune/common/gmpfield.hh> 14 #include<dune/common/exceptions.hh> 15 #include<dune/common/fvector.hh> 16 #include<dune/common/array.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> 49 template<
int k,
int d>
struct PkSize;
51 template<
int j,
int k,
int d>
59 template<
int k,
int d>
70 template<
int k,
int d>
103 template<
int k,
int d>
126 template<
int k,
int d>
177 if (d==1)
return k+1;
179 for (
int j=0; j<=k; j++)
201 alpha[
dim]=k; count++;
203 if (count==i)
return;
206 for (
int m=k-1; m>=0; m--)
211 if (count==i)
return;
222 for (
int j=0; j<d; j++) alpha[j] = 0;
224 for (
int k=1; k<25; k++)
231 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
238 for (
int i=0; i<d; ++i)
247 for (
int j = 0; j<d; ++j) {
248 alpha[j] = i % (k+1);
260 template <BasisType basisType>
266 template <
int k,
int d>
274 template <
int k,
int d>
297 template <
int k,
int d>
305 template <
int k,
int d>
340 for (
long i=k+1; i<=n; i++) nominator *= i;
342 for (
long i=2; i<=n-k; i++) denominator *= i;
343 return nominator/denominator;
348 for (
long i=n-k+1; i<=n; i++) nominator *= i;
350 for (
long i=2; i<=k; i++) denominator *= i;
351 return nominator/denominator;
361 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
368 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
374 template<
typename ComputationFieldType,
int d>
381 ComputationFieldType result(1.0);
382 for (
int i=0; i<d; i++)
384 ComputationFieldType exponent(a[i]+1);
385 result = result/exponent;
393 template<
typename ComputationFieldType>
400 ComputationFieldType one(1.0);
401 ComputationFieldType exponent0(a[0]+1);
402 return one/exponent0;
408 template<
typename ComputationFieldType>
415 ComputationFieldType sum(0.0);
416 for (
int k=0; k<=a[1]+1; k++)
420 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
421 ComputationFieldType denom(a[0]+k+1);
422 sum = sum + (nom/denom);
424 ComputationFieldType denom(a[1]+1);
431 template<
typename ComputationFieldType>
438 ComputationFieldType sumk(0.0);
439 for (
int k=0; k<=a[2]+1; k++)
443 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
444 ComputationFieldType denom(a[1]+k+1);
445 sumk = sumk + (nom/denom);
447 ComputationFieldType sumj(0.0);
448 for (
int j=0; j<=a[1]+a[2]+2; j++)
452 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
453 ComputationFieldType denom(a[0]+j+1);
454 sumj = sumj + (nom/denom);
456 ComputationFieldType denom(a[2]+1);
457 return sumk*sumj/denom;
466 template<
typename F,
int d>
469 template<
typename X,
typename A>
473 for (
int i=0; i<a[d]; i++)
482 template<
typename X,
typename A>
486 for (
int i=0; i<a[0]; i++)
521 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
532 : coeffs(new LowprecMat)
534 for (
int i=0; i<d; ++i)
535 gradcoeffs[i].reset(
new LowprecMat());
537 for (
int i=0; i<n; i++)
550 : coeffs(new LowprecMat)
552 for (
int i=0; i<d; ++i)
553 gradcoeffs[i].reset(
new LowprecMat());
555 for (
int i=0; i<n; i++)
568 return Traits::size(l,d);
572 template<
typename Po
int,
typename Result>
575 std::fill(r.begin(),r.end(),0.0);
576 for (
int j=0; j<n; ++j)
579 for (
int i=j; i<n; ++i)
580 r[i] += (*coeffs)[i][j] * monomial_value;
585 template<
typename Po
int,
typename Result>
588 std::fill(r.begin(),r.end(),0.0);
590 for (
int j=0; j<n; ++j)
593 for (
int i=j; i<n; ++i)
594 for (
int s=0;
s<d; ++
s)
595 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
600 template<
typename Po
int,
typename Result>
604 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
606 for (
int i=0; i<Traits::size(l,d); i++)
609 for (
int j=0; j<=i; j++)
616 template<
typename Po
int,
typename Result>
620 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
622 for (
int i=0; i<Traits::size(l,d); i++)
625 for (
int s=0;
s<d;
s++)
628 for (
int j=0; j<=i; j++)
631 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
637 std::array<std::shared_ptr<MultiIndex<d> >,n> alpha;
638 std::shared_ptr<LowprecMat> coeffs;
639 std::array<std::shared_ptr<LowprecMat>,d > gradcoeffs;
642 void orthonormalize()
657 for (
int s=0;
s<d;
s++)
658 for (
int i=0; i<n; i++)
659 for (
int j=0; j<=i; j++)
660 (*gradcoeffs[
s])[i][j] = 0;
661 for (
int i=0; i<n; i++)
662 for (
int j=0; j<=i; j++)
663 for (
int s=0; s<d; s++)
664 if ((*alpha[j])[s]>0)
667 FieldType factor = beta[
s];
669 int l = invert_index(beta);
670 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
689 for (
int i=0; i<n; i++)
692 for (
int j=0; j<d; j++)
693 if (a[j]!=(*alpha[i])[j]) found=
false;
696 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
702 HighprecMat *
p =
new HighprecMat();
706 for (
int i=0; i<n; i++)
707 for (
int j=0; j<n; j++)
709 c[i][j] = ComputationFieldType(1.0);
711 c[i][j] = ComputationFieldType(0.0);
715 for (
int i=0; i<n; i++)
718 ComputationFieldType bi[n];
721 for (
int j=0; j<i; j++)
724 bi[j] = ComputationFieldType(0.0);
725 for (
int l=0; l<=j; l++)
728 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
729 bi[j] = bi[j] + c[j][l]*integrator.
integrate(a);
731 for (
int l=0; l<=j; l++)
732 c[i][l] = c[i][l] - bi[j]*c[j][l];
736 ComputationFieldType s2(0.0);
738 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
740 for (
int j=0; j<i; j++)
741 s2 = s2 - bi[j]*bi[j];
742 ComputationFieldType
s(1.0);
745 for (
int l=0; l<=i; l++)
750 for (
int i=0; i<n; i++)
751 for (
int j=0; j<n; j++)
752 (*coeffs)[i][j] = c[i][j];
764 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
770 Dune::GeometryType gt;
773 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 0>
Traits;
781 unsigned int size ()
const {
return n; }
785 std::vector<typename Traits::RangeType>& out)
const {
793 std::vector<typename Traits::JacobianType>& out)
const {
803 Dune::GeometryType
type ()
const {
return gt; }
806 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
812 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
816 std::size_t
size ()
const {
return n; }
824 std::vector<Dune::LocalKey> li;
836 template<
typename F,
typename C>
840 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
842 typedef typename LB::Traits::RangeType RangeType;
843 const int d = LB::Traits::dimDomain;
844 const Dune::QuadratureRule<RealFieldType,d>&
845 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
849 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
852 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
853 it=rule.begin(); it!=rule.end(); ++it)
856 typename LB::Traits::DomainType x;
858 for (
int i=0; i<d; i++) x[i] = it->position()[i];
862 std::vector<RangeType> phi(LB::n);
863 lb.evaluateFunction(it->position(),phi);
866 for (
int i=0; i<LB::n; i++)
867 out[i] += y*phi[i]*it->weight();
872 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
875 Dune::GeometryType gt;
880 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
885 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
890 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
894 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
909 return interpolation;
912 Dune::GeometryType
type ()
const {
return gt; }
923 #endif // DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH Integrate monomials over the reference element.
Definition: l2orthonormal.hh:362
OPBLocalFiniteElement()
Definition: l2orthonormal.hh:884
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:531
MultiIndex()
Definition: l2orthonormal.hh:167
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:799
static const int dim
Definition: adaptivity.hh:83
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:882
Definition: l2orthonormal.hh:765
compute
Definition: l2orthonormal.hh:467
Definition: l2orthonormal.hh:49
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:811
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:185
const std::string s
Definition: function.hh:830
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:819
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:816
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:366
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:837
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:522
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:902
static int size(int k, int d)
Definition: l2orthonormal.hh:314
OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:893
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:220
Definition: l2orthonormal.hh:55
Definition: l2orthonormal.hh:807
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:245
unsigned int size() const
Definition: l2orthonormal.hh:781
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:586
int qk_size(int k, int d)
Definition: l2orthonormal.hh:235
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:784
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:527
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:331
Definition: l2orthonormal.hh:257
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:320
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:470
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 >, 0 > Traits
Definition: l2orthonormal.hh:773
int size(int l)
Definition: l2orthonormal.hh:566
static int size(int k, int d)
Definition: l2orthonormal.hh:282
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:436
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:833
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:897
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:779
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:413
Definition: l2orthonormal.hh:873
Definition: l2orthonormal.hh:104
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:907
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:889
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:483
double DefaultComputationalFieldType
Definition: l2orthonormal.hh:42
Definition: l2orthonormal.hh:127
Definition: l2orthonormal.hh:828
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:379
Dune::GeometryType type() const
Definition: l2orthonormal.hh:803
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:549
int pk_size(int k, int d)
Definition: l2orthonormal.hh:174
Dune::GeometryType type() const
Definition: l2orthonormal.hh:912
Definition: l2orthonormal.hh:261
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:398
Definition: l2orthonormal.hh:257
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:914
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:117
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:601
OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:776
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:195
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:288
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:573
Definition: l2orthonormal.hh:163
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:792
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:617
BasisType
Definition: l2orthonormal.hh:256
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:528
Definition: l2orthonormal.hh:52