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 
25  template<class D, class R>
27  {
29  class Pk1dLocalBasis
30  {
31  Dune::GeometryType gt; // store geometry type for the basis
32  std::size_t k; // polynomial degree
33  std::size_t n; // the number of basis functions
34  std::vector<R> s; // Lagrange points on the reference interval
35  public:
36  typedef Dune::LocalBasisTraits<D,1,Dune::FieldVector<D,1>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,1>, 1> Traits;
37 
39  Pk1dLocalBasis (std::size_t k_) : gt(Dune::GeometryType::cube,1), k(k_), n(k_+1), s(n)
40  {
41  for (std::size_t i=0; i<=k; i++) s[i] = (1.0*i)/k;
42  }
43 
45  std::size_t size () const { return n; }
46 
48  inline void evaluateFunction (const typename Traits::DomainType& in,
49  std::vector<typename Traits::RangeType>& out) const {
50  out.resize(n);
51  for (std::size_t i=0; i<=k; i++)
52  {
53  out[i] = 1.0;
54  for (std::size_t j=0; j<=k; j++)
55  if (i!=j) out[i] *= (in[0]-s[j])/(s[i]-s[j]);
56  }
57  }
58 
60  inline void
61  evaluateJacobian (const typename Traits::DomainType& in,
62  std::vector<typename Traits::JacobianType>& out) const {
63  out.resize(n);
64  for (std::size_t i=0; i<=k; i++) // derivative of basis function i
65  {
66  out[i][0][0] = 0.0;
67  R factor = 1.0;
68  R denominator = 1.0;
69  for (std::size_t j=0; j<=k; j++)
70  {
71  if (j==i) continue; // treat factor (x-s_j)
72  denominator *= s[i]-s[j];
73  R a=1.0; // product of remaining factors (might be empty)
74  for (std::size_t l=j+1; l<=k; l++)
75  {
76  if (l==i) continue;
77  a *= in[0]-s[l];
78  }
79  out[i][0][0] += factor*a;
80  factor *= in[0]-s[j];
81  }
82  out[i][0][0] /= denominator;
83  }
84  }
85 
87  unsigned int order () const {
88  return k;
89  }
90 
92  Dune::GeometryType type () const { return gt; }
93  };
94 
96  class Pk1dLocalCoefficients
97  {
98  public:
99  Pk1dLocalCoefficients (std::size_t k_) : k(k_), n(k_+1), li(k_+1) {
100  li[0] = Dune::LocalKey(0,1,0);
101  for (int i=1; i<int(k); i++) li[i] = Dune::LocalKey(0,0,i-1);
102  li[k] = Dune::LocalKey(1,1,0);
103  }
104 
106  std::size_t size () const { return n; }
107 
109  const Dune::LocalKey& localKey (int i) const {
110  return li[i];
111  }
112 
113  private:
114  std::size_t k; // polynomial degree
115  std::size_t n; // the number of basis functions
116  std::vector<Dune::LocalKey> li; // assignment of basis function to subentities
117  };
118 
120  template<typename LB>
121  class Pk1dLocalInterpolation
122  {
123  public:
124  Pk1dLocalInterpolation (std::size_t k_) : k(k_), n(k_+1) {}
125 
127  template<typename F, typename C>
128  void interpolate (const F& f, std::vector<C>& out) const
129  {
130  out.resize(n);
131  typename LB::Traits::DomainType x;
132  typename LB::Traits::RangeType y;
133  for (int i=0; i<=int(k); i++)
134  {
135  x[0] = (1.0*i)/k; // the point to evaluate
136  f.evaluate(x,y);
137  out[i] = y[0];
138  }
139  }
140  private:
141  std::size_t k; // polynomial degree
142  std::size_t n; // the number of basis functions
143  };
144 
145  Dune::GeometryType gt;
146  Pk1dLocalBasis basis;
147  Pk1dLocalCoefficients coefficients;
148  Pk1dLocalInterpolation<Pk1dLocalBasis> interpolation;
149 
150  public:
151  typedef Dune::LocalFiniteElementTraits<Pk1dLocalBasis,
152  Pk1dLocalCoefficients,
153  Pk1dLocalInterpolation<Pk1dLocalBasis> > Traits;
154 
155  Pk1dLocalFiniteElement (std::size_t k)
156  : gt(Dune::GeometryType::cube,1), basis(k), coefficients(k), interpolation(k)
157  {}
158 
159  const typename Traits::LocalBasisType& localBasis () const
160  {
161  return basis;
162  }
163 
164  const typename Traits::LocalCoefficientsType& localCoefficients () const
165  {
166  return coefficients;
167  }
168 
169  const typename Traits::LocalInterpolationType& localInterpolation () const
170  {
171  return interpolation;
172  }
173 
174  Dune::GeometryType type () const { return gt; }
175 
177  return new Pk1dLocalFiniteElement(*this);
178  }
179  };
180 }
181 
182 #endif // DUNE_PDELAB_FINITEELEMENT_PK1D_HH
const std::string s
Definition: function.hh:830
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: pk1d.hh:164
Pk1dLocalFiniteElement * clone() const
Definition: pk1d.hh:176
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:204
Dune::LocalFiniteElementTraits< Pk1dLocalBasis, Pk1dLocalCoefficients, Pk1dLocalInterpolation< Pk1dLocalBasis > > Traits
Definition: pk1d.hh:153
Define the Pk Lagrange basis functions in 1d on the reference interval.
Definition: pk1d.hh:26
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Pk1dLocalFiniteElement(std::size_t k)
Definition: pk1d.hh:155
const Traits::LocalBasisType & localBasis() const
Definition: pk1d.hh:159
const Traits::LocalInterpolationType & localInterpolation() const
Definition: pk1d.hh:169
Dune::GeometryType type() const
Definition: pk1d.hh:174