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