dune-localfunctions 3.0-git
qklocalcoefficients.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
4#ifndef DUNE_LOCALFUNCTIONS_QK_LOCALCOEFFICIENTS_HH
5#define DUNE_LOCALFUNCTIONS_QK_LOCALCOEFFICIENTS_HH
6
7#include <array>
8#include <cassert>
9#include <vector>
10
11#include <dune/common/exceptions.hh>
12#include <dune/common/power.hh>
13
15
16namespace Dune
17{
23 template<int k, int d>
25
26 // Return i as a d-digit number in the (k+1)-nary system
27 static std::array<unsigned int,d> multiindex (unsigned int i)
28 {
29 std::array<unsigned int,d> alpha;
30 for (int j=0; j<d; j++)
31 {
32 alpha[j] = i % (k+1);
33 i = i/(k+1);
34 }
35 return alpha;
36 }
37
39 void setup1d(std::vector<unsigned int>& subEntity)
40 {
41 // Special-handling for piecewise constant elements
42 if (k==0)
43 {
44 subEntity[0] = 0;
45 return;
46 }
47
48 unsigned lastIndex=0;
49
50 /* edge and vertex numbering
51 0----0----1
52 */
53
54 // edge (0)
55 subEntity[lastIndex++] = 0; // corner 0
56 for (unsigned i = 0; i < k - 1; ++i)
57 subEntity[lastIndex++] = 0; // inner dofs of element (0)
58
59 subEntity[lastIndex++] = 1; // corner 1
60
61 assert((StaticPower<k+1,d>::power==lastIndex));
62 }
63
64 void setup2d(std::vector<unsigned int>& subEntity)
65 {
66 // Special-handling for piecewise constant elements
67 if (k==0)
68 {
69 subEntity[0] = 0;
70 return;
71 }
72
73 unsigned lastIndex=0;
74
75 // LocalKey: entity number, entity codim, dof indices within each entity
76 /* edge and vertex numbering
77 2----3----3
78 | |
79 | |
80 0 1
81 | |
82 | |
83 0----2----1
84 */
85
86 // lower edge (2)
87 subEntity[lastIndex++] = 0; // corner 0
88 for (unsigned i = 0; i < k - 1; ++i)
89 subEntity[lastIndex++] = 2; // inner dofs of lower edge (2)
90
91 subEntity[lastIndex++] = 1; // corner 1
92
93 // iterate from bottom to top over inner edge dofs
94 for (unsigned e = 0; e < k - 1; ++e) {
95 subEntity[lastIndex++] = 0; // left edge (0)
96 for (unsigned i = 0; i < k - 1; ++i)
97 subEntity[lastIndex++] = 0; // face dofs
98 subEntity[lastIndex++] = 1; // right edge (1)
99 }
100
101 // upper edge (3)
102 subEntity[lastIndex++] = 2; // corner 2
103 for (unsigned i = 0; i < k - 1; ++i)
104 subEntity[lastIndex++] = 3; // inner dofs of upper edge (3)
105
106 subEntity[lastIndex++] = 3; // corner 3
107
108 assert((StaticPower<k+1,d>::power==lastIndex));
109 }
110
111
112
113 void setup3d(std::vector<unsigned int>& subEntity)
114 {
115 // Special-handling for piecewise constant elements
116 if (k==0)
117 {
118 subEntity[0] = 0;
119 return;
120 }
121
122 unsigned lastIndex=0;
123#ifndef NDEBUG
124 const unsigned numIndices = StaticPower<k+1,d>::power;
125 const unsigned numFaceIndices = StaticPower<k+1,d-1>::power;
126#endif
127 const unsigned numInnerEdgeDofs = k-1;
128
129 // LocalKey: entity number, entity codim, dof indices within each entity
130 /* edge and vertex numbering
131
132 6---(11)--7 6---------7
133 /| /| /| (5) /|
134 (8)| (9)| / | top / |
135 / (2) / (3) / |(3)bac/k |
136 4---(10)--5 | 4---------5 |
137 | | | | left|(0)| |(1)|right
138 | 2--(7)|---3 | 2-----|---3
139 (0) / (1) / |(2)front | /
140 |(4) |(5) | / (4) | /
141 |/ |/ |/ bottom |/
142 0---(6)---1 0---------1
143 */
144
145 // bottom face (4)
146 lastIndex=0;
147 // lower edge (6)
148 subEntity[lastIndex++] = 0; // corner 0
149 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
150 subEntity[lastIndex++] = 6; // inner dofs of lower edge (6)
151
152 subEntity[lastIndex++] = 1; // corner 1
153
154 // iterate from bottom to top over inner edge dofs
155 for (unsigned e = 0; e < numInnerEdgeDofs; ++e) {
156 subEntity[lastIndex++] = 4; // left edge (4)
157 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
158 subEntity[lastIndex++] = 4; // inner face dofs
159 subEntity[lastIndex++] = 5; // right edge (5)
160 }
161
162 // upper edge (7)
163 subEntity[lastIndex++] = 2; // corner 2
164 for (unsigned i = 0; i < k - 1; ++i)
165 subEntity[lastIndex++] = 7; // inner dofs of upper edge (7)
166 subEntity[lastIndex++] = 3; // corner 3
167
168 assert(numFaceIndices==lastIndex); // added 1 face so far
170
172 for(unsigned f = 0; f < numInnerEdgeDofs; ++f) {
173
174 // lower edge (connecting edges 0 and 1)
175 subEntity[lastIndex++] = 0; // dof on edge 0
176 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
177 subEntity[lastIndex++] = 2; // dof in front face
178 subEntity[lastIndex++] = 1; // dof on edge 1
179
180 // iterate from bottom to top over inner edge dofs
181 for (unsigned e = 0; e < numInnerEdgeDofs; ++e) {
182 subEntity[lastIndex++] = 0; // on left face (0)
183 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
184 subEntity[lastIndex++] = 0; // volume dofs
185 subEntity[lastIndex++] = 1; // right face (1)
186 }
187
188 // upper edge (connecting edges 0 and 1)
189 subEntity[lastIndex++] = 2; // dof on edge 2
190 for (unsigned i = 0; i < numInnerEdgeDofs; ++i)
191 subEntity[lastIndex++] = 3; // dof on rear face (3)
192 subEntity[lastIndex++] = 3; // dof on edge 3
193
194 assert(lastIndex==(f+1+1)*numFaceIndices);
195 }
196
198 // lower edge (10)
199 subEntity[lastIndex++] = 4; // corner 4
200 for (unsigned i = 0; i < k - 1; ++i)
201 subEntity[lastIndex++] = 10; // inner dofs on lower edge (10)
202 subEntity[lastIndex++] = 5; // corner 5
203
204 // iterate from bottom to top over inner edge dofs
205 for (unsigned e = 0; e < k - 1; ++e) {
206 subEntity[lastIndex++] = 8; // left edge (8)
207 for (unsigned i = 0; i < k - 1; ++i)
208 subEntity[lastIndex++] = 5; // face dofs
209 subEntity[lastIndex++] = 9; // right edge (9)
210 }
211
212 // upper edge (11)
213 subEntity[lastIndex++] = 6; // corner 6
214 for (unsigned i = 0; i < k - 1; ++i)
215 subEntity[lastIndex++] = 11; // inner dofs of upper edge (11)
216 subEntity[lastIndex++] = 7; // corner 7
217
218 assert(numIndices==lastIndex);
219 }
220
221 public:
223 QkLocalCoefficients () : li(StaticPower<k+1,d>::power)
224 {
225 // Set up array of codimension-per-dof-number
226 std::vector<unsigned int> codim(li.size());
227
228 for (std::size_t i=0; i<codim.size(); i++) {
229 codim[i] = 0;
230 if (k==0)
231 continue;
232 // Codimension gets increased by 1 for each coordinate direction
233 // where dof is on boundary
234 std::array<unsigned int,d> mIdx = multiindex(i);
235 for (int j=0; j<d; j++)
236 if (mIdx[j]==0 or mIdx[j]==k)
237 codim[i]++;
238 }
239
240 // Set up index vector (the index of the dof in the set of dofs of a given subentity)
241 // Algorithm: the 'index' has the same ordering as the dof number 'i'.
242 // To make it consecutive we interpret 'i' in the (k+1)-adic system, omit all digits
243 // that correspond to axes where the dof is on the element boundary, and transform the
244 // rest to the (k-1)-adic system.
245 std::vector<unsigned int> index(size());
246
247 for (std::size_t i=0; i<size(); i++) {
248
249 index[i] = 0;
250
251 std::array<unsigned int,d> mIdx = multiindex(i);
252
253 for (int j=d-1; j>=0; j--)
254 if (mIdx[j]>0 and mIdx[j]<k)
255 index[i] = (k-1)*index[i] + (mIdx[j]-1);
256
257 }
258
259 // Set up entity and dof numbers for each (supported) dimension separately
260 std::vector<unsigned int> subEntity(li.size());
261
262 if (k==1) { // We can handle the first-order case in any dimension
263
264 for (std::size_t i=0; i<size(); i++)
265 subEntity[i] = i;
266
267 } else if (d==1) {
268
269 setup1d(subEntity);
270
271 } else if (d==2) {
272
273 setup2d(subEntity);
274
275 } else if (d==3) {
276
277 setup3d(subEntity);
278
279 } else
280 DUNE_THROW(Dune::NotImplemented, "QkLocalCoefficients for k==" << k << " and d==" << d);
281
282 for (size_t i=0; i<li.size(); i++)
283 li[i] = LocalKey(subEntity[i], codim[i], index[i]);
284 }
285
287 std::size_t size () const
288 {
289 return StaticPower<k+1,d>::power;
290 }
291
293 const LocalKey& localKey (std::size_t i) const
294 {
295 return li[i];
296 }
297
298 private:
299 std::vector<LocalKey> li;
300 };
301
302}
303
304#endif
Definition brezzidouglasmarini1cube2d.hh:14
Describe position of one degree of freedom.
Definition localkey.hh:21
Attaches a shape function to an entity.
Definition qklocalcoefficients.hh:24
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition qklocalcoefficients.hh:293
QkLocalCoefficients()
Default constructor.
Definition qklocalcoefficients.hh:223
std::size_t size() const
number of coefficients
Definition qklocalcoefficients.hh:287