dune-vtk  0.2
lagrangepoints.impl.hh
Go to the documentation of this file.
1 #pragma once
2 
3 #include <cassert>
4 #include <array>
5 
6 #include <dune/common/exceptions.hh>
7 #include <dune/common/version.hh>
8 #include <dune/geometry/type.hh>
9 #include <dune/localfunctions/lagrange/equidistantpoints.hh>
10 
11 namespace Dune {
12 namespace Vtk {
13 namespace Impl {
14 
58 template <class K>
59  template <class Points>
60 void LagrangePointSetBuilder<K,0>::operator() (GeometryType gt, int /*order*/, Points& points) const
61 {
62  assert(gt.isVertex());
63  points.push_back(LP{Vec{},Key{0,0,0}});
64 }
65 
66 
67 template <class K>
68  template <class Points>
69 void LagrangePointSetBuilder<K,1>::operator() (GeometryType gt, int order, Points& points) const
70 {
71  assert(gt.isLine());
72 
73  // Vertex nodes
74  points.push_back(LP{Vec{0.0}, Key{0,dim,0}});
75  points.push_back(LP{Vec{1.0}, Key{1,dim,0}});
76 
77  if (order > 1) {
78  // Inner nodes
79  Vec p{0.0};
80  for (int i = 0; i < order-1; i++)
81  {
82  p[0] += 1.0 / order;
83  points.push_back(LP{p,Key{0,dim-1,(unsigned int)(i)}});
84  }
85  }
86 }
87 
88 
89 template <class K>
90  template <class Points>
91 void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt, int order, Points& points) const
92 {
93 #if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
94  std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
95 #else
96  std::size_t nPoints = numLagrangePoints(gt, order);
97 #endif
98 
99  if (gt.isTriangle())
100  buildTriangle(nPoints, order, points);
101  else if (gt.isQuadrilateral())
102  buildQuad(nPoints, order, points);
103  else {
104  DUNE_THROW(Dune::NotImplemented,
105  "Lagrange points not yet implemented for this GeometryType.");
106  }
107 
108  assert(points.size() == nPoints);
109 }
110 
111 
112 // Construct the point set in a triangle element.
113 // Loop from the outside to the inside
114 template <class K>
115  template <class Points>
116 void LagrangePointSetBuilder<K,2>::buildTriangle (std::size_t nPoints, int order, Points& points) const
117 {
118  points.reserve(nPoints);
119 
120  const int nVertexDOFs = 3;
121  const int nEdgeDOFs = 3 * (order-1);
122 
123  static const unsigned int vertexPerm[3] = {0,1,2};
124  static const unsigned int edgePerm[3] = {0,2,1};
125 
126  auto calcKey = [&](int p) -> Key
127  {
128  if (p < nVertexDOFs) {
129  return Key{vertexPerm[p], dim, 0};
130  }
131  else if (p < nVertexDOFs+nEdgeDOFs) {
132  unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
133  unsigned int index = (p - nVertexDOFs) % (order-1);
134  return Key{edgePerm[entityIndex], dim-1, index};
135  }
136  else {
137  unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
138  return Key{0, dim-2, index};
139  }
140  };
141 
142  std::array<int,3> bindex;
143 
144  K order_d = K(order);
145  for (std::size_t p = 0; p < nPoints; ++p) {
146  barycentricIndex(p, bindex, order);
147  Vec point{bindex[0] / order_d, bindex[1] / order_d};
148  points.push_back(LP{point, calcKey(p)});
149  }
150 }
151 
152 
153 // "Barycentric index" is a triplet of integers, each running from 0 to
154 // <Order>. It is the index of a point on the triangle in barycentric
155 // coordinates.
156 template <class K>
157 void LagrangePointSetBuilder<K,2>::barycentricIndex (int index, std::array<int,3>& bindex, int order)
158 {
159  int max = order;
160  int min = 0;
161 
162  // scope into the correct triangle
163  while (index != 0 && index >= 3 * order)
164  {
165  index -= 3 * order;
166  max -= 2;
167  min++;
168  order -= 3;
169  }
170 
171  // vertex DOFs
172  if (index < 3)
173  {
174  bindex[index] = bindex[(index + 1) % 3] = min;
175  bindex[(index + 2) % 3] = max;
176  }
177  // edge DOFs
178  else
179  {
180  index -= 3;
181  int dim = index / (order - 1);
182  int offset = (index - dim * (order - 1));
183  bindex[(dim + 1) % 3] = min;
184  bindex[(dim + 2) % 3] = (max - 1) - offset;
185  bindex[dim] = (min + 1) + offset;
186  }
187 }
188 
189 
190 // Construct the point set in the quad element
191 // 1. build equispaced points with index tuple (i,j)
192 // 2. map index tuple to DOF index and LocalKey
193 template <class K>
194  template <class Points>
195 void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints, int order, Points& points) const
196 {
197  points.resize(nPoints);
198 
199  std::array<int,2> orders{order, order};
200  std::array<Vec,4> nodes{Vec{0., 0.}, Vec{1., 0.}, Vec{1., 1.}, Vec{0., 1.}};
201 
202  for (int n = 0; n <= orders[1]; ++n) {
203  for (int m = 0; m <= orders[0]; ++m) {
204  // int idx = pointIndexFromIJ(m,n,orders);
205 
206  const K r = K(m) / orders[0];
207  const K s = K(n) / orders[1];
208  Vec p = K(1.0-r) * (nodes[3] * s + nodes[0] * K(1.0 - s)) +
209  r * (nodes[2] * s + nodes[1] * K(1.0 - s));
210 
211  auto [idx,key] = calcQuadKey(m,n,orders);
212  points[idx] = LP{p, key};
213  // points[idx] = LP{p, calcQuadKey(n,m,orders)};
214  }
215  }
216 }
217 
218 
219 // Obtain the VTK DOF index of the node (i,j) in the quad element
220 // and construct a LocalKey
221 template <class K>
222 std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
223 LagrangePointSetBuilder<K,2>::calcQuadKey (int i, int j, std::array<int,2> order)
224 {
225  const bool ibdy = (i == 0 || i == order[0]);
226  const bool jbdy = (j == 0 || j == order[1]);
227  const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0); // How many boundaries do we lie on at once?
228 
229  int dof = 0;
230  unsigned int entityIndex = 0;
231  unsigned int index = 0;
232 
233  if (nbdy == 2) // Vertex DOF
234  {
235  dof = (i ? (j ? 2 : 1) : (j ? 3 : 0));
236  entityIndex = (j ? (i ? 3 : 2) : (i ? 1 : 0));
237  return std::make_pair(dof,Key{entityIndex, dim, 0});
238  }
239 
240  int offset = 4;
241  if (nbdy == 1) // Edge DOF
242  {
243  if (!ibdy) {
244  dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
245  entityIndex = j ? 3 : 2;
246  index = i-1;
247  }
248  else if (!jbdy) {
249  dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
250  entityIndex = i ? 1 : 0;
251  index = j-1;
252  }
253  return std::make_pair(dof, Key{entityIndex, dim-1, index});
254  }
255 
256  offset += 2 * (order[0]-1 + order[1]-1);
257 
258  // nbdy == 0: Face DOF
259  dof = offset + (i - 1) + (order[0]-1) * ((j - 1));
260  Key innerKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
261  return std::make_pair(dof, Key{0, dim-2, innerKey.index()});
262 }
263 
264 
265 // Lagrange points on 3d geometries
266 template <class K>
267  template <class Points>
268 void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt, unsigned int order, Points& points) const
269 {
270 #if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
271  std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
272 #else
273  std::size_t nPoints = numLagrangePoints(gt, order);
274 #endif
275 
276  if (gt.isTetrahedron())
277  buildTetra(nPoints, order, points);
278  else if (gt.isHexahedron())
279  buildHex(nPoints, order, points);
280  else {
281  DUNE_THROW(Dune::NotImplemented,
282  "Lagrange points not yet implemented for this GeometryType.");
283  }
284 
285  assert(points.size() == nPoints);
286 }
287 
288 
289 // Construct the point set in the tetrahedron element
290 // 1. construct barycentric (index) coordinates
291 // 2. obtains the DOF index, LocalKey and actual coordinate from barycentric index
292 template <class K>
293  template <class Points>
294 void LagrangePointSetBuilder<K,3>::buildTetra (std::size_t nPoints, int order, Points& points) const
295 {
296  points.reserve(nPoints);
297 
298  const int nVertexDOFs = 4;
299  const int nEdgeDOFs = 6 * (order-1);
300  const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
301 
302  static const unsigned int vertexPerm[4] = {0,1,2,3};
303  static const unsigned int edgePerm[6] = {0,2,1,3,4,5};
304  static const unsigned int facePerm[4] = {1,2,0,3};
305 
306  auto calcKey = [&](int p) -> Key
307  {
308  if (p < nVertexDOFs) {
309  return Key{vertexPerm[p], dim, 0};
310  }
311  else if (p < nVertexDOFs+nEdgeDOFs) {
312  unsigned int entityIndex = (p - nVertexDOFs) / (order-1);
313  unsigned int index = (p - nVertexDOFs) % (order-1);
314  return Key{edgePerm[entityIndex], dim-1, index};
315  }
316  else if (p < nVertexDOFs+nEdgeDOFs+nFaceDOFs) {
317  unsigned int index = (p - (nVertexDOFs + nEdgeDOFs)) % ((order-1)*(order-2)/2);
318  unsigned int entityIndex = (p - (nVertexDOFs + nEdgeDOFs)) / ((order-1)*(order-2)/2);
319  return Key{facePerm[entityIndex], dim-2, index};
320  }
321  else {
322  unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
323  return Key{0, dim-3, index};
324  }
325  };
326 
327  std::array<int,4> bindex;
328 
329  K order_d = K(order);
330  for (std::size_t p = 0; p < nPoints; ++p) {
331  barycentricIndex(p, bindex, order);
332  Vec point{bindex[0] / order_d, bindex[1] / order_d, bindex[2] / order_d};
333  points.push_back(LP{point, calcKey(p)});
334  }
335 }
336 
337 
338 // "Barycentric index" is a set of 4 integers, each running from 0 to
339 // <Order>. It is the index of a point in the tetrahedron in barycentric
340 // coordinates.
341 template <class K>
342 void LagrangePointSetBuilder<K,3>::barycentricIndex (int p, std::array<int,4>& bindex, int order)
343 {
344  const int nVertexDOFs = 4;
345  const int nEdgeDOFs = 6 * (order-1);
346 
347  static const int edgeVertices[6][2] = {{0,1}, {1,2}, {2,0}, {0,3}, {1,3}, {2,3}};
348  static const int linearVertices[4][4] = {{0,0,0,1}, {1,0,0,0}, {0,1,0,0}, {0,0,1,0}};
349  static const int vertexMaxCoords[4] = {3,0,1,2};
350  static const int faceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}};
351  static const int faceMinCoord[4] = {1,3,0,2};
352 
353  int max = order;
354  int min = 0;
355 
356  // scope into the correct tetra
357  while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
358  {
359  p -= 2 * (order * order + 1);
360  max -= 3;
361  min++;
362  order -= 4;
363  }
364 
365  // vertex DOFs
366  if (p < nVertexDOFs)
367  {
368  for (int coord = 0; coord < 4; ++coord)
369  bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
370  }
371  // edge DOFs
372  else if (p < nVertexDOFs+nEdgeDOFs)
373  {
374  int edgeId = (p - nVertexDOFs) / (order-1);
375  int vertexId = (p - nVertexDOFs) % (order-1);
376  for (int coord = 0; coord < 4; ++coord)
377  {
378  bindex[coord] = min +
379  (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
380  linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
381  }
382  }
383  // face DOFs
384  else
385  {
386  int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
387  int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
388 
389  std::array<int,3> projectedBIndex;
390  if (order == 3)
391  projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
392  else
393  LagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
394 
395  for (int i = 0; i < 3; i++)
396  bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
397 
398  bindex[faceMinCoord[faceId]] = min;
399  }
400 }
401 
402 
403 // Construct the point set in the hexahedral element
404 // 1. build equispaced points with index tuple (i,j,k)
405 // 2. map index tuple to DOF index and LocalKey
406 template <class K>
407  template <class Points>
408 void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints, int order, Points& points) const
409 {
410  points.resize(nPoints);
411 
412  std::array<int,3> orders{order, order, order};
413  std::array<Vec,8> nodes{Vec{0., 0., 0.}, Vec{1., 0., 0.}, Vec{1., 1., 0.}, Vec{0., 1., 0.},
414  Vec{0., 0., 1.}, Vec{1., 0., 1.}, Vec{1., 1., 1.}, Vec{0., 1., 1.}};
415 
416  for (int p = 0; p <= orders[2]; ++p) {
417  for (int n = 0; n <= orders[1]; ++n) {
418  for (int m = 0; m <= orders[0]; ++m) {
419  const K r = K(m) / orders[0];
420  const K s = K(n) / orders[1];
421  const K t = K(p) / orders[2];
422  Vec point = K(1.0-r) * ((nodes[3] * K(1.0-t) + nodes[7] * t) * s + (nodes[0] * K(1.0-t) + nodes[4] * t) * K(1.0-s)) +
423  r * ((nodes[2] * K(1.0-t) + nodes[6] * t) * s + (nodes[1] * K(1.0-t) + nodes[5] * t) * K(1.0-s));
424 
425  auto [idx,key] = calcHexKey(m,n,p,orders);
426  points[idx] = LP{point, key};
427  }
428  }
429  }
430 }
431 
432 
433 // Obtain the VTK DOF index of the node (i,j,k) in the hexahedral element
434 template <class K>
435 std::pair<int,typename LagrangePointSetBuilder<K,3>::Key>
436 LagrangePointSetBuilder<K,3>::calcHexKey (int i, int j, int k, std::array<int,3> order)
437 {
438  const bool ibdy = (i == 0 || i == order[0]);
439  const bool jbdy = (j == 0 || j == order[1]);
440  const bool kbdy = (k == 0 || k == order[2]);
441  const int nbdy = (ibdy ? 1 : 0) + (jbdy ? 1 : 0) + (kbdy ? 1 : 0); // How many boundaries do we lie on at once?
442 
443  int dof = 0;
444  unsigned int entityIndex = 0;
445  unsigned int index = 0;
446 
447  if (nbdy == 3) // Vertex DOF
448  {
449  dof = (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
450  entityIndex = (i ? 1 : 0) + (j ? 2 : 0) + (k ? 4 : 0);
451  return std::make_pair(dof, Key{entityIndex, dim, 0});
452  }
453 
454  int offset = 8;
455  if (nbdy == 2) // Edge DOF
456  {
457  entityIndex = (k ? 8 : 4);
458  if (!ibdy)
459  { // On i axis
460  dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
461  index = i;
462  entityIndex += (i ? 1 : 0);
463  }
464  else if (!jbdy)
465  { // On j axis
466  dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
467  index = j;
468  entityIndex += (j ? 3 : 2);
469  }
470  else
471  { // !kbdy, On k axis
472  offset += 4 * (order[0]-1) + 4 * (order[1]-1);
473  dof = (k - 1) + (order[2]-1) * (i ? (j ? 3 : 1) : (j ? 2 : 0)) + offset;
474  index = k;
475  entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
476  }
477  return std::make_pair(dof, Key{entityIndex, dim-1, index});
478  }
479 
480  offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
481  if (nbdy == 1) // Face DOF
482  {
483  Key faceKey;
484  if (ibdy) // On i-normal face
485  {
486  dof = (j - 1) + ((order[1]-1) * (k - 1)) + (i ? (order[1]-1) * (order[2]-1) : 0) + offset;
487  entityIndex = (i ? 1 : 0);
488  faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(j-1,k-1,{order[1]-2, order[2]-2}).second;
489  }
490  else {
491  offset += 2 * (order[1] - 1) * (order[2] - 1);
492  if (jbdy) // On j-normal face
493  {
494  dof = (i - 1) + ((order[0]-1) * (k - 1)) + (j ? (order[2]-1) * (order[0]-1) : 0) + offset;
495  entityIndex = (j ? 3 : 2);
496  faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,k-1,{order[0]-2, order[2]-2}).second;
497  }
498  else
499  { // kbdy, On k-normal face
500  offset += 2 * (order[2]-1) * (order[0]-1);
501  dof = (i - 1) + ((order[0]-1) * (j - 1)) + (k ? (order[0]-1) * (order[1]-1) : 0) + offset;
502  entityIndex = (k ? 5 : 4);
503  faceKey = LagrangePointSetBuilder<K,2>::calcQuadKey(i-1,j-1,{order[0]-2, order[1]-2}).second;
504  }
505  }
506  return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
507  }
508 
509  // nbdy == 0: Body DOF
510  offset += 2 * ((order[1]-1) * (order[2]-1) + (order[2]-1) * (order[0]-1) + (order[0]-1) * (order[1]-1));
511  dof = offset + (i - 1) + (order[0]-1) * ((j - 1) + (order[1]-1) * ((k - 1)));
512  Key innerKey = LagrangePointSetBuilder<K,3>::calcHexKey(i-1,j-1,k-1,{order[0]-2, order[1]-2, order[2]-2}).second;
513  return std::make_pair(dof, Key{0, dim-3, innerKey.index()});
514 }
515 
516 }}} // end namespace Dune::Vtk::Impl
Definition: writer.hh:13