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