DUNE PDELab (2.8)

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#ifndef DUNE_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH
4#define DUNE_HIERARCHICAL_PRISM_P2_LOCALBASIS_HH
5
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
17namespace Dune
18{
19 template<class D, class R>
20 class HierarchicalPrismP2LocalBasis
21 {
22 public:
24 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>, Dune::FieldMatrix<R,1,3> > Traits;
25
27 unsigned int size () const
28 {
29 return 18;
30 }
31
33 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
65 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
147 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
229 unsigned int order() const
230 {
231 return 2;
232 }
233
234 };
235}
236#endif
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:69
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:95
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
Dune
Dune namespace.
Definition: alignedallocator.hh:11
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)