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_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH
4 #define DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH
5 
6 #include <array>
7 #include <bitset>
8 #include <numeric>
9 #include <vector>
10 
11 #include <dune/common/fmatrix.hh>
12 
13 #include "../../common/localbasis.hh"
14 
15 namespace Dune
16 {
17   /**
18    * \ingroup LocalBasisImplementation
19    * \brief First order Brezzi-Douglas-Marini shape functions on quadrilaterals.
20    *
21    * \tparam D Type to represent the field in the domain.
22    * \tparam R Type to represent the field in the range.
23    *
24    * \nosubgrouping
25    */
26   template<class D, class R>
27   class BDM2Simplex2DLocalBasis
28   {
29 
30   public:
31     typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,
32         R,2,Dune::FieldVector<R,2>,
33         Dune::FieldMatrix<R,2,2> > Traits;
34 
35     //! \brief Standard constructor
BDM2Simplex2DLocalBasis()36     BDM2Simplex2DLocalBasis()
37     {
38       for (size_t i=0; i<3; i++)
39         sign_[i] = 1.0;
40     }
41 
42     /**
43      * \brief Make set number s, where 0 <= s < 8
44      *
45      * \param s Edge orientation indicator
46      */
BDM2Simplex2DLocalBasis(std::bitset<3> s)47     BDM2Simplex2DLocalBasis(std::bitset<3> s)
48     {
49       for (size_t i=0; i<3; i++)
50         sign_[i] = s[i] ? -1.0 : 1.0;
51     }
52 
53     //! \brief number of shape functions
size() const54     unsigned int size() const
55     {
56       return 12;
57     }
58 
59     /**
60      * \brief Evaluate all shape functions
61      *
62      * \param in Position
63      * \param out return value
64      */
evaluateFunction(const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const65     inline void evaluateFunction(const typename Traits::DomainType& in,
66                                  std::vector<typename Traits::RangeType>& out) const
67     {
68       out.resize(size());
69 
70       out[0][0] = sign_[0]*(-2*in[0]*in[1] + in[0]*in[0]);
71       out[0][1] = sign_[0]*(-1 + 6*in[1] -2*in[0]*in[1] - 5*in[1]*in[1]);
72 
73       out[1][0] = 1.5*in[0] + 3*in[0]*in[1] - 4.5*in[0]*in[0];
74       out[1][1] = -3 + 6*in[0] + 10.5*in[1] - 15*in[0]*in[1] - 7.5*in[1]*in[1];
75 
76       out[2][0] = sign_[0]*(-7.5*in[0] + 5*in[0]*in[1] + 12.5*in[0]*in[0]);
77       out[2][1] = sign_[0]*(-5 + 30*in[0] + 7.5*in[1] - 25*in[0]*in[1] - 30*in[0]*in[0] - 2.5*in[1]*in[1]);
78 
79 
80 
81       out[3][0] = sign_[1]*(-1 + 6*in[0] - 2*in[0]*in[1] - 5*in[0]*in[0]);
82       out[3][1] = sign_[1]*(-2*in[0]*in[1] + in[1]*in[1]);
83 
84       out[4][0] = 3 - 10.5*in[0] - 6*in[1] + 15*in[0]*in[1] + 7.5*in[0]*in[0];
85       out[4][1] = -1.5*in[1] - 3*in[0]*in[1] + 4.5*in[1]*in[1];
86 
87       out[5][0] = sign_[1]*(-5 + 7.5*in[0] + 30*in[1] - 25*in[0]*in[1] - 2.5*in[0]*in[0] - 30*in[1]*in[1]);
88       out[5][1] = sign_[1]*(-7.5*in[1] + 5*in[0]*in[1] + 12.5*in[1]*in[1]);
89 
90 
91 
92       out[6][0] = sign_[2]*(-3*in[0] + 4*in[0]*in[1] + 4*in[0]*in[0]);
93       out[6][1] = sign_[2]*(-3*in[1] + 4*in[0]*in[1] + 4*in[1]*in[1]);
94 
95       out[7][0] = -3*in[0] + 6*in[0]*in[0];
96       out[7][1] = 3*in[1] - 6*in[1]*in[1];
97 
98       out[8][0] = sign_[2]*(-10*in[0]*in[1] + 5*in[0]*in[0]);
99       out[8][1] = sign_[2]*(-10*in[0]*in[1] + 5*in[1]*in[1]);
100 
101 
102 
103       out[9][0] = 18*in[0] - 12*in[0]*in[1] - 18*in[0]*in[0];
104       out[9][1] = 6*in[1] - 12*in[0]*in[1] - 6*in[1]*in[1];
105 
106       out[10][0] = 6*in[0] - 12*in[0]*in[1] - 6*in[0]*in[0];
107       out[10][1] = 18*in[1] - 12*in[0]*in[1] - 18*in[1]*in[1];
108 
109       out[11][0] = 90*in[0] - 180*in[0]*in[1] - 90*in[0]*in[0];
110       out[11][1] = -90*in[1] + 180*in[0]*in[1] + 90*in[1]*in[1];
111     }
112 
113     /**
114      * \brief Evaluate Jacobian of all shape functions
115      *
116      * \param in Position
117      * \param out return value
118      */
evaluateJacobian(const typename Traits::DomainType & in,std::vector<typename Traits::JacobianType> & out) const119     inline void evaluateJacobian(const typename Traits::DomainType& in,
120                                  std::vector<typename Traits::JacobianType>& out) const
121     {
122       out.resize(size());
123 
124       out[0][0][0] = sign_[0]*(-2*in[1] + 2*in[0]);
125       out[0][0][1] = sign_[0]*(-2*in[0]);
126 
127       out[0][1][0] = sign_[0]*(-2*in[1]);
128       out[0][1][1] = sign_[0]*(6 -2*in[0] - 10*in[1]);
129 
130 
131       out[1][0][0] = 1.5 + 3*in[1] - 9*in[0];
132       out[1][0][1] = 3*in[0];
133 
134       out[1][1][0] = 6 - 15*in[1];
135       out[1][1][1] = 10.5 - 15*in[0] - 15*in[1];
136 
137 
138       out[2][0][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]);
139       out[2][0][1] = sign_[0]*(5*in[0]);
140 
141       out[2][1][0] = sign_[0]*(30 - 25*in[1] - 60*in[0]);
142       out[2][1][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]);
143 
144 
145 
146       out[3][0][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]);
147       out[3][0][1] = sign_[1]*(-2*in[0]);
148 
149       out[3][1][0] = sign_[1]*(-2*in[1]);
150       out[3][1][1] = sign_[1]*(-2*in[0] + 2*in[1]);
151 
152 
153       out[4][0][0] = -10.5 + 15*in[1] + 15*in[0];
154       out[4][0][1] = -6 + 15*in[0];
155 
156       out[4][1][0] = -3*in[1];
157       out[4][1][1] = -1.5 - 3*in[0] + 9*in[1];
158 
159 
160       out[5][0][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]);
161       out[5][0][1] = sign_[1]*(30 - 25*in[0] - 60*in[1]);
162 
163       out[5][1][0] = sign_[1]*(5*in[1]);
164       out[5][1][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]);
165 
166 
167 
168       out[6][0][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]);
169       out[6][0][1] = sign_[2]*(4*in[0]);
170 
171       out[6][1][0] = sign_[2]*(4*in[1]);
172       out[6][1][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]);
173 
174 
175       out[7][0][0] = -3 + 12*in[0];
176       out[7][0][1] = 0;
177 
178       out[7][1][0] = 0;
179       out[7][1][1] = 3 - 12*in[1];
180 
181 
182       out[8][0][0] = sign_[2]*(-10*in[1] + 10*in[0]);
183       out[8][0][1] = sign_[2]*(-10*in[0]);
184 
185       out[8][1][0] = sign_[2]*(-10*in[1]);
186       out[8][1][1] = sign_[2]*(-10*in[0] + 10*in[1]);
187 
188 
189       out[9][0][0] = 18 - 12*in[1] - 36*in[0];
190       out[9][0][1] = -12*in[0];
191 
192       out[9][1][0] = -12*in[1];
193       out[9][1][1] = 6 - 12*in[0] - 12*in[1];
194 
195       out[10][0][0] = 6 - 12*in[1] - 12*in[0];
196       out[10][0][1] = -12*in[0];
197 
198       out[10][1][0] = -12*in[1];
199       out[10][1][1] = 18 - 12*in[0] - 36*in[1];
200 
201       out[11][0][0] = 90 - 180*in[1] - 180*in[0];
202       out[11][0][1] = -180*in[0];
203 
204       out[11][1][0] = 180*in[1];
205       out[11][1][1] = -90 + 180*in[0] + 180*in[1];
206     }
207 
208     //! \brief Evaluate partial derivatives of all shape functions
partial(const std::array<unsigned int,2> & order,const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const209     void partial (const std::array<unsigned int, 2>& order,
210                   const typename Traits::DomainType& in,         // position
211                   std::vector<typename Traits::RangeType>& out) const      // return value
212     {
213       auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
214       if (totalOrder == 0) {
215         evaluateFunction(in, out);
216       } else if (totalOrder == 1) {
217         out.resize(size());
218         auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
219 
220         switch (direction) {
221         case 0:
222           out[0][0] = sign_[0]*(-2*in[1] + 2*in[0]);
223           out[0][1] = sign_[0]*(-2*in[1]);
224 
225           out[1][0] = 1.5 + 3*in[1] - 9*in[0];
226           out[1][1] = 6 - 15*in[1];
227 
228           out[2][0] = sign_[0]*(-7.5 + 5*in[1] + 25*in[0]);
229           out[2][1] = sign_[0]*(30 - 25*in[1] - 60*in[0]);
230 
231           out[3][0] = sign_[1]*(6 - 2*in[1] - 10*in[0]);
232           out[3][1] = sign_[1]*(-2*in[1]);
233 
234           out[4][0] = -10.5 + 15*in[1] + 15*in[0];
235           out[4][1] = -3*in[1];
236 
237           out[5][0] = sign_[1]*(7.5 - 25*in[1] - 5*in[0]);
238           out[5][1] = sign_[1]*(5*in[1]);
239 
240           out[6][0] = sign_[2]*(-3 + 4*in[1] + 8*in[0]);
241           out[6][1] = sign_[2]*(4*in[1]);
242 
243           out[7][0] = -3 + 12*in[0];
244           out[7][1] = 0;
245 
246           out[8][0] = sign_[2]*(-10*in[1] + 10*in[0]);
247           out[8][1] = sign_[2]*(-10*in[1]);
248 
249           out[9][0] = 18 - 12*in[1] - 36*in[0];
250           out[9][1] = -12*in[1];
251 
252           out[10][0] = 6 - 12*in[1] - 12*in[0];
253           out[10][1] = -12*in[1];
254 
255           out[11][0] = 90 - 180*in[1] - 180*in[0];
256           out[11][1] = 180*in[1];
257           break;
258         case 1:
259           out[0][0] = sign_[0]*(-2*in[0]);
260           out[0][1] = sign_[0]*(6 -2*in[0] - 10*in[1]);
261 
262           out[1][0] = 3*in[0];
263           out[1][1] = 10.5 - 15*in[0] - 15*in[1];
264 
265           out[2][0] = sign_[0]*(5*in[0]);
266           out[2][1] = sign_[0]*(7.5 - 25*in[0] - 5*in[1]);
267 
268           out[3][0] = sign_[1]*(-2*in[0]);
269           out[3][1] = sign_[1]*(-2*in[0] + 2*in[1]);
270 
271           out[4][0] = -6 + 15*in[0];
272           out[4][1] = -1.5 - 3*in[0] + 9*in[1];
273 
274           out[5][0] = sign_[1]*(30 - 25*in[0] - 60*in[1]);
275           out[5][1] = sign_[1]*(-7.5 + 5*in[0] + 25*in[1]);
276 
277           out[6][0] = sign_[2]*(4*in[0]);
278           out[6][1] = sign_[2]*(-3 + 4*in[0] + 8*in[1]);
279 
280           out[7][0] = 0;
281           out[7][1] = 3 - 12*in[1];
282 
283           out[8][0] = sign_[2]*(-10*in[0]);
284           out[8][1] = sign_[2]*(-10*in[0] + 10*in[1]);
285 
286           out[9][0] = -12*in[0];
287           out[9][1] = 6 - 12*in[0] - 12*in[1];
288 
289           out[10][0] = -12*in[0];
290           out[10][1] = 18 - 12*in[0] - 36*in[1];
291 
292           out[11][0] = -180*in[0];
293           out[11][1] = -90 + 180*in[0] + 180*in[1];
294           break;
295         default:
296           DUNE_THROW(RangeError, "Component out of range.");
297         }
298       } else {
299         DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
300       }
301     }
302 
303     //! \brief Polynomial order of the shape functions
order() const304     unsigned int order() const
305     {
306       return 2; // TODO: check whether this is not order 3
307     }
308 
309   private:
310     std::array<R,3> sign_;
311   };
312 } // end namespace Dune
313 #endif // DUNE_LOCALFUNCTIONS_BREZZIDOUGLASMARINI2_SIMPLEX2D_LOCALBASIS_HH
314