dune-localfunctions  2.2.1
raviartthomas02dlocalinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_RT02DLOCALINTERPOLATION_HH
2 #define DUNE_RT02DLOCALINTERPOLATION_HH
3 
4 #include <cmath>
5 #include <vector>
6 #include <dune/common/exceptions.hh>
7 
8 namespace Dune
9 {
10  template<class LB>
12  {
13  public:
14 
17  {
18  sign0 = sign1 = sign2 = 1.0;
19  }
20 
22  RT02DLocalInterpolation (unsigned int s)
23  {
24  sign0 = sign1 = sign2 = 1.0;
25  if (s&1) sign0 *= -1.0;
26  if (s&2) sign1 *= -1.0;
27  if (s&4) sign2 *= -1.0;
28  m0[0] = 0.5; m0[1] = 0.0;
29  m1[0] = 0.0; m1[1] = 0.5;
30  m2[0] = 0.5; m2[1] = 0.5;
31  n0[0] = 0.0; n0[1] = -1.0;
32  n1[0] = -1.0; n1[1] = 0.0;
33  n2[0] = 1.0/sqrt(2.0); n2[1] = 1.0/sqrt(2.0);
34  c0 = ( 0.5*n0[0] - 1.0*n0[1]);
35  c1 = (-1.0*n1[0] + 0.5*n1[1]);
36  c2 = ( 0.5*n2[0] + 0.5*n2[1]);
37  }
38 
39  template<typename F, typename C>
40  void interpolate (const F& f, std::vector<C>& out) const
41  {
42  // f gives v*outer normal at a point on the edge!
43  typename F::Traits::RangeType y;
44 
45  out.resize(3);
46 
47  f.evaluate(m0,y); out[0] = (y[0]*n0[0]+y[1]*n0[1])*sign0/c0;
48  f.evaluate(m1,y); out[1] = (y[0]*n1[0]+y[1]*n1[1])*sign1/c1;
49  f.evaluate(m2,y); out[2] = (y[0]*n2[0]+y[1]*n2[1])*sign2/c2;
50  }
51 
52  private:
53  typename LB::Traits::RangeFieldType sign0,sign1,sign2;
54  typename LB::Traits::DomainType m0,m1,m2;
55  typename LB::Traits::DomainType n0,n1,n2;
56  typename LB::Traits::RangeFieldType c0,c1,c2;
57  };
58 }
59 
60 #endif