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>
59 template <
class Po
ints>
60void LagrangePointSetBuilder<K,0>::operator() (GeometryType gt,
int , Points& points)
const
62 assert(gt.isVertex());
63 points.push_back(LP{Vec{},Key{0,0,0}});
68 template <
class Po
ints>
69void LagrangePointSetBuilder<K,1>::operator() (GeometryType gt,
int order, Points& points)
const
74 points.push_back(LP{Vec{0.0}, Key{0,dim,0}});
75 points.push_back(LP{Vec{1.0}, Key{1,dim,0}});
80 for (
int i = 0; i < order-1; i++)
83 points.push_back(LP{p,Key{0,dim-1,(
unsigned int)(i)}});
90 template <
class Po
ints>
91void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt,
int order, Points& points)
const
93 std::size_t nPoints = numLagrangePoints(gt, order);
96 buildTriangle(nPoints, order, points);
97 else if (gt.isQuadrilateral())
98 buildQuad(nPoints, order, points);
100 DUNE_THROW(Dune::NotImplemented,
101 "Lagrange points not yet implemented for this GeometryType.");
104 assert(points.size() == nPoints);
111 template <
class Po
ints>
112void LagrangePointSetBuilder<K,2>::buildTriangle (std::size_t nPoints,
int order, Points& points)
const
114 points.reserve(nPoints);
116 const int nVertexDOFs = 3;
117 const int nEdgeDOFs = 3 * (order-1);
119 static const unsigned int vertexPerm[3] = {0,1,2};
120 static const unsigned int edgePerm[3] = {0,2,1};
122 auto calcKey = [&](
int p) -> Key
124 if (p < nVertexDOFs) {
125 return Key{vertexPerm[p], dim, 0};
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};
133 unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
134 return Key{0, dim-2, index};
138 std::array<int,3> bindex;
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)});
153void LagrangePointSetBuilder<K,2>::barycentricIndex (
int index, std::array<int,3>& bindex,
int order)
159 while (index != 0 && index >= 3 * order)
170 bindex[index] = bindex[(index + 1) % 3] = min;
171 bindex[(index + 2) % 3] = max;
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;
190 template <
class Po
ints>
191void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints,
int order, Points& points)
const
193 points.resize(nPoints);
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.}};
198 for (
int n = 0; n <= orders[1]; ++n) {
199 for (
int m = 0; m <= orders[0]; ++m) {
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));
207 auto [idx,key] = calcQuadKey(m,n,orders);
208 points[idx] = LP{p, key};
218std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
219LagrangePointSetBuilder<K,2>::calcQuadKey (
int i,
int j, std::array<int,2> order)
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);
226 unsigned int entityIndex = 0;
227 unsigned int index = 0;
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});
240 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
241 entityIndex = j ? 3 : 2;
245 dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
246 entityIndex = i ? 1 : 0;
249 return std::make_pair(dof, Key{entityIndex, dim-1, index});
252 offset += 2 * (order[0]-1 + order[1]-1);
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()});
263 template <
class Po
ints>
264void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt,
unsigned int order, Points& points)
const
266 std::size_t nPoints = numLagrangePoints(gt, order);
268 if (gt.isTetrahedron())
269 buildTetra(nPoints, order, points);
270 else if (gt.isHexahedron())
271 buildHex(nPoints, order, points);
273 DUNE_THROW(Dune::NotImplemented,
274 "Lagrange points not yet implemented for this GeometryType.");
277 assert(points.size() == nPoints);
285 template <
class Po
ints>
286void LagrangePointSetBuilder<K,3>::buildTetra (std::size_t nPoints,
int order, Points& points)
const
288 points.reserve(nPoints);
290 const int nVertexDOFs = 4;
291 const int nEdgeDOFs = 6 * (order-1);
292 const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
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};
298 auto calcKey = [&](
int p) -> Key
300 if (p < nVertexDOFs) {
301 return Key{vertexPerm[p], dim, 0};
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};
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};
314 unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
315 return Key{0, dim-3, index};
319 std::array<int,4> bindex;
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)});
334void LagrangePointSetBuilder<K,3>::barycentricIndex (
int p, std::array<int,4>& bindex,
int order)
336 const int nVertexDOFs = 4;
337 const int nEdgeDOFs = 6 * (order-1);
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};
349 while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
351 p -= 2 * (order * order + 1);
360 for (
int coord = 0; coord < 4; ++coord)
361 bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
364 else if (p < nVertexDOFs+nEdgeDOFs)
366 int edgeId = (p - nVertexDOFs) / (order-1);
367 int vertexId = (p - nVertexDOFs) % (order-1);
368 for (
int coord = 0; coord < 4; ++coord)
370 bindex[coord] = min +
371 (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
372 linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
378 int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
379 int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
381 std::array<int,3> projectedBIndex;
383 projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
385 LagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
387 for (
int i = 0; i < 3; i++)
388 bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
390 bindex[faceMinCoord[faceId]] = min;
399 template <
class Po
ints>
400void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints,
int order, Points& points)
const
402 points.resize(nPoints);
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.}};
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));
417 auto [idx,key] = calcHexKey(m,n,p,orders);
418 points[idx] = LP{point, key};
427std::pair<int,typename LagrangePointSetBuilder<K,3>::Key>
428LagrangePointSetBuilder<K,3>::calcHexKey (
int i,
int j,
int k, std::array<int,3> order)
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);
436 unsigned int entityIndex = 0;
437 unsigned int index = 0;
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});
449 entityIndex = (k ? 8 : 4);
452 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
454 entityIndex += (i ? 1 : 0);
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;
460 entityIndex += (j ? 3 : 2);
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;
467 entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
469 return std::make_pair(dof, Key{entityIndex, dim-1, index});
472 offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
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;
483 offset += 2 * (order[1] - 1) * (order[2] - 1);
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;
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;
498 return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
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()});