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