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_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH
4 #define DUNE_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH
5 
6 /** \file
7     \brief Hierarchical prism p2 shape functions for the simplex
8  */
9 
10 #include <numeric>
11 
12 #include <dune/common/fvector.hh>
13 #include <dune/common/fmatrix.hh>
14 
15 #include <dune/localfunctions/common/localbasis.hh>
16 
17 namespace Dune
18 {
19   template<class D, class R>
20   class HierarchicalPrismP2LocalBasis
21   {
22   public:
23     //! \brief export type traits for function signature
24     typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>, Dune::FieldMatrix<R,1,3> > Traits;
25 
26     //! \brief number of shape functions
size() const27     unsigned int size () const
28     {
29       return 18;
30     }
31 
32     //! \brief Evaluate all shape functions
evaluateFunction(const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const33     void evaluateFunction (const typename Traits::DomainType& in,
34                            std::vector<typename Traits::RangeType> & out) const
35     {
36       out.resize(18);
37 
38       out[0]=(1.0-in[0]-in[1])*(1.0-in[2]);
39       out[1]= in[0]*(1-in[2]);
40       out[2]=in[1]*(1-in[2]);
41       out[3]=in[2]*(1.0-in[0]-in[1]);
42       out[4]=in[0]*in[2];
43       out[5]=in[1]*in[2];
44 
45       //edges
46       out[6]=2*(1.0-in[0]-in[1])*(0.5-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]);
47       out[7]=2*in[0]*(-0.5+in[0])*(4*in[2]-4*in[2]*in[2]);
48       out[8]=2*in[1]*(-0.5+in[1])*(4*in[2]-4*in[2]*in[2]);
49       out[9]=4*in[0]*(1-in[0]-in[1])*(1-3*in[2]+2*in[2]*in[2]);
50       out[10]=4*in[1]*(1-in[0]-in[1])*(1-3*in[2]+2*in[2]*in[2]);
51       out[11]=4*in[0]*in[1]*(1-3*in[2]+2*in[2]*in[2]);
52       out[12]=4*in[0]*(1-in[0]-in[1])*(-in[2]+2*in[2]*in[2]);
53       out[13]=4*in[1]*(1-in[0]-in[1])*(-in[2]+2*in[2]*in[2]);
54       out[14]=4*in[0]*in[1]*(-in[2]+2*in[2]*in[2]);
55 
56       //faces
57       out[15]=4*in[0]*(1-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]);
58       out[16]=4*in[1]*(1-in[0]-in[1])*(4*in[2]-4*in[2]*in[2]);
59       out[17]=4*in[0]*in[1]*(4*in[2]-4*in[2]*in[2]);
60     }
61 
62 
63 
64     //! \brief Evaluate Jacobian of all shape functions
evaluateJacobian(const typename Traits::DomainType & in,std::vector<typename Traits::JacobianType> & out) const65     void evaluateJacobian (const typename Traits::DomainType& in,     //position
66                            std::vector<typename Traits::JacobianType>& out) const  //return  value
67     {
68       out.resize(18);
69 
70       //vertices
71       out[0][0][0] = in[2]-1;
72       out[0][0][1] = in[2]-1;
73       out[0][0][2] = in[0]+in[1]-1;
74 
75       out[1][0][0] = 1-in[2];
76       out[1][0][1] = 0;
77       out[1][0][2] =-in[0];
78 
79       out[2][0][0] = 0;
80       out[2][0][1] = 1-in[2];
81       out[2][0][2] = -in[1];
82 
83       out[3][0][0] = -in[2];
84       out[3][0][1] = -in[2];
85       out[3][0][2] = 1-in[0]-in[1];
86 
87       out[4][0][0] = in[2];
88       out[4][0][1] = 0;
89       out[4][0][2] = in[0];
90 
91       out[5][0][0] = 0;
92       out[5][0][1] = in[2];
93       out[5][0][2] = in[1];
94 
95       //edges
96       out[6][0][0] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]);
97       out[6][0][1] = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]);
98       out[6][0][2] = 2*(1-in[0]-in[1])*(0.5-in[0]-in[1])*(4-8*in[2]);
99 
100       out[7][0][0] = (-1+4*in[0])*(4*in[2]-4*in[2]*in[2]);
101       out[7][0][1] = 0;
102       out[7][0][2] = 2*in[0]*(-0.5+in[0])*(4-8*in[2]);
103 
104       out[8][0][0] = 0;
105       out[8][0][1] = (-1+4*in[1])*(4*in[2]-4*in[2]*in[2]);
106       out[8][0][2] = 2*in[1]*(-0.5+in[1])*(4-8*in[2]);
107 
108       out[9][0][0] = (4-8*in[0]-4*in[1])*(1-3*in[2]+2*in[2]*in[2]);
109       out[9][0][1] = -4*in[0]*(1-3*in[2]+2*in[2]*in[2]);
110       out[9][0][2] = 4*in[0]*(1-in[0]-in[1])*(-3+4*in[2]);
111 
112       out[10][0][0] = (-4*in[1])*(1-3*in[2]+2*in[2]*in[2]);
113       out[10][0][1] = (4-4*in[0]-8*in[1])*(1-3*in[2]+2*in[2]*in[2]);
114       out[10][0][2] = 4*in[1]*(1-in[0]-in[1])*(-3+4*in[2]);
115 
116       out[11][0][0] = 4*in[1]*(1-3*in[2]+2*in[2]*in[2]);
117       out[11][0][1] = 4*in[0]*(1-3*in[2]+2*in[2]*in[2]);
118       out[11][0][2] = 4*in[0]*in[1]*(-3+4*in[2]);
119 
120       out[12][0][0] = (4-8*in[0]-4*in[1])*(-in[2]+2*in[2]*in[2]);
121       out[12][0][1] = (-4*in[0])*(-in[2]+2*in[2]*in[2]);
122       out[12][0][2] = 4*in[0]*(1-in[0]-in[1])*(-1+4*in[2]);
123 
124       out[13][0][0] = -4*in[1]*(-in[2]+2*in[2]*in[2]);
125       out[13][0][1] = (4-4*in[0]-8*in[1])*(-in[2]+2*in[2]*in[2]);
126       out[13][0][2] = 4*in[1]*(1-in[0]-in[1])*(-1+4*in[2]);
127 
128       out[14][0][0] = 4*in[1]*(-in[2]+2*in[2]*in[2]);
129       out[14][0][1] = 4*in[0]*(-in[2]+2*in[2]*in[2]);
130       out[14][0][2] = 4*in[0]*in[1]*(-1+4*in[2]);
131 
132       //faces
133       out[15][0][0] = (4-8*in[0]-4*in[1])*(4*in[2]-4*in[2]*in[2]);
134       out[15][0][1] = -4*in[0]*(4*in[2]-4*in[2]*in[2]);
135       out[15][0][2] = 4*in[0]*(1-in[0]-in[1])*(4-8*in[2]);
136 
137       out[16][0][0] = -4*in[1]*(4*in[2]-4*in[2]*in[2]);
138       out[16][0][1] = (4-4*in[0]-8*in[1])*(4*in[2]-4*in[2]*in[2]);
139       out[16][0][2] = 4*in[1]*(1-in[0]-in[1])*(4-8*in[2]);
140 
141       out[17][0][0] = 4*in[1]*(4*in[2]-4*in[2]*in[2]);
142       out[17][0][1] = 4*in[0]*(4*in[2]-4*in[2]*in[2]);
143       out[17][0][2] = 4*in[0]*in[1]*(4-8*in[2]);
144     }
145 
146     //! \brief Evaluate partial derivatives of all shape functions
partial(const std::array<unsigned int,3> & order,const typename Traits::DomainType & in,std::vector<typename Traits::RangeType> & out) const147     void partial (const std::array<unsigned int, 3>& order,
148                   const typename Traits::DomainType& in,         // position
149                   std::vector<typename Traits::RangeType>& out) const      // return value
150     {
151       auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
152       if (totalOrder == 0) {
153         evaluateFunction(in, out);
154       } else if (totalOrder == 1) {
155         out.resize(size());
156         auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
157 
158         switch (direction) {
159         case 0:
160           out[0]  = in[2]-1;
161           out[1]  = 1-in[2];
162           out[2]  = 0;
163           out[3]  = -in[2];
164           out[4]  = in[2];
165           out[5]  = 0;
166           out[6]  = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]);
167           out[7]  = (-1+4*in[0])*(4*in[2]-4*in[2]*in[2]);
168           out[8]  = 0;
169           out[9]  = (4-8*in[0]-4*in[1])*(1-3*in[2]+2*in[2]*in[2]);
170           out[10] = (-4*in[1])*(1-3*in[2]+2*in[2]*in[2]);
171           out[11] = 4*in[1]*(1-3*in[2]+2*in[2]*in[2]);
172           out[12] = (4-8*in[0]-4*in[1])*(-in[2]+2*in[2]*in[2]);
173           out[13] = -4*in[1]*(-in[2]+2*in[2]*in[2]);
174           out[14] = 4*in[1]*(-in[2]+2*in[2]*in[2]);
175           out[15] = (4-8*in[0]-4*in[1])*(4*in[2]-4*in[2]*in[2]);
176           out[16] = -4*in[1]*(4*in[2]-4*in[2]*in[2]);
177           out[17] = 4*in[1]*(4*in[2]-4*in[2]*in[2]);
178           break;
179         case 1:
180           out[0]  = in[2]-1;
181           out[1]  = 0;
182           out[2]  = 1-in[2];
183           out[3]  = -in[2];
184           out[4]  = 0;
185           out[5]  = in[2];
186           out[6]  = (-3+4*in[0]+4*in[1])*(4*in[2]-4*in[2]*in[2]);
187           out[7]  = 0;
188           out[8]  = (-1+4*in[1])*(4*in[2]-4*in[2]*in[2]);
189           out[9]  = -4*in[0]*(1-3*in[2]+2*in[2]*in[2]);
190           out[10] = (4-4*in[0]-8*in[1])*(1-3*in[2]+2*in[2]*in[2]);
191           out[11] = 4*in[0]*(1-3*in[2]+2*in[2]*in[2]);
192           out[12] = (-4*in[0])*(-in[2]+2*in[2]*in[2]);
193           out[13] = (4-4*in[0]-8*in[1])*(-in[2]+2*in[2]*in[2]);
194           out[14] = 4*in[0]*(-in[2]+2*in[2]*in[2]);
195           out[15] = -4*in[0]*(4*in[2]-4*in[2]*in[2]);
196           out[16] = (4-4*in[0]-8*in[1])*(4*in[2]-4*in[2]*in[2]);
197           out[17] = 4*in[0]*(4*in[2]-4*in[2]*in[2]);
198           break;
199         case 2:
200           out[0]  = in[0]+in[1]-1;
201           out[1]  =-in[0];
202           out[2]  = -in[1];
203           out[3]  = 1-in[0]-in[1];
204           out[4]  = in[0];
205           out[5]  = in[1];
206           out[6]  = 2*(1-in[0]-in[1])*(0.5-in[0]-in[1])*(4-8*in[2]);
207           out[7]  = 2*in[0]*(-0.5+in[0])*(4-8*in[2]);
208           out[8]  = 2*in[1]*(-0.5+in[1])*(4-8*in[2]);
209           out[9]  = 4*in[0]*(1-in[0]-in[1])*(-3+4*in[2]);
210           out[10] = 4*in[1]*(1-in[0]-in[1])*(-3+4*in[2]);
211           out[11] = 4*in[0]*in[1]*(-3+4*in[2]);
212           out[12] = 4*in[0]*(1-in[0]-in[1])*(-1+4*in[2]);
213           out[13] = 4*in[1]*(1-in[0]-in[1])*(-1+4*in[2]);
214           out[14] = 4*in[0]*in[1]*(-1+4*in[2]);
215           out[15] = 4*in[0]*(1-in[0]-in[1])*(4-8*in[2]);
216           out[16] = 4*in[1]*(1-in[0]-in[1])*(4-8*in[2]);
217           out[17] = 4*in[0]*in[1]*(4-8*in[2]);
218           break;
219         default:
220           DUNE_THROW(RangeError, "Component out of range.");
221         }
222       } else {
223         DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
224       }
225     }
226 
227     /** \brief Polynomial order of the shape functions
228      */
order() const229     unsigned int order() const
230     {
231       return 2;
232     }
233 
234   };
235 }
236 #endif
237