dune-pdelab  2.4.1
pk1dbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_NPDE_PK1D_HH
2 #define DUNE_NPDE_PK1D_HH
3 
4 #include<vector>
5 #include<iostream>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/fmatrix.hh>
8 #include<dune/common/exceptions.hh>
9 
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/referenceelements.hh>
12 #include<dune/geometry/quadraturerules.hh>
13 
14 #include<dune/localfunctions/common/localbasis.hh>
15 #include<dune/localfunctions/common/localkey.hh>
16 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
17 
19 
20 namespace Dune {
21 
27  template<class D, class R>
29  {
31  class Pk1dLocalBasis
32  {
33  Dune::GeometryType gt; // store geometry type for the basis
34  std::size_t k; // polynomial degree
35  std::size_t n; // the number of basis functions
36  std::vector<R> s; // Lagrange points on the reference interval
37  public:
38  typedef Dune::LocalBasisTraits<D,1,Dune::FieldVector<D,1>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,1>, 1> Traits;
39 
41  Pk1dLocalBasis (std::size_t k_) : gt(Dune::GeometryType::cube,1), k(k_), n(k_+1), s(n)
42  {
43  for (std::size_t i=0; i<=k; i++) s[i] = (1.0*i)/k;
44  }
45 
47  std::size_t size () const { return n; }
48 
50  inline void evaluateFunction (const typename Traits::DomainType& in,
51  std::vector<typename Traits::RangeType>& out) const {
52  out.resize(n);
53  for (std::size_t i=0; i<=k; i++)
54  {
55  out[i] = 1.0;
56  for (std::size_t j=0; j<=k; j++)
57  if (i!=j) out[i] *= (in[0]-s[j])/(s[i]-s[j]);
58  }
59  }
60 
62  inline void
63  evaluateJacobian (const typename Traits::DomainType& in,
64  std::vector<typename Traits::JacobianType>& out) const {
65  out.resize(n);
66  for (std::size_t i=0; i<=k; i++) // derivative of basis function i
67  {
68  out[i][0][0] = 0.0;
69  R factor = 1.0;
70  R denominator = 1.0;
71  for (std::size_t j=0; j<=k; j++)
72  {
73  if (j==i) continue; // treat factor (x-s_j)
74  denominator *= s[i]-s[j];
75  R a=1.0; // product of remaining factors (might be empty)
76  for (std::size_t l=j+1; l<=k; l++)
77  {
78  if (l==i) continue;
79  a *= in[0]-s[l];
80  }
81  out[i][0][0] += factor*a;
82  factor *= in[0]-s[j];
83  }
84  out[i][0][0] /= denominator;
85  }
86  }
87 
89  unsigned int order () const {
90  return k;
91  }
92 
94  Dune::GeometryType type () const { return gt; }
95  };
96 
98  class Pk1dLocalCoefficients
99  {
100  public:
101  Pk1dLocalCoefficients (std::size_t k_) : k(k_), n(k_+1), li(k_+1) {
102  li[0] = Dune::LocalKey(0,1,0);
103  for (int i=1; i<int(k); i++) li[i] = Dune::LocalKey(0,0,i-1);
104  li[k] = Dune::LocalKey(1,1,0);
105  }
106 
108  std::size_t size () const { return n; }
109 
111  const Dune::LocalKey& localKey (int i) const {
112  return li[i];
113  }
114 
115  private:
116  std::size_t k; // polynomial degree
117  std::size_t n; // the number of basis functions
118  std::vector<Dune::LocalKey> li; // assignment of basis function to subentities
119  };
120 
122  template<typename LB>
123  class Pk1dLocalInterpolation
124  {
125  public:
126  Pk1dLocalInterpolation (std::size_t k_) : k(k_), n(k_+1) {}
127 
129  template<typename F, typename C>
130  void interpolate (const F& f, std::vector<C>& out) const
131  {
132  out.resize(n);
133  typename LB::Traits::DomainType x;
134  typename LB::Traits::RangeType y;
135  for (int i=0; i<=int(k); i++)
136  {
137  x[0] = (1.0*i)/k; // the point to evaluate
138  f.evaluate(x,y);
139  out[i] = y[0];
140  }
141  }
142  private:
143  std::size_t k; // polynomial degree
144  std::size_t n; // the number of basis functions
145  };
146 
147  Dune::GeometryType gt;
148  Pk1dLocalBasis basis;
149  Pk1dLocalCoefficients coefficients;
150  Pk1dLocalInterpolation<Pk1dLocalBasis> interpolation;
151 
152  public:
153  typedef Dune::LocalFiniteElementTraits<Pk1dLocalBasis,
154  Pk1dLocalCoefficients,
155  Pk1dLocalInterpolation<Pk1dLocalBasis> > Traits;
156 
157  Pk1dLocalFiniteElement (std::size_t k)
158  : gt(Dune::GeometryType::cube,1), basis(k), coefficients(k), interpolation(k)
159  {}
160 
161  const typename Traits::LocalBasisType& localBasis () const
162  {
163  return basis;
164  }
165 
166  const typename Traits::LocalCoefficientsType& localCoefficients () const
167  {
168  return coefficients;
169  }
170 
171  const typename Traits::LocalInterpolationType& localInterpolation () const
172  {
173  return interpolation;
174  }
175 
176  Dune::GeometryType type () const { return gt; }
177 
179  return new Pk1dLocalFiniteElement(*this);
180  }
181  };
182 
183  namespace PDELab {
184 
190  template<class D, class R>
192  : public Dune::PDELab::SimpleLocalFiniteElementMap< Pk1dLocalFiniteElement<D,R> >
193  {
194  public:
197  , _k(k)
198  {}
199 
200  bool fixedSize() const
201  {
202  return true;
203  }
204 
205  bool hasDOFs(int codim) const
206  {
207  switch (codim)
208  {
209  case 0:
210  return _k != 1;
211  case 1:
212  return _k > 0;
213  }
214  return false;
215  }
216 
217  std::size_t size(GeometryType gt) const
218  {
219  if (gt.isVertex())
220  return _k > 0 ? 1 : 0;
221  if (gt.isLine())
222  return _k > 0 ? _k - 1 : 1;
223  return 0;
224  }
225 
226  std::size_t maxLocalSize() const
227  {
228  return _k + 1;
229  }
230 
231  private:
232  const std::size_t _k;
233  };
234  }
235 }
236 #endif
Define the Pk Lagrange basis functions in 1d on the reference interval.
Definition: pk1dbasis.hh:28
Dune::LocalFiniteElementTraits< Pk1dLocalBasis, Pk1dLocalCoefficients, Pk1dLocalInterpolation< Pk1dLocalBasis > > Traits
Definition: pk1dbasis.hh:155
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:191
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: pk1dbasis.hh:166
Definition: adaptivity.hh:27
Pk1dLocalFiniteElementMap(std::size_t k)
Definition: pk1dbasis.hh:195
const Traits::LocalBasisType & localBasis() const
Definition: pk1dbasis.hh:161
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:107
Dune::GeometryType type() const
Definition: pk1dbasis.hh:176
std::size_t size(GeometryType gt) const
Definition: pk1dbasis.hh:217
bool fixedSize() const
Definition: pk1dbasis.hh:200
FiniteElementMap for the Pk basis in 1d.
Definition: pk1dbasis.hh:191
const std::string s
Definition: function.hh:1102
Pk1dLocalFiniteElement * clone() const
Definition: pk1dbasis.hh:178
bool hasDOFs(int codim) const
Definition: pk1dbasis.hh:205
Pk1dLocalFiniteElement(std::size_t k)
Definition: pk1dbasis.hh:157
std::size_t maxLocalSize() const
Definition: pk1dbasis.hh:226
const Traits::LocalInterpolationType & localInterpolation() const
Definition: pk1dbasis.hh:171