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>
60 void 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>
69 void 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>
91 void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt,
int order, Points& points)
const
93 #if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
94 std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
96 std::size_t nPoints = numLagrangePoints(gt, order);
100 buildTriangle(nPoints, order, points);
101 else if (gt.isQuadrilateral())
102 buildQuad(nPoints, order, points);
104 DUNE_THROW(Dune::NotImplemented,
105 "Lagrange points not yet implemented for this GeometryType.");
108 assert(points.size() == nPoints);
115 template <
class Po
ints>
116 void LagrangePointSetBuilder<K,2>::buildTriangle (std::size_t nPoints,
int order, Points& points)
const
118 points.reserve(nPoints);
120 const int nVertexDOFs = 3;
121 const int nEdgeDOFs = 3 * (order-1);
123 static const unsigned int vertexPerm[3] = {0,1,2};
124 static const unsigned int edgePerm[3] = {0,2,1};
126 auto calcKey = [&](
int p) -> Key
128 if (p < nVertexDOFs) {
129 return Key{vertexPerm[p], dim, 0};
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};
137 unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
138 return Key{0, dim-2, index};
142 std::array<int,3> bindex;
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)});
157 void LagrangePointSetBuilder<K,2>::barycentricIndex (
int index, std::array<int,3>& bindex,
int order)
163 while (index != 0 && index >= 3 * order)
174 bindex[index] = bindex[(index + 1) % 3] = min;
175 bindex[(index + 2) % 3] = max;
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;
194 template <
class Po
ints>
195 void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints,
int order, Points& points)
const
197 points.resize(nPoints);
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.}};
202 for (
int n = 0; n <= orders[1]; ++n) {
203 for (
int m = 0; m <= orders[0]; ++m) {
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));
211 auto [idx,key] = calcQuadKey(m,n,orders);
212 points[idx] = LP{p, key};
222 std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
223 LagrangePointSetBuilder<K,2>::calcQuadKey (
int i,
int j, std::array<int,2> order)
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);
230 unsigned int entityIndex = 0;
231 unsigned int index = 0;
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});
244 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
245 entityIndex = j ? 3 : 2;
249 dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
250 entityIndex = i ? 1 : 0;
253 return std::make_pair(dof, Key{entityIndex, dim-1, index});
256 offset += 2 * (order[0]-1 + order[1]-1);
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()});
267 template <
class Po
ints>
268 void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt,
unsigned int order, Points& points)
const
270 #if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
271 std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
273 std::size_t nPoints = numLagrangePoints(gt, order);
276 if (gt.isTetrahedron())
277 buildTetra(nPoints, order, points);
278 else if (gt.isHexahedron())
279 buildHex(nPoints, order, points);
281 DUNE_THROW(Dune::NotImplemented,
282 "Lagrange points not yet implemented for this GeometryType.");
285 assert(points.size() == nPoints);
293 template <
class Po
ints>
294 void LagrangePointSetBuilder<K,3>::buildTetra (std::size_t nPoints,
int order, Points& points)
const
296 points.reserve(nPoints);
298 const int nVertexDOFs = 4;
299 const int nEdgeDOFs = 6 * (order-1);
300 const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
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};
306 auto calcKey = [&](
int p) -> Key
308 if (p < nVertexDOFs) {
309 return Key{vertexPerm[p], dim, 0};
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};
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};
322 unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
323 return Key{0, dim-3, index};
327 std::array<int,4> bindex;
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)});
342 void LagrangePointSetBuilder<K,3>::barycentricIndex (
int p, std::array<int,4>& bindex,
int order)
344 const int nVertexDOFs = 4;
345 const int nEdgeDOFs = 6 * (order-1);
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};
357 while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
359 p -= 2 * (order * order + 1);
368 for (
int coord = 0; coord < 4; ++coord)
369 bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
372 else if (p < nVertexDOFs+nEdgeDOFs)
374 int edgeId = (p - nVertexDOFs) / (order-1);
375 int vertexId = (p - nVertexDOFs) % (order-1);
376 for (
int coord = 0; coord < 4; ++coord)
378 bindex[coord] = min +
379 (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
380 linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
386 int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
387 int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
389 std::array<int,3> projectedBIndex;
391 projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
393 LagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
395 for (
int i = 0; i < 3; i++)
396 bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
398 bindex[faceMinCoord[faceId]] = min;
407 template <
class Po
ints>
408 void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints,
int order, Points& points)
const
410 points.resize(nPoints);
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.}};
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));
425 auto [idx,key] = calcHexKey(m,n,p,orders);
426 points[idx] = LP{point, key};
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)
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);
444 unsigned int entityIndex = 0;
445 unsigned int index = 0;
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});
457 entityIndex = (k ? 8 : 4);
460 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
462 entityIndex += (i ? 1 : 0);
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;
468 entityIndex += (j ? 3 : 2);
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;
475 entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
477 return std::make_pair(dof, Key{entityIndex, dim-1, index});
480 offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
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;
491 offset += 2 * (order[1] - 1) * (order[2] - 1);
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;
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;
506 return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
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()});