dune-localfunctions  2.3.1
dualp1localbasis.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 #ifndef DUNE_DUAL_P1_LOCALBASIS_HH
4 #define DUNE_DUAL_P1_LOCALBASIS_HH
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
9 
10 namespace Dune
11 {
23  template<class D, class R, int dim>
25  {
26  public:
28  typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
29  Dune::FieldMatrix<R,1,dim> > Traits;
30 
32  unsigned int size () const
33  {
34  return dim+1;
35  }
36 
38  inline void evaluateFunction (const typename Traits::DomainType& in,
39  std::vector<typename Traits::RangeType>& out) const
40  {
41  // evaluate P1 basis functions
42  std::vector<typename Traits::RangeType> p1Values(size());
43 
44  p1Values[0] = 1.0;
45 
46  for (int i=0; i<dim; i++) {
47  p1Values[0] -= in[i];
48  p1Values[i+1] = in[i];
49  }
50 
51  // compute dual basis function values as a linear combination of the Lagrange values
52  out.resize(size());
53 
54  for (int i=0; i<=dim; i++) {
55  out[i] = (dim+1)*p1Values[i];
56  for (int j=0; j<i; j++)
57  out[i] -= p1Values[j];
58 
59  for (int j=i+1; j<=dim; j++)
60  out[i] -= p1Values[j];
61  }
62  }
63 
65  inline void
66  evaluateJacobian (const typename Traits::DomainType& in,
67  std::vector<typename Traits::JacobianType>& out) const
68  {
69  // evaluate P1 jacobians
70  std::vector<typename Traits::JacobianType> p1Jacs(size());
71 
72  for (int i=0; i<dim; i++)
73  p1Jacs[0][0][i] = -1;
74 
75  for (int i=0; i<dim; i++)
76  for (int j=0; j<dim; j++)
77  p1Jacs[i+1][0][j] = (i==j);
78 
79  // compute dual basis jacobians as linear combination of the Lagrange jacobians
80  out.resize(size());
81 
82  for (size_t i=0; i<=dim; i++) {
83  out[i][0] = 0;
84  out[i][0].axpy((dim+1),p1Jacs[i][0]);
85 
86  for (size_t j=0; j<i; j++)
87  out[i][0] -= p1Jacs[j][0];
88 
89  for (int j=i+1; j<=dim; j++)
90  out[i][0] -= p1Jacs[j][0];
91  }
92  }
93 
95  unsigned int order () const
96  {
97  return 1;
98  }
99  };
100 }
101 #endif