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>
67 for (
int j=0; j<=k; j++)
83 template<
int k,
int d>
99 alpha[
dim]=k; count++;
101 if (count==i)
return;
104 for (
int m=k-1; m>=0; m--)
109 if (count==i)
return;
120 for (
int j=0; j<d; j++) alpha[j] = 0;
122 for (
int k=1; k<25; k++)
129 DUNE_THROW(Dune::NotImplemented,
"Polynomial degree greater than 24 in pk_multiindex");
136 for (
int i=0; i<d; ++i)
145 for (
int j = 0; j<d; ++j) {
146 alpha[j] = i % (k+1);
158 template <BasisType basisType>
164 template <
int k,
int d>
172 template <
int k,
int d>
195 template <
int k,
int d>
203 template <
int k,
int d>
238 for (
long i=k+1; i<=n; i++) nominator *= i;
240 for (
long i=2; i<=n-k; i++) denominator *= i;
241 return nominator/denominator;
246 for (
long i=n-k+1; i<=n; i++) nominator *= i;
248 for (
long i=2; i<=k; i++) denominator *= i;
249 return nominator/denominator;
259 template<
typename ComputationFieldType, Dune::GeometryType::BasicType bt,
int d>
266 DUNE_THROW(Dune::NotImplemented,
"non-specialized version of MonomalIntegrator called. Please implement.");
272 template<
typename ComputationFieldType,
int d>
279 ComputationFieldType result(1.0);
280 for (
int i=0; i<d; i++)
282 ComputationFieldType exponent(a[i]+1);
283 result = result/exponent;
291 template<
typename ComputationFieldType>
298 ComputationFieldType one(1.0);
299 ComputationFieldType exponent0(a[0]+1);
300 return one/exponent0;
306 template<
typename ComputationFieldType>
313 ComputationFieldType sum(0.0);
314 for (
int k=0; k<=a[1]+1; k++)
318 ComputationFieldType nom(sign*
binomial(a[1]+1,k));
319 ComputationFieldType denom(a[0]+k+1);
320 sum = sum + (nom/denom);
322 ComputationFieldType denom(a[1]+1);
329 template<
typename ComputationFieldType>
336 ComputationFieldType sumk(0.0);
337 for (
int k=0; k<=a[2]+1; k++)
341 ComputationFieldType nom(sign*
binomial(a[2]+1,k));
342 ComputationFieldType denom(a[1]+k+1);
343 sumk = sumk + (nom/denom);
345 ComputationFieldType sumj(0.0);
346 for (
int j=0; j<=a[1]+a[2]+2; j++)
350 ComputationFieldType nom(sign*
binomial(a[1]+a[2]+2,j));
351 ComputationFieldType denom(a[0]+j+1);
352 sumj = sumj + (nom/denom);
354 ComputationFieldType denom(a[2]+1);
355 return sumk*sumj/denom;
364 template<
typename F,
int d>
367 template<
typename X,
typename A>
371 for (
int i=0; i<a[d]; i++)
380 template<
typename X,
typename A>
384 for (
int i=0; i<a[0]; i++)
419 template<
typename FieldType,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
432 for (
int i=0; i<d; ++i)
435 for (
int i=0; i<
n; i++)
450 for (
int i=0; i<d; ++i)
453 for (
int i=0; i<
n; i++)
466 return Traits::size(l,d);
470 template<
typename Po
int,
typename Result>
473 std::fill(r.begin(),r.end(),0.0);
474 for (
int j=0; j<
n; ++j)
477 for (
int i=j; i<
n; ++i)
478 r[i] += (*coeffs)[i][j] * monomial_value;
483 template<
typename Po
int,
typename Result>
486 std::fill(r.begin(),r.end(),0.0);
488 for (
int j=0; j<
n; ++j)
491 for (
int i=j; i<
n; ++i)
492 for (
int s=0;
s<d; ++
s)
493 r[i][0][
s] += (*gradcoeffs[
s])[i][j]*monomial_value;
498 template<
typename Po
int,
typename Result>
502 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
504 for (
int i=0; i<Traits::size(l,d); i++)
507 for (
int j=0; j<=i; j++)
514 template<
typename Po
int,
typename Result>
518 DUNE_THROW(Dune::RangeError,
"l>k in OrthonormalPolynomialBasis::evaluateFunction");
520 for (
int i=0; i<Traits::size(l,d); i++)
523 for (
int s=0;
s<d;
s++)
526 for (
int j=0; j<=i; j++)
529 for (
int s=0;
s<d;
s++) r[i][0][
s] = sum[
s];
535 std::array<std::shared_ptr<MultiIndex<d> >,
n> alpha;
536 std::shared_ptr<LowprecMat> coeffs;
537 std::array<std::shared_ptr<LowprecMat>,d > gradcoeffs;
540 void orthonormalize()
555 for (
int s=0;
s<d;
s++)
556 for (
int i=0; i<
n; i++)
557 for (
int j=0; j<=i; j++)
558 (*gradcoeffs[
s])[i][j] = 0;
559 for (
int i=0; i<
n; i++)
560 for (
int j=0; j<=i; j++)
561 for (
int s=0;
s<d;
s++)
562 if ((*alpha[j])[
s]>0)
565 FieldType factor = beta[
s];
567 int l = invert_index(beta);
568 (*gradcoeffs[
s])[i][l] += (*coeffs)[i][j]*factor;
585 int invert_index (MultiIndex<d>& a)
587 for (
int i=0; i<
n; i++)
590 for (
int j=0; j<d; j++)
591 if (a[j]!=(*alpha[i])[j]) found=
false;
594 DUNE_THROW(Dune::RangeError,
"index not found in invertindex");
604 for (
int i=0; i<
n; i++)
605 for (
int j=0; j<
n; j++)
607 c[i][j] = ComputationFieldType(1.0);
609 c[i][j] = ComputationFieldType(0.0);
612 MonomialIntegrator<ComputationFieldType,bt,d> integrator;
613 for (
int i=0; i<
n; i++)
616 ComputationFieldType bi[
n];
619 for (
int j=0; j<i; j++)
622 bi[j] = ComputationFieldType(0.0);
623 for (
int l=0; l<=j; l++)
626 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
627 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
629 for (
int l=0; l<=j; l++)
630 c[i][l] = c[i][l] - bi[j]*c[j][l];
634 ComputationFieldType s2(0.0);
636 for (
int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
637 s2 = s2 + integrator.integrate(a);
638 for (
int j=0; j<i; j++)
639 s2 = s2 - bi[j]*bi[j];
640 ComputationFieldType
s(1.0);
643 for (
int l=0; l<=i; l++)
648 for (
int i=0; i<
n; i++)
649 for (
int j=0; j<
n; j++)
650 (*coeffs)[i][j] = c[i][j];
662 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
668 Dune::GeometryType gt;
671 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> >
Traits;
674 DUNE_NO_DEPRECATED_BEGIN
681 DUNE_NO_DEPRECATED_END
683 unsigned int size ()
const {
return n; }
687 std::vector<typename Traits::RangeType>& out)
const {
695 std::vector<typename Traits::JacobianType>& out)
const {
701 void partial(
const std::array<unsigned int, Traits::dimDomain>& DUNE_UNUSED(
order),
702 const typename Traits::DomainType& DUNE_UNUSED(in),
703 std::vector<typename Traits::RangeType>& DUNE_UNUSED(out))
const {
704 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
705 if (totalOrder == 0) {
708 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
717 Dune::GeometryType
type ()
const {
return gt; }
720 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
726 for (
int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
730 std::size_t
size ()
const {
return n; }
738 std::vector<Dune::LocalKey> li;
750 template<
typename F,
typename C>
754 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
756 typedef typename LB::Traits::RangeType RangeType;
757 const int d = LB::Traits::dimDomain;
758 const Dune::QuadratureRule<RealFieldType,d>&
759 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
763 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
766 for (
typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
767 it=rule.begin(); it!=rule.end(); ++it)
770 typename LB::Traits::DomainType x;
772 for (
int i=0; i<d; i++) x[i] = it->position()[i];
776 std::vector<RangeType> phi(LB::n);
777 lb.evaluateFunction(it->position(),phi);
780 for (
int i=0; i<LB::n; i++)
781 out[i] += y*phi[i]*it->weight();
786 template<
class D,
class R,
int k,
int d, Dune::GeometryType::BasicType bt,
typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
789 Dune::GeometryType gt;
794 typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
798 DUNE_NO_DEPRECATED_BEGIN
801 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
806 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
809 DUNE_NO_DEPRECATED_END
812 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
827 return interpolation;
837 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:133
BasisType
Definition: l2orthonormal.hh:154
@ Qk
Definition: l2orthonormal.hh:155
@ Pk
Definition: l2orthonormal.hh:155
double DefaultComputationalFieldType
Definition: l2orthonormal.hh:43
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:229
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:143
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:73
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:118
constexpr int pk_size(int k, int d)
Definition: l2orthonormal.hh:62
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:93
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:119
Definition: l2orthonormal.hh:52
MultiIndex()
Definition: l2orthonormal.hh:55
Definition: l2orthonormal.hh:85
@ value
Definition: l2orthonormal.hh:87
Definition: l2orthonormal.hh:159
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:186
static int size(int k, int d)
Definition: l2orthonormal.hh:180
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:218
static int size(int k, int d)
Definition: l2orthonormal.hh:212
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:261
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:264
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:277
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:296
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:311
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:334
compute
Definition: l2orthonormal.hh:366
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:368
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:381
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:421
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:499
@ n
Definition: l2orthonormal.hh:424
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:429
int size(int l)
Definition: l2orthonormal.hh:464
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:515
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:471
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:447
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:484
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:426
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:425
Definition: l2orthonormal.hh:664
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:713
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:694
DUNE_NO_DEPRECATED_END unsigned int size() const
Definition: l2orthonormal.hh:683
@ n
Definition: l2orthonormal.hh:672
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:686
void partial(const std::array< unsigned int, Traits::dimDomain > &DUNE_UNUSED(order), const typename Traits::DomainType &DUNE_UNUSED(in), std::vector< typename Traits::RangeType > &DUNE_UNUSED(out)) const
Evaluate partial derivative of all shape functions.
Definition: l2orthonormal.hh:701
Dune::GeometryType type() const
Definition: l2orthonormal.hh:717
DUNE_NO_DEPRECATED_BEGIN OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:676
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: l2orthonormal.hh:671
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:679
Definition: l2orthonormal.hh:722
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:730
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:725
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:733
Definition: l2orthonormal.hh:743
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:751
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:747
Definition: l2orthonormal.hh:788
std::size_t size() const
Definition: l2orthonormal.hh:832
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:805
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:815
Dune::GeometryType type() const
Definition: l2orthonormal.hh:837
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:796
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:825
DUNE_NO_DEPRECATED_BEGIN OPBLocalFiniteElement()
Definition: l2orthonormal.hh:800
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:820
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:839
DUNE_NO_DEPRECATED_END OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:811
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139