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_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_SIMPLEX2D_LOCALINTERPOLATION_HH
4 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_SIMPLEX2D_LOCALINTERPOLATION_HH
5 
6 #include <vector>
7 
8 #include <dune/geometry/quadraturerules.hh>
9 #include <dune/localfunctions/common/localinterpolation.hh>
10 
11 namespace Dune
12 {
13   /**
14    * \ingroup LocalInterpolationImplementation
15    * \brief First order Brezzi-Douglas-Marini shape functions on the reference triangle.
16    *
17    * \tparam LB corresponding LocalBasis giving traits
18    *
19    * \nosubgrouping
20    */
21   template<class LB>
22   class BDM1Simplex2DLocalInterpolation
23   {
24 
25   public:
26     //! \brief Standard constructor
BDM1Simplex2DLocalInterpolation()27     BDM1Simplex2DLocalInterpolation ()
28     {
29       sign0 = sign1 = sign2 = 1.0;
30     }
31 
32     /**
33      * \brief Make set number s, where 0 <= s < 8
34      *
35      * \param s Edge orientation indicator
36      */
BDM1Simplex2DLocalInterpolation(unsigned int s)37     BDM1Simplex2DLocalInterpolation (unsigned int s)
38     {
39       using std::sqrt;
40       sign0 = sign1 = sign2 = 1.0;
41       if (s & 1)
42       {
43         sign0 = -1.0;
44       }
45       if (s & 2)
46       {
47         sign1 = -1.0;
48       }
49       if (s & 4)
50       {
51         sign2 = -1.0;
52       }
53 
54       n0[0] =  0.0;
55       n0[1] = -1.0;
56       n1[0] = -1.0;
57       n1[1] =  0.0;
58       n2[0] =  1.0/sqrt(2.0);
59       n2[1] =  1.0/sqrt(2.0);
60       c0 =  0.5*n0[0] - 1.0*n0[1];
61       c1 = -1.0*n1[0] + 0.5*n1[1];
62       c2 =  0.5*n2[0] + 0.5*n2[1];
63     }
64 
65     /**
66      * \brief Interpolate a given function with shape functions
67      *
68      * \tparam F Function type for function which should be interpolated
69      * \tparam C Coefficient type
70      * \param ff function which should be interpolated
71      * \param out return value, vector of coefficients
72      */
73     template<typename F, typename C>
interpolate(const F & ff,std::vector<C> & out) const74     void interpolate (const F& ff, std::vector<C>& out) const
75     {
76       // f gives v*outer normal at a point on the edge!
77       typedef typename LB::Traits::RangeFieldType Scalar;
78 
79       auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
80 
81       out.resize(6);
82       fill(out.begin(), out.end(), 0.0);
83 
84       const int qOrder = 4;
85       const Dune::QuadratureRule<Scalar,1>& rule = Dune::QuadratureRules<Scalar,1>::rule(Dune::GeometryTypes::simplex(1), qOrder);
86 
87       for (typename Dune::QuadratureRule<Scalar,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
88       {
89         Scalar qPos = it->position();
90         typename LB::Traits::DomainType localPos;
91 
92         localPos[0] = qPos;
93         localPos[1] = 0.0;
94         auto y = f(localPos);
95         out[0] += (y[0]*n0[0] + y[1]*n0[1])*it->weight()*sign0/c0;
96         out[3] += (y[0]*n0[0] + y[1]*n0[1])*(2.0*qPos - 1.0)*it->weight()/c0;
97 
98         localPos[0] = 0.0;
99         localPos[1] = qPos;
100         y = f(localPos);
101         out[1] += (y[0]*n1[0] + y[1]*n1[1])*it->weight()*sign1/c1;
102         out[4] += (y[0]*n1[0] + y[1]*n1[1])*(1.0 - 2.0*qPos)*it->weight()/c1;
103 
104         localPos[0] = 1.0 - qPos;
105         localPos[1] = qPos;
106         y = f(localPos);
107         out[2] += (y[0]*n2[0] + y[1]*n2[1])*it->weight()*sign2/c2;
108         out[5] += (y[0]*n2[0] + y[1]*n2[1])*(2.0*qPos - 1.0)*it->weight()/c2;
109       }
110     }
111 
112   private:
113     typename LB::Traits::RangeFieldType sign0,sign1,sign2;
114     typename LB::Traits::DomainType n0,n1,n2;
115     typename LB::Traits::RangeFieldType c0,c1,c2;
116   };
117 }
118 
119 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI1_SIMPLEX2D_LOCALINTERPOLATION_HH
120