6 #include <dune/common/exceptions.hh>
7 #include <dune/geometry/type.hh>
8 #include <dune/localfunctions/lagrange/equidistantpoints.hh>
58 template <
class Po
ints>
59 void LagrangePointSetBuilder<K,0>::operator() (GeometryType gt,
int , Points& points)
const
61 assert(gt.isVertex());
62 points.push_back(LP{Vec{},Key{0,0,0}});
67 template <
class Po
ints>
68 void LagrangePointSetBuilder<K,1>::operator() (GeometryType gt,
int order, Points& points)
const
73 points.push_back(LP{Vec{0.0}, Key{0,dim,0}});
74 points.push_back(LP{Vec{1.0}, Key{1,dim,0}});
79 for (
unsigned int i = 0; i < order-1; i++)
82 points.push_back(LP{p,Key{0,dim-1,i}});
89 template <
class Po
ints>
90 void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt,
int order, Points& points)
const
92 std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
95 buildTriangle(nPoints, order, points);
96 else if (gt.isQuadrilateral())
97 buildQuad(nPoints, order, points);
99 DUNE_THROW(Dune::NotImplemented,
100 "Lagrange points not yet implemented for this GeometryType.");
103 assert(points.size() == nPoints);
110 template <
class Po
ints>
111 void LagrangePointSetBuilder<K,2>::buildTriangle (std::size_t nPoints,
int order, Points& points)
const
113 points.reserve(nPoints);
115 const int nVertexDOFs = 3;
116 const int nEdgeDOFs = 3 * (order-1);
118 static const unsigned int vertexPerm[3] = {0,1,2};
119 static const unsigned int edgePerm[3] = {0,2,1};
121 auto calcKey = [&](
int p) -> Key
123 if (p < nVertexDOFs) {
124 return Key{vertexPerm[p], dim, 0};
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};
132 unsigned int index = p - (nVertexDOFs + nEdgeDOFs);
133 return Key{0, dim-2, index};
137 std::array<int,3> bindex;
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)});
152 void LagrangePointSetBuilder<K,2>::barycentricIndex (
int index, std::array<int,3>& bindex,
int order)
158 while (index != 0 && index >= 3 * order)
169 bindex[index] = bindex[(index + 1) % 3] = min;
170 bindex[(index + 2) % 3] = max;
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;
189 template <
class Po
ints>
190 void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints,
int order, Points& points)
const
192 points.resize(nPoints);
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.}};
197 for (
int n = 0; n <= orders[1]; ++n) {
198 for (
int m = 0; m <= orders[0]; ++m) {
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));
206 auto [idx,key] = calcQuadKey(m,n,orders);
207 points[idx] = LP{p, key};
217 std::pair<int,typename LagrangePointSetBuilder<K,2>::Key>
218 LagrangePointSetBuilder<K,2>::calcQuadKey (
int i,
int j, std::array<int,2> order)
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);
225 unsigned int entityIndex = 0;
226 unsigned int index = 0;
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});
239 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + offset;
240 entityIndex = j ? 3 : 2;
244 dof = (j - 1) + (i ? order[0]-1 : 2 * (order[0]-1) + order[1]-1) + offset;
245 entityIndex = i ? 1 : 0;
248 return std::make_pair(dof, Key{entityIndex, dim-1, index});
251 offset += 2 * (order[0]-1 + order[1]-1);
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()});
262 template <
class Po
ints>
263 void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt,
unsigned int order, Points& points)
const
265 std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
267 if (gt.isTetrahedron())
268 buildTetra(nPoints, order, points);
269 else if (gt.isHexahedron())
270 buildHex(nPoints, order, points);
272 DUNE_THROW(Dune::NotImplemented,
273 "Lagrange points not yet implemented for this GeometryType.");
276 assert(points.size() == nPoints);
284 template <
class Po
ints>
285 void LagrangePointSetBuilder<K,3>::buildTetra (std::size_t nPoints,
int order, Points& points)
const
287 points.reserve(nPoints);
289 const int nVertexDOFs = 4;
290 const int nEdgeDOFs = 6 * (order-1);
291 const int nFaceDOFs = 4 * (order-1)*(order-2)/2;
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};
297 auto calcKey = [&](
int p) -> Key
299 if (p < nVertexDOFs) {
300 return Key{vertexPerm[p], dim, 0};
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};
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};
313 unsigned int index = p - (nVertexDOFs + nEdgeDOFs + nFaceDOFs);
314 return Key{0, dim-3, index};
318 std::array<int,4> bindex;
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)});
333 void LagrangePointSetBuilder<K,3>::barycentricIndex (std::size_t p, std::array<int,4>& bindex,
int order)
335 const int nVertexDOFs = 4;
336 const int nEdgeDOFs = 6 * (order-1);
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};
348 while (p >= 2 * (order * order + 1) && p != 0 && order > 3)
350 p -= 2 * (order * order + 1);
359 for (
int coord = 0; coord < 4; ++coord)
360 bindex[coord] = (coord == vertexMaxCoords[p] ? max : min);
363 else if (p < nVertexDOFs+nEdgeDOFs)
365 int edgeId = (p - nVertexDOFs) / (order-1);
366 int vertexId = (p - nVertexDOFs) % (order-1);
367 for (
int coord = 0; coord < 4; ++coord)
369 bindex[coord] = min +
370 (linearVertices[edgeVertices[edgeId][0]][coord] * (max - min - 1 - vertexId) +
371 linearVertices[edgeVertices[edgeId][1]][coord] * (1 + vertexId));
377 int faceId = (p - (nVertexDOFs+nEdgeDOFs)) / ((order-2)*(order-1)/2);
378 int vertexId = (p - (nVertexDOFs+nEdgeDOFs)) % ((order-2)*(order-1)/2);
380 std::array<int,3> projectedBIndex;
382 projectedBIndex[0] = projectedBIndex[1] = projectedBIndex[2] = 0;
384 LagrangePointSetBuilder<K,2>::barycentricIndex(vertexId, projectedBIndex, order-3);
386 for (
int i = 0; i < 3; i++)
387 bindex[faceBCoords[faceId][i]] = (min + 1 + projectedBIndex[i]);
389 bindex[faceMinCoord[faceId]] = min;
398 template <
class Po
ints>
399 void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints,
int order, Points& points)
const
401 points.resize(nPoints);
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.}};
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));
416 auto [idx,key] = calcHexKey(m,n,p,orders);
417 points[idx] = LP{point, key};
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)
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);
435 unsigned int entityIndex = 0;
436 unsigned int index = 0;
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});
448 entityIndex = (k ? 8 : 4);
451 dof = (i - 1) + (j ? order[0]-1 + order[1]-1 : 0) + (k ? 2 * (order[0]-1 + order[1]-1) : 0) + offset;
453 entityIndex += (i ? 1 : 0);
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;
459 entityIndex += (j ? 3 : 2);
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;
466 entityIndex = (i ? 1 : 0) + (j ? 2 : 0);
468 return std::make_pair(dof, Key{entityIndex, dim-1, index});
471 offset += 4 * (order[0]-1 + order[1]-1 + order[2]-1);
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;
482 offset += 2 * (order[1] - 1) * (order[2] - 1);
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;
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;
497 return std::make_pair(dof, Key{entityIndex, dim-2, faceKey.index()});
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()});
Definition: datacollectorinterface.hh:9