dune-localfunctions  2.4.1
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 
16 namespace 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  // lower edge (2)
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
std::size_t size() const
number of coefficients
Definition: qklocalcoefficients.hh:287
const LocalKey & localKey(std::size_t i) const
get i&#39;th index
Definition: qklocalcoefficients.hh:293
Describe position of one degree of freedom.
Definition: localkey.hh:21
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
QkLocalCoefficients()
Default constructor.
Definition: qklocalcoefficients.hh:223
Attaches a shape function to an entity.
Definition: qklocalcoefficients.hh:24