dune-localfunctions  2.3.0
dualq1localbasis.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_Q1_LOCALBASIS_HH
4 #define DUNE_DUAL_Q1_LOCALBASIS_HH
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 
10 
11 namespace Dune
12 {
22  template<class D, class R, int dim>
24  {
25  public:
26  typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>,
27  Dune::FieldMatrix<R,1,dim> > Traits;
28 
29  void setCoefficients(const Dune::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)>& coefficients)
30  {
31  coefficients_ = coefficients;
32  }
33 
35  unsigned int size () const
36  {
37  return 1<<dim;
38  }
39 
41  inline void evaluateFunction (const typename Traits::DomainType& in,
42  std::vector<typename Traits::RangeType>& out) const
43  {
44  // compute q1 values
45  std::vector<typename Traits::RangeType> q1Values(size());
46 
47  for (size_t i=0; i<size(); i++) {
48 
49  q1Values[i] = 1;
50 
51  for (int j=0; j<dim; j++)
52  // if j-th bit of i is set multiply with in[j], else with 1-in[j]
53  q1Values[i] *= (i & (1<<j)) ? in[j] : 1-in[j];
54 
55  }
56 
57  // compute the dual values by using that they are linear combinations of q1 functions
58  out.resize(size());
59  for (size_t i=0; i<size(); i++)
60  out[i] = 0;
61 
62  for (size_t i=0; i<size(); i++)
63  for (size_t j=0; j<size(); j++)
64  out[i] += coefficients_[i][j]*q1Values[j];
65 
66 
67  }
68 
70  inline void
71  evaluateJacobian (const typename Traits::DomainType& in, // position
72  std::vector<typename Traits::JacobianType>& out) const // return value
73  {
74  // compute q1 jacobians
75  std::vector<typename Traits::JacobianType> q1Jacs(size());
76 
77  // Loop over all shape functions
78  for (size_t i=0; i<size(); i++) {
79 
80  // Loop over all coordinate directions
81  for (int j=0; j<dim; j++) {
82 
83  // Initialize: the overall expression is a product
84  // if j-th bit of i is set to -1, else 1
85  q1Jacs[i][0][j] = (i & (1<<j)) ? 1 : -1;
86 
87  for (int k=0; k<dim; k++) {
88 
89  if (j!=k)
90  // if k-th bit of i is set multiply with in[j], else with 1-in[j]
91  q1Jacs[i][0][j] *= (i & (1<<k)) ? in[k] : 1-in[k];
92 
93  }
94 
95  }
96 
97  }
98 
99  // compute the dual jacobians by using that they are linear combinations of q1 functions
100  out.resize(size());
101  for (size_t i=0; i<size(); i++)
102  out[i] = 0;
103 
104  for (size_t i=0; i<size(); i++)
105  for (size_t j=0; j<size(); j++)
106  out[i].axpy(coefficients_[i][j],q1Jacs[j]);
107 
108  }
109 
111  unsigned int order () const
112  {
113  return 1;
114  }
115 
116  private:
117  Dune::array<Dune::FieldVector<R, (1<<dim)> ,(1<<dim)> coefficients_;
118  };
119 }
120 #endif
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:39
LocalBasisTraits< D, dim, Dune::FieldVector< D, dim >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, dim > > Traits
Definition: dualq1localbasis.hh:27
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: dualq1localbasis.hh:71
unsigned int order() const
Polynomial order of the shape functions.
Definition: dualq1localbasis.hh:111
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: dualq1localbasis.hh:41
void setCoefficients(const Dune::array< Dune::FieldVector< R,(1<< dim)>,(1<< dim)> &coefficients)
Definition: dualq1localbasis.hh:29
Dual Lagrange shape functions of order 1 on the reference cube.
Definition: dualq1localbasis.hh:23
unsigned int size() const
number of shape functions
Definition: dualq1localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:51