4 #ifndef DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
5 #define DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
13 #include<dune/common/fvector.hh>
14 #include<dune/common/fmatrix.hh>
15 #include<dune/common/gmpfield.hh>
16 #include<dune/common/exceptions.hh>
17 #include<dune/common/fvector.hh>
19 #include<dune/geometry/referenceelements.hh>
20 #include<dune/geometry/quadraturerules.hh>
21 #include<dune/geometry/type.hh>
23 #include<dune/localfunctions/common/localbasis.hh>
24 #include<dune/localfunctions/common/localkey.hh>
25 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
26 #include<dune/localfunctions/common/localinterpolation.hh>
68 for (
int j=0; j<=k; j++)
84 template<
int k,
int d>
100 alpha[
dim]=k; count++;
102 if (count==i)
return;
105 for (
int m=k-1; m>=0; m--)
110 if (count==i)
return;
121 for (
int j=0; j<d; j++) alpha[j] = 0;
123 for (
int k=1; k<25; k++)
130 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
137 for (
int i=0; i<d; ++i)
146 for (
int j = 0; j<d; ++j) {
147 alpha[j] = i % (k+1);
159 template <BasisType basisType>
165 template <
int k,
int d>
173 template <
int k,
int d>
196 template <
int k,
int d>
204 template <
int k,
int d>
239 for (
long i=k+1; i<=n; i++) nominator *= i;
241 for (
long i=2; i<=n-k; i++) denominator *= i;
242 return nominator/denominator;
247 for (
long i=n-k+1; i<=n; i++) nominator *= i;
249 for (
long i=2; i<=k; i++) denominator *= i;
250 return nominator/denominator;
260 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
267 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
273 template<
typename ComputationFieldType,
int d>
280 ComputationFieldType result(1.0);
281 for (
int i=0; i<d; i++)
283 ComputationFieldType exponent(a[i]+1);
284 result = result/exponent;
292 template<
typename ComputationFieldType>
299 ComputationFieldType one(1.0);
300 ComputationFieldType exponent0(a[0]+1);
301 return one/exponent0;
307 template<
typename ComputationFieldType>
314 ComputationFieldType sum(0.0);
315 for (
int k=0; k<=a[1]+1; k++)
319 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
320 ComputationFieldType denom(a[0]+k+1);
321 sum = sum + (nom/denom);
323 ComputationFieldType denom(a[1]+1);
330 template<
typename ComputationFieldType>
337 ComputationFieldType sumk(0.0);
338 for (
int k=0; k<=a[2]+1; k++)
342 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
343 ComputationFieldType denom(a[1]+k+1);
344 sumk = sumk + (nom/denom);
346 ComputationFieldType sumj(0.0);
347 for (
int j=0; j<=a[1]+a[2]+2; j++)
351 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
352 ComputationFieldType denom(a[0]+j+1);
353 sumj = sumj + (nom/denom);
355 ComputationFieldType denom(a[2]+1);
356 return sumk*sumj/denom;
365 template<
typename F,
int d>
368 template<
typename X,
typename A>
372 for (
int i=0; i<a[d]; i++)
381 template<
typename X,
typename A>
385 for (
int i=0; i<a[0]; i++)
420 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
433 for (
int i=0; i<d; ++i)
436 for (
int i=0; i<
n; i++)
451 for (
int i=0; i<d; ++i)
454 for (
int i=0; i<
n; i++)
467 return Traits::size(l,d);
471 template<
typename Po
int,
typename Result>
474 std::fill(r.begin(),r.end(),0.0);
475 for (
int j=0; j<
n; ++j)
478 for (
int i=j; i<
n; ++i)
479 r[i] += (*coeffs)[i][j] * monomial_value;
484 template<
typename Po
int,
typename Result>
487 std::fill(r.begin(),r.end(),0.0);
489 for (
int j=0; j<
n; ++j)
492 for (
int i=j; i<
n; ++i)
493 for (
int s=0;
s<d; ++
s)
494 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
499 template<
typename Po
int,
typename Result>
503 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
505 for (
int i=0; i<Traits::size(l,d); i++)
508 for (
int j=0; j<=i; j++)
515 template<
typename Po
int,
typename Result>
519 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
521 for (
int i=0; i<Traits::size(l,d); i++)
524 for (
int s=0;
s<d;
s++)
527 for (
int j=0; j<=i; j++)
530 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
536 std::array<std::shared_ptr<MultiIndex<d> >,
n> alpha;
537 std::shared_ptr<LowprecMat> coeffs;
538 std::array<std::shared_ptr<LowprecMat>,d > gradcoeffs;
541 void orthonormalize()
556 for (
int s=0;
s<d;
s++)
557 for (
int i=0; i<
n; i++)
558 for (
int j=0; j<=i; j++)
559 (*gradcoeffs[
s])[i][j] = 0;
560 for (
int i=0; i<
n; i++)
561 for (
int j=0; j<=i; j++)
562 for (
int s=0;
s<d;
s++)
563 if ((*alpha[j])[
s]>0)
566 FieldType factor = beta[
s];
568 int l = invert_index(beta);
569 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
586 int invert_index (MultiIndex<d>& a)
588 for (
int i=0; i<
n; i++)
591 for (
int j=0; j<d; j++)
592 if (a[j]!=(*alpha[i])[j]) found=
false;
595 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
605 for (
int i=0; i<
n; i++)
606 for (
int j=0; j<
n; j++)
608 c[i][j] = ComputationFieldType(1.0);
610 c[i][j] = ComputationFieldType(0.0);
613 MonomialIntegrator<ComputationFieldType,bt,d> integrator;
614 for (
int i=0; i<
n; i++)
617 ComputationFieldType bi[
n];
620 for (
int j=0; j<i; j++)
623 bi[j] = ComputationFieldType(0.0);
624 for (
int l=0; l<=j; l++)
627 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
628 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
630 for (
int l=0; l<=j; l++)
631 c[i][l] = c[i][l] - bi[j]*c[j][l];
635 ComputationFieldType s2(0.0);
637 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
638 s2 = s2 + integrator.integrate(a);
639 for (
int j=0; j<i; j++)
640 s2 = s2 - bi[j]*bi[j];
641 ComputationFieldType
s(1.0);
644 for (
int l=0; l<=i; l++)
649 for (
int i=0; i<
n; i++)
650 for (
int j=0; j<
n; j++)
651 (*coeffs)[i][j] = c[i][j];
663 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
669 Dune::GeometryType gt;
672 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
675 DUNE_NO_DEPRECATED_BEGIN
682 DUNE_NO_DEPRECATED_END
684 unsigned int size ()
const {
return n; }
688 std::vector<typename Traits::RangeType>& out)
const {
696 std::vector<typename Traits::JacobianType>& out)
const {
703 const typename Traits::DomainType& in,
704 std::vector<typename Traits::RangeType>& out)
const {
705 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
706 if (totalOrder == 0) {
709 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
718 Dune::GeometryType
type ()
const {
return gt; }
721 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
727 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
731 std::size_t
size ()
const {
return n; }
739 std::vector<Dune::LocalKey> li;
751 template<
typename F,
typename C>
755 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
757 typedef typename LB::Traits::RangeType RangeType;
758 const int d = LB::Traits::dimDomain;
759 const Dune::QuadratureRule<RealFieldType,d>&
760 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
762 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
766 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
769 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
770 it=rule.begin(); it!=rule.end(); ++it)
773 typename LB::Traits::DomainType x;
775 for (
int i=0; i<d; i++) x[i] = it->position()[i];
779 std::vector<RangeType> phi(LB::n);
780 lb.evaluateFunction(it->position(),phi);
783 for (
int i=0; i<LB::n; i++)
784 out[i] += y*phi[i]*it->weight();
789 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
792 Dune::GeometryType gt;
797 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
801 DUNE_NO_DEPRECATED_BEGIN
804 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
809 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
812 DUNE_NO_DEPRECATED_END
815 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
830 return interpolation;
840 Dune::GeometryType
type ()
const {
return gt; }
static const int dim
Definition: adaptivity.hh:84
const std::string s
Definition: function.hh:843
const P & p
Definition: constraints.hh:148
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
constexpr int qk_size(int k, int d)
Definition: l2orthonormal.hh:134
BasisType
Definition: l2orthonormal.hh:155
@ Qk
Definition: l2orthonormal.hh:156
@ Pk
Definition: l2orthonormal.hh:156
double DefaultComputationalFieldType
Definition: l2orthonormal.hh:44
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:230
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:144
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:74
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:119
constexpr int pk_size(int k, int d)
Definition: l2orthonormal.hh:63
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:94
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:121
Definition: l2orthonormal.hh:53
MultiIndex()
Definition: l2orthonormal.hh:56
Definition: l2orthonormal.hh:86
@ value
Definition: l2orthonormal.hh:88
Definition: l2orthonormal.hh:160
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:187
static int size(int k, int d)
Definition: l2orthonormal.hh:181
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:219
static int size(int k, int d)
Definition: l2orthonormal.hh:213
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:262
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:265
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:278
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:297
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:312
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:335
compute
Definition: l2orthonormal.hh:367
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:369
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:382
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:422
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:500
@ n
Definition: l2orthonormal.hh:425
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:430
int size(int l)
Definition: l2orthonormal.hh:465
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:516
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:472
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:448
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:485
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:427
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:426
Definition: l2orthonormal.hh:665
@ n
Definition: l2orthonormal.hh:673
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:714
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:695
DUNE_NO_DEPRECATED_END unsigned int size() const
Definition: l2orthonormal.hh:684
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:687
void partial(const std::array< unsigned int, Traits::dimDomain > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivative of all shape functions.
Definition: l2orthonormal.hh:702
Dune::GeometryType type() const
Definition: l2orthonormal.hh:718
DUNE_NO_DEPRECATED_BEGIN OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:677
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: l2orthonormal.hh:672
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:680
Definition: l2orthonormal.hh:723
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:731
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:726
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:734
Definition: l2orthonormal.hh:744
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:748
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:752
Definition: l2orthonormal.hh:791
std::size_t size() const
Definition: l2orthonormal.hh:835
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:808
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:818
Dune::GeometryType type() const
Definition: l2orthonormal.hh:840
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:799
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:828
DUNE_NO_DEPRECATED_BEGIN OPBLocalFiniteElement()
Definition: l2orthonormal.hh:803
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:823
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:842
DUNE_NO_DEPRECATED_END OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:814
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139