4#ifndef DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
5#define DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
19#include<dune/geometry/referenceelements.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>
42 typedef GMPField<512> DefaultComputationalFieldType;
44 typedef double DefaultComputationalFieldType;
52 class MultiIndex :
public std::array<int,d>
56 MultiIndex () :
std::array<int,d>()
63 constexpr int pk_size (
int k,
int d)
68 for (
int j=0; j<=k; j++)
69 n += pk_size(k-j,d-1);
74 inline int pk_size_exact (
int k,
int d)
79 return pk_size(k,d)-pk_size(k-1,d);
84 template<
int k,
int d>
94 void pk_enumerate_multiindex (MultiIndex<d>& alpha,
int& count,
int norm,
int dim,
int k,
int i)
100 alpha[dim]=k; count++;
102 if (count==i)
return;
105 for (
int m=k-1; m>=0; m--)
109 pk_enumerate_multiindex(alpha,count,norm+m,dim+1,k-m,i);
110 if (count==i)
return;
119 void pk_multiindex (
int i, MultiIndex<d>& alpha)
121 for (
int j=0; j<d; j++) alpha[j] = 0;
123 for (
int k=1; k<25; k++)
124 if (i>=pk_size(k-1,d) && i<pk_size(k,d))
127 pk_enumerate_multiindex(alpha,count,0,0,k,i-pk_size(k-1,d));
134 constexpr int qk_size (
int k,
int d)
137 for (
int i=0; i<d; ++i)
144 void qk_multiindex (
int i,
int k, MultiIndex<d>& alpha)
146 for (
int j = 0; j<d; ++j) {
147 alpha[j] = i % (k+1);
159 template <BasisType basisType>
163 struct BasisTraits<BasisType::Pk>
165 template <
int k,
int d>
173 template <
int k,
int d>
181 static int size(
int k,
int d)
187 static void multiindex(
int i,
int k, MultiIndex<d>& alpha)
189 pk_multiindex(i,alpha);
194 struct BasisTraits<BasisType::Qk>
196 template <
int k,
int d>
204 template <
int k,
int d>
213 static int size(
int k,
int d)
219 static void multiindex(
int i,
int k, MultiIndex<d>& alpha)
221 return qk_multiindex(i,k,alpha);
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>
265 ComputationFieldType
integrate (
const MultiIndex<d>& a)
const
273 template<
typename ComputationFieldType,
int d>
278 ComputationFieldType
integrate (
const MultiIndex<d>& a)
const
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>
297 ComputationFieldType
integrate (
const MultiIndex<1>& a)
const
299 ComputationFieldType one(1.0);
300 ComputationFieldType exponent0(a[0]+1);
301 return one/exponent0;
307 template<
typename ComputationFieldType>
312 ComputationFieldType
integrate (
const MultiIndex<2>& a)
const
314 ComputationFieldType sum(0.0);
315 for (
int k=0; k<=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>
335 ComputationFieldType
integrate (
const MultiIndex<3>& a)
const
337 ComputationFieldType sumk(0.0);
338 for (
int k=0; k<=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++)
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>
369 static F compute (
const X& x,
const A& a)
372 for (
int i=0; i<a[d]; i++)
381 template<
typename X,
typename A>
382 static F compute (
const X& x,
const A& 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>
423 typedef BasisTraits<basisType> Traits;
425 enum{ n = Traits::template Size<k,d>::value };
433 for (
int i=0; i<d; ++i)
436 for (
int i=0; i<n; i++)
438 alpha[i].reset(
new MultiIndex<d>());
439 Traits::multiindex(i,k,*alpha[i]);
451 for (
int i=0; i<d; ++i)
454 for (
int i=0; i<n; i++)
456 alpha[i].reset(
new MultiIndex<d>());
457 Traits::multiindex(i,k,*alpha[i]);
467 return Traits::size(l,d);
471 template<
typename Po
int,
typename Result>
472 void evaluateFunction (
const Point& x, Result& r)
const
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>
485 void evaluateJacobian (
const Point& x, Result& r)
const
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>
500 void evaluateFunction (
int l,
const Point& x, Result& r)
const
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>
516 void evaluateJacobian (
int l,
const Point& x, Result& r)
const
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)
565 MultiIndex<d> beta = *alpha[j];
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;
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);
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];
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>
666 typedef PB::BasisTraits<basisType> BasisTraits;
672 typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,
Dune::FieldVector<R,1>,
Dune::FieldMatrix<R,1,d> > Traits;
673 enum{ n = BasisTraits::template Size<k,d>::value };
677 OPBLocalBasis (
int order_) : opb(),
gt(bt,d) {}
680 OPBLocalBasis (
int order_,
const LFE & lfe) : opb(lfe),
gt(bt,d) {}
684 unsigned int size ()
const {
return n; }
688 std::vector<typename Traits::RangeType>& out)
const {
696 std::vector<typename Traits::JacobianType>& out)
const {
702 void partial(
const std::array<unsigned int, Traits::dimDomain>& order,
704 std::vector<typename Traits::RangeType>& out)
const {
706 if (totalOrder == 0) {
707 evaluateFunction(in, out);
709 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
714 unsigned int order ()
const {
715 return BasisTraits::template Order<k,d>::value;
721 template<
int k,
int d, PB::BasisType basisType = PB::BasisType::Pk>
722 class OPBLocalCoefficients
724 enum{ n = PB::BasisTraits<basisType>::template Size<k,d>::value };
726 OPBLocalCoefficients (
int order_) : li(n) {
731 std::size_t
size ()
const {
return n; }
739 std::vector<Dune::LocalKey> li;
743 class OPBLocalInterpolation
748 OPBLocalInterpolation (
const LB& lb_,
int order_) : lb(lb_) {}
751 template<
typename F,
typename C>
752 void interpolate (
const F& ff, std::vector<C>& out)
const
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;
762 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
766 for (
int i=0; i<LB::n; i++) out[i] = 0.0;
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>
790 class OPBLocalFiniteElement
793 OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> basis;
794 OPBLocalCoefficients<k,d,basisType> coefficients;
795 OPBLocalInterpolation<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> > interpolation;
798 OPBLocalCoefficients<k,d,basisType>,
799 OPBLocalInterpolation<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> > > Traits;
803 OPBLocalFiniteElement ()
804 :
gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
808 explicit OPBLocalFiniteElement (
const LFE & lfe)
809 :
gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
814 OPBLocalFiniteElement (
const OPBLocalFiniteElement & other)
815 :
gt(other.
gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
830 return interpolation;
835 std::size_t
size()
const
842 OPBLocalFiniteElement* clone ()
const {
843 return new OPBLocalFiniteElement(*
this);
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:122
@ simplex
Simplicial element in any nonnegative dimension.
Definition: type.hh:121
Describe position of one degree of freedom.
Definition: localkey.hh:24
Default exception for dummy implementations.
Definition: exceptions.hh:263
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
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:262
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:265
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:422
Definition: polynomialbasis.hh:65
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: polynomialbasis.hh:119
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:127
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
Default exception class for range errors.
Definition: exceptions.hh:254
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Wrapper for the GNU multiprecision (GMP) library.
#define DUNE_NO_DEPRECATED_END
Ignore deprecation warnings (end)
Definition: deprecated.hh:38
#define DUNE_NO_DEPRECATED_BEGIN
Ignore deprecation warnings (start)
Definition: deprecated.hh:32
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:131
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
compute
Definition: l2orthonormal.hh:367
A unique label for each type of element that can occur in a grid.