dune-pdelab  2.5-dev
pk1d.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 // Pk in one dimension with k as runtime variable
5 
6 #ifndef DUNE_PDELAB_FINITEELEMENT_PK1D_HH
7 #define DUNE_PDELAB_FINITEELEMENT_PK1D_HH
8 
9 #include <vector>
10 
11 #include <dune/common/fmatrix.hh>
12 #include <dune/geometry/type.hh>
13 
14 #include<dune/localfunctions/common/localbasis.hh>
15 #include<dune/localfunctions/common/localkey.hh>
16 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
17 
18 namespace Dune {
19 
28  template<class D, class R>
30  {
32  class Pk1dLocalBasis
33  {
34  Dune::GeometryType gt; // store geometry type for the basis
35  std::size_t k; // polynomial degree
36  std::size_t n; // the number of basis functions
37  std::vector<R> s; // Lagrange points on the reference interval
38  public:
39  typedef Dune::LocalBasisTraits<D,1,Dune::FieldVector<D,1>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,1>, 1> Traits;
40 
42  Pk1dLocalBasis (std::size_t k_) : gt(Dune::GeometryType::cube,1), k(k_), n(k_+1), s(n)
43  {
44  for (std::size_t i=0; i<=k; i++) s[i] = (1.0*i)/k;
45  }
46 
48  std::size_t size () const { return n; }
49 
51  inline void evaluateFunction (const typename Traits::DomainType& in,
52  std::vector<typename Traits::RangeType>& out) const {
53  out.resize(n);
54  for (std::size_t i=0; i<=k; i++)
55  {
56  out[i] = 1.0;
57  for (std::size_t j=0; j<=k; j++)
58  if (i!=j) out[i] *= (in[0]-s[j])/(s[i]-s[j]);
59  }
60  }
61 
63  inline void
64  evaluateJacobian (const typename Traits::DomainType& in,
65  std::vector<typename Traits::JacobianType>& out) const {
66  out.resize(n);
67  for (std::size_t i=0; i<=k; i++) // derivative of basis function i
68  {
69  out[i][0][0] = 0.0;
70  R factor = 1.0;
71  R denominator = 1.0;
72  for (std::size_t j=0; j<=k; j++)
73  {
74  if (j==i) continue; // treat factor (x-s_j)
75  denominator *= s[i]-s[j];
76  R a=1.0; // product of remaining factors (might be empty)
77  for (std::size_t l=j+1; l<=k; l++)
78  {
79  if (l==i) continue;
80  a *= in[0]-s[l];
81  }
82  out[i][0][0] += factor*a;
83  factor *= in[0]-s[j];
84  }
85  out[i][0][0] /= denominator;
86  }
87  }
88 
90  unsigned int order () const {
91  return k;
92  }
93 
95  Dune::GeometryType type () const { return gt; }
96  };
97 
99  class Pk1dLocalCoefficients
100  {
101  public:
102  Pk1dLocalCoefficients (std::size_t k_) : k(k_), n(k_+1), li(k_+1) {
103  li[0] = Dune::LocalKey(0,1,0);
104  for (int i=1; i<int(k); i++) li[i] = Dune::LocalKey(0,0,i-1);
105  li[k] = Dune::LocalKey(1,1,0);
106  }
107 
109  std::size_t size () const { return n; }
110 
112  const Dune::LocalKey& localKey (int i) const {
113  return li[i];
114  }
115 
116  private:
117  std::size_t k; // polynomial degree
118  std::size_t n; // the number of basis functions
119  std::vector<Dune::LocalKey> li; // assignment of basis function to subentities
120  };
121 
123  template<typename LB>
124  class Pk1dLocalInterpolation
125  {
126  public:
127  Pk1dLocalInterpolation (std::size_t k_) : k(k_), n(k_+1) {}
128 
130  template<typename F, typename C>
131  void interpolate (const F& f, std::vector<C>& out) const
132  {
133  out.resize(n);
134  typename LB::Traits::DomainType x;
135  typename LB::Traits::RangeType y;
136  for (int i=0; i<=int(k); i++)
137  {
138  x[0] = (1.0*i)/k; // the point to evaluate
139  f.evaluate(x,y);
140  out[i] = y[0];
141  }
142  }
143  private:
144  std::size_t k; // polynomial degree
145  std::size_t n; // the number of basis functions
146  };
147 
148  Dune::GeometryType gt;
149  Pk1dLocalBasis basis;
150  Pk1dLocalCoefficients coefficients;
151  Pk1dLocalInterpolation<Pk1dLocalBasis> interpolation;
152 
153  public:
154  typedef Dune::LocalFiniteElementTraits<Pk1dLocalBasis,
155  Pk1dLocalCoefficients,
156  Pk1dLocalInterpolation<Pk1dLocalBasis> > Traits;
157 
158  Pk1dLocalFiniteElement (std::size_t k)
159  : gt(Dune::GeometryType::cube,1), basis(k), coefficients(k), interpolation(k)
160  {}
161 
162  const typename Traits::LocalBasisType& localBasis () const
163  {
164  return basis;
165  }
166 
167  const typename Traits::LocalCoefficientsType& localCoefficients () const
168  {
169  return coefficients;
170  }
171 
172  const typename Traits::LocalInterpolationType& localInterpolation () const
173  {
174  return interpolation;
175  }
176 
177  Dune::GeometryType type () const { return gt; }
178 
180  return new Pk1dLocalFiniteElement(*this);
181  }
182  };
183 }
184 
185 #endif // DUNE_PDELAB_FINITEELEMENT_PK1D_HH
Define the Pk Lagrange basis functions in 1d on the reference interval.
Definition: pk1d.hh:29
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: pk1d.hh:167
Dune::LocalFiniteElementTraits< Pk1dLocalBasis, Pk1dLocalCoefficients, Pk1dLocalInterpolation< Pk1dLocalBasis > > Traits
Definition: pk1d.hh:156
Dune::GeometryType type() const
Definition: pk1d.hh:177
const Traits::LocalInterpolationType & localInterpolation() const
Definition: pk1d.hh:172
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
const Traits::LocalBasisType & localBasis() const
Definition: pk1d.hh:162
Pk1dLocalFiniteElement * clone() const
Definition: pk1d.hh:179
Pk1dLocalFiniteElement(std::size_t k)
Definition: pk1d.hh:158
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:204
const std::string s
Definition: function.hh:830