MiniTensor Version of the Day
Loading...
Searching...
No Matches
MiniTensor_Geometry.h
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MiniTensor Package
4//
5// Copyright 2016 NTESS and the MiniTensor contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#if !defined(MiniTensor_Geometry_h)
11#define MiniTensor_Geometry_h
12
13#include <vector>
14#include "MiniTensor_Tensor.h"
15
16namespace minitensor {
17
20
24namespace ELEMENT{
25
36
37} //namespace ELEMENT
38
42template<typename T, Index N>
44T
45length(Vector<T, N> const & p0, Vector<T, N> const & p1);
46
50template<typename T, Index N>
52T
53area(Vector<T, N> const & p0, Vector<T, N> const & p1,
54 Vector<T, N> const & p2);
55
65template<typename T, Index N>
67T
68area(Vector<T, N> const & p0, Vector<T, N> const & p1,
69 Vector<T, N> const & p2, Vector<T, N> const & p3);
70
74template<typename T, Index N>
76T
77volume(Vector<T, N> const & p0, Vector<T, N> const & p1,
78 Vector<T, N> const & p2, Vector<T, N> const & p3);
79
89template<typename T, Index N>
91T
92volume(Vector<T, N> const & p0, Vector<T, N> const & p1,
93 Vector<T, N> const & p2, Vector<T, N> const & p3,
94 Vector<T, N> const & p4, Vector<T, N> const & p5,
95 Vector<T, N> const & p6, Vector<T, N> const & p7);
96
103template<typename T, Index N>
104Vector<T, N>
105centroid(std::vector<Vector<T, N>> const & points);
106
112template<typename T, Index N>
114Vector<T, N>
115normal(Vector<T, N> const & p0,
116 Vector<T, N> const & p1,
117 Vector<T, N> const & p2);
118
125template<typename T, Index N>
127bool
129 Vector<T, N> const & p,
130 Vector<T, N> const & p0,
131 Vector<T, N> const & p1,
132 Vector<T, N> const & p2,
133 T const tolerance = 0);
134
142template<typename T, typename I, Index N>
144std::pair<Vector<T, N>, Vector<T, N>>
145bounding_box(I start, I end);
146
147template<typename T, typename I>
149std::pair<Vector<T, DYNAMIC>, Vector<T, DYNAMIC>>
150bounding_box(I start, I end);
151
159template<typename T, Index N>
161bool
162in_box(
163 Vector<T, N> const & p,
164 Vector<T, N> const & min,
165 Vector<T, N> const & max);
166
173template<typename T, Index N>
175Vector<T, N>
177 Vector<T, N> const & min,
178 Vector<T, N> const & max);
179
185template<typename T, Index N>
187bool
189 Vector<T, N> const & p,
190 Vector<T, N> const & p0,
191 Vector<T, N> const & p1,
192 Vector<T, N> const & p2,
193 Vector<T, N> const & p3,
194 T const tolerance = 0);
195
202template<typename T, Index N>
204bool
206 Vector<T, N> const & p,
207 Vector<T, N> const & p0,
208 Vector<T, N> const & p1,
209 Vector<T, N> const & p2,
210 Vector<T, N> const & p3,
211 Vector<T, N> const & p4,
212 Vector<T, N> const & p5,
213 Vector<T, N> const & p6,
214 Vector<T, N> const & p7,
215 T const tolerance = 0);
216
223template<typename T, Index N>
224typename std::vector<Vector<T, N>>::size_type
225closest_point(Vector<T, N> const & p, std::vector<Vector<T, N>> const & n);
226
233template<typename T, typename Iterator>
235T
236median(Iterator begin, Iterator end);
237
248template<typename T, Index N>
250Vector<T, N>
252 Vector<T, dimension_const<N, 2>::value> & xi,
253 Vector<T, N> const & p0,
254 Vector<T, N> const & p1,
255 Vector<T, N> const & p2,
256 Vector<T, N> const & p3);
257
267template<typename T, Index N>
269Vector<T, N>
271 Vector<T, dimension_const<N, 3>::value> & xi,
272 Vector<T, N> const & p0,
273 Vector<T, N> const & p1,
274 Vector<T, N> const & p2);
275
290template<typename T, Index N>
292Vector<T, N>
294 Vector<T, dimension_const<N, 3>::value> & xi,
295 Vector<T, N> const & p0,
296 Vector<T, N> const & p1,
297 Vector<T, N> const & p2,
298 Vector<T, N> const & p3,
299 Vector<T, N> const & p4,
300 Vector<T, N> const & p5,
301 Vector<T, N> const & p6,
302 Vector<T, N> const & p7);
303
314template<typename T, Index N>
316Vector<T, N>
318 Vector<T, dimension_const<N, 4>::value> & xi,
319 Vector<T, N> const & p0,
320 Vector<T, N> const & p1,
321 Vector<T, N> const & p2,
322 Vector<T, N> const & p3);
323
332template<typename T, Index M, Index N>
334Vector<T, N>
336 ELEMENT::Type element_type,
337 Vector<T, M> & xi,
338 std::vector<Vector<T, N>> const & v);
339
346template<typename T, Index N>
347std::vector< std::vector<T>>
348distance_matrix(std::vector<Vector<T, N>> const & points);
349
356template<typename T>
357std::vector<T>
358minimum_distances(std::vector< std::vector<T>> const & distances);
359
366find_type(Index const dimension, Index const number_nodes);
367
371template<typename T, Index N>
475
479template<typename T, Index N>
583
587template<typename T, Index N>
690
694template<typename T, Index N>
798
802template<typename T, Index N>
907
911template<typename T, Index N>
913{
914
915public:
916
922
931 Vector<T, N> const & lower,
932 Vector<T, N> const & upper,
933 Vector<Index, N> const & points_per_dimension);
934
939 template<typename Visitor>
941 void
942 traverse(Visitor & visitor) const;
943
944private:
945
951
957
963
964};
965
966} // namespace minitensor
967
968namespace minitensor {
969
970//
971// Helper functions for determining the type of element
972//
973namespace {
974
977find_type_1D(Index const nodes)
978{
979 switch (nodes) {
980 case 2: return ELEMENT::SEGMENTAL;
981 default: return ELEMENT::UNKNOWN;
982 }
983}
984
987find_type_2D(Index const nodes)
988{
989 switch (nodes) {
990 case 3: return ELEMENT::TRIANGULAR;
991 case 4: return ELEMENT::QUADRILATERAL;
992 default: return ELEMENT::UNKNOWN;
993 }
994}
995
998find_type_3D(Index const nodes)
999{
1000 switch (nodes) {
1001 case 4: return ELEMENT::TETRAHEDRAL;
1002 case 8: return ELEMENT::HEXAHEDRAL;
1003 default: return ELEMENT::UNKNOWN;
1004 }
1005}
1006
1007} // anonymous namespace
1008
1009
1010//
1011//
1012//
1015find_type(Index const dimension, Index const number_nodes)
1016{
1017
1019 type = ELEMENT::UNKNOWN;
1020
1021 switch (dimension) {
1022
1023 case 1:
1024 type = find_type_1D(number_nodes);
1025 break;
1026
1027 case 2:
1028 type = find_type_2D(number_nodes);
1029 break;
1030
1031 case 3:
1032 type = find_type_3D(number_nodes);
1033 break;
1034
1035 default:
1036 break;
1037
1038 }
1039
1040 if (type == ELEMENT::UNKNOWN) {
1041 MT_ERROR_EXIT("Unknown element type.");
1042 }
1043
1044 return type;
1045}
1046
1047//
1048// Constructor for SphericalParametrization
1049//
1050template<typename T, Index N>
1053 Tensor4<T, N> const & A) : tangent_(A)
1054{
1055#if defined(KOKKOS_ENABLE_CUDA)
1056 minimum_ = DBL_MAX;
1057 maximum_ = DBL_MIN;
1058#else
1059 minimum_ = std::numeric_limits<T>::max();
1060 maximum_ = std::numeric_limits<T>::min();
1061#endif
1062 return;
1063}
1064
1065//
1066// Normal vector for SphericalParametrization
1067//
1068template<typename T, Index N>
1072 Vector<T, dimension_const<N, 2>::value> const & parameters
1073) const
1074{
1075 T const &
1076 phi = parameters(0);
1077
1078 T const &
1079 theta = parameters(1);
1080
1082 normal(sin(phi) * sin(theta), cos(phi), sin(phi) * cos(theta));
1083
1084 return normal;
1085}
1086
1087//
1088// Evaluation for SphericalParametrization
1089//
1090template<typename T, Index N>
1092void
1094 Vector<T, dimension_const<N, 2>::value> const & parameters
1095)
1096{
1097 assert(parameters.get_dimension() == 2);
1098
1099 Vector<T, N> const
1100 normal = get_normal(parameters);
1101
1102 // Localization tensor
1103 Tensor<T, N> const
1104 Q = dot2(normal, dot(tangent_, normal));
1105
1106 T const
1107 determinant = det(Q);
1108
1109 if (determinant < minimum_) {
1110 minimum_ = determinant;
1111 arg_minimum_ = parameters;
1112 }
1113
1114 if (determinant > maximum_) {
1115 maximum_ = determinant;
1116 arg_maximum_ = parameters;
1117 }
1118
1119 return;
1120}
1121
1122//
1123// Constructor for StereographicParametrization
1124//
1125template<typename T, Index N>
1128 Tensor4<T, N> const & A) : tangent_(A)
1129{
1130#if defined(KOKKOS_ENABLE_CUDA)
1131 minimum_ = DBL_MAX;
1132 maximum_ = DBL_MIN;
1133#else
1134 minimum_ = std::numeric_limits<T>::max();
1135 maximum_ = std::numeric_limits<T>::min();
1136#endif
1137 return;
1138}
1139
1140//
1141// Normal vector for StereographicParametrization
1142//
1143template<typename T, Index N>
1147 Vector<T, dimension_const<N, 2>::value> const & parameters
1148) const
1149{
1150 T const &
1151 x = parameters(0);
1152
1153 T const &
1154 y = parameters(1);
1155
1156 T const
1157 r2 = x * x + y * y;
1158
1160 normal(2.0 * x, 2.0 * y, r2 - 1.0);
1161
1162 normal /= (r2 + 1.0);
1163
1164 return normal;
1165}
1166
1167//
1168// Evaluation for StereographicParametrization
1169//
1170template<typename T, Index N>
1172void
1174 Vector<T, dimension_const<N, 2>::value> const & parameters
1175)
1176{
1177 assert(parameters.get_dimension() == 2);
1178
1179 Vector<T, N> const
1180 normal = get_normal(parameters);
1181
1182 // Localization tensor
1183 Tensor<T, N> const
1184 Q = dot2(normal, dot(tangent_, normal));
1185
1186 T const
1187 determinant = det(Q);
1188
1189 if (determinant < minimum_) {
1190 minimum_ = determinant;
1191 arg_minimum_ = parameters;
1192 }
1193
1194 if (determinant > maximum_) {
1195 maximum_ = determinant;
1196 arg_maximum_ = parameters;
1197 }
1198
1199 return;
1200}
1201
1202//
1203// Constructor for ProjectiveParametrization
1204//
1205template<typename T, Index N>
1208 Tensor4<T, N> const & A) : tangent_(A)
1209{
1210#if defined(KOKKOS_ENABLE_CUDA)
1211 minimum_ = DBL_MAX;
1212 maximum_ = DBL_MIN;
1213#else
1214 minimum_ = std::numeric_limits<T>::max();
1215 maximum_ = std::numeric_limits<T>::min();
1216#endif
1217 return;
1218}
1219
1220//
1221// Normal vector for ProjectiveParametrization
1222//
1223template<typename T, Index N>
1227 Vector<T, dimension_const<N, 3>::value> const & parameters
1228) const
1229{
1230 T const &
1231 x = parameters(0);
1232
1233 T const &
1234 y = parameters(1);
1235
1236 T const &
1237 z = parameters(2);
1238
1240 normal(x, y, z);
1241
1242 T const
1243 n = norm(normal);
1244
1245 if (n > 0.0) {
1246 normal /= n;
1247 } else {
1248 normal = Vector<T, N>(1.0, 1.0, 1.0);
1249 }
1250
1251 return normal;
1252}
1253
1254//
1255// Evaluation for ProjectiveParametrization
1256//
1257template<typename T, Index N>
1259void
1261 Vector<T, dimension_const<N, 3>::value> const & parameters
1262)
1263{
1264 assert(parameters.get_dimension() == 3);
1265
1266 Vector<T, N> const
1267 normal = get_normal(parameters);
1268
1269 // Localization tensor
1270 Tensor<T, N> const
1271 Q = dot2(normal, dot(tangent_, normal));
1272
1273 T const
1274 determinant = det(Q);
1275
1276 if (determinant < minimum_) {
1277 minimum_ = determinant;
1278 arg_minimum_ = parameters;
1279 }
1280
1281 if (determinant > maximum_) {
1282 maximum_ = determinant;
1283 arg_maximum_ = parameters;
1284 }
1285
1286 return;
1287}
1288
1289//
1290// Constructor for TangentParametrization
1291//
1292template<typename T, Index N>
1295 Tensor4<T, N> const & A) : tangent_(A)
1296{
1297#if defined(KOKKOS_ENABLE_CUDA)
1298 minimum_ = DBL_MAX;
1299 maximum_ = DBL_MIN;
1300#else
1301 minimum_ = std::numeric_limits<T>::max();
1302 maximum_ = std::numeric_limits<T>::min();
1303#endif
1304 return;
1305}
1306
1307//
1308// Normal vector for TangentParametrization
1309//
1310template<typename T, Index N>
1314 Vector<T, dimension_const<N, 2>::value> const & parameters
1315) const
1316{
1317 T const &
1318 x = parameters(0);
1319
1320 T const &
1321 y = parameters(1);
1322
1323 T const
1324 r = std::sqrt(x * x + y * y);
1325
1328
1329 if (r > 0.0) {
1330 normal(0) = x * sin(r) / r;
1331 normal(1) = y * sin(r) / r;
1332 normal(2) = cos(r);
1333 } else {
1334 normal(2) = 1.0;
1335 }
1336
1337 return normal;
1338}
1339
1340//
1341// Evaluation for TangentParametrization
1342//
1343template<typename T, Index N>
1345void
1347 Vector<T, dimension_const<N, 2>::value> const & parameters
1348)
1349{
1350 assert(parameters.get_dimension() == 2);
1351
1352 Vector<T, N> const
1353 normal = get_normal(parameters);
1354
1355 // Localization tensor
1356 Tensor<T, N> const
1357 Q = dot2(normal, dot(tangent_, normal));
1358
1359 T const
1360 determinant = det(Q);
1361
1362 if (determinant < minimum_) {
1363 minimum_ = determinant;
1364 arg_minimum_ = parameters;
1365 }
1366
1367 if (determinant > maximum_) {
1368 maximum_ = determinant;
1369 arg_maximum_ = parameters;
1370 }
1371
1372 return;
1373}
1374
1375//
1376// Constructor for CartesianParametrization
1377//
1378template<typename T, Index N>
1381 Tensor4<T, N> const & A) : tangent_(A)
1382{
1383#if defined(KOKKOS_ENABLE_CUDA)
1384 minimum_ = DBL_MAX;
1385 maximum_ = DBL_MIN;
1386#else
1387 minimum_ = std::numeric_limits<T>::max();
1388 maximum_ = std::numeric_limits<T>::min();
1389 return;
1390#endif
1391}
1392
1393//
1394// Normal vector for CartesianParametrization
1395//
1396template<typename T, Index N>
1400 Vector<T, dimension_const<N, 3>::value> const & parameters
1401) const
1402{
1403 T const &
1404 x = parameters(0);
1405
1406 T const &
1407 y = parameters(1);
1408
1409 T const
1410 z = parameters(2);
1411
1412 Vector<T, N> const
1413 normal(x, y, z);
1414
1415 return normal;
1416}
1417
1418//
1419// Evaluation for CartesianParametrization
1420//
1421template<typename T, Index N>
1423void
1425 Vector<T, dimension_const<N, 3>::value> const & parameters
1426)
1427{
1429 normal = get_normal(parameters);
1430
1431 // Localization tensor
1432 Tensor<T, N> const
1433 Q = dot2(normal, dot(tangent_, normal));
1434
1435 T const
1436 determinant = det(Q);
1437
1438 if (determinant < minimum_) {
1439 minimum_ = determinant;
1440 arg_minimum_ = parameters;
1441 }
1442
1443 if (determinant > maximum_) {
1444 maximum_ = determinant;
1445 arg_maximum_ = parameters;
1446 }
1447
1448 return;
1449}
1450
1451//
1452// Constructor for ParametricGrid
1453//
1454template<typename T, Index N>
1457 Vector<T, N> const & lower,
1458 Vector<T, N> const & upper,
1459 Vector<Index, N> const & points_per_dimension)
1460{
1461 assert(lower.get_dimension() == upper.get_dimension());
1462 assert(lower.get_dimension() == points_per_dimension.get_dimension());
1463
1464 lower_ = lower;
1465 upper_ = upper;
1466 points_per_dimension_ = points_per_dimension;
1467
1468 return;
1469}
1470
1471//
1472// Traverse the grid and apply the visitor to each point.
1473//
1474template<typename T, Index N>
1475template<typename Visitor>
1477void
1478ParametricGrid<T, N>::traverse(Visitor & visitor) const
1479{
1480 // Loop over the grid
1481 Index const
1482 number_parameters = lower_.get_dimension();
1483
1484 LongIndex
1485 total_number_points = 1;
1486
1487 for (Index dimension = 0; dimension < number_parameters; ++dimension) {
1488 total_number_points *= points_per_dimension_(dimension);
1489 }
1490
1491 Vector<LongIndex, N> steps(number_parameters, Filler::ONES);
1492
1493 for (Index dimension = 1; dimension < number_parameters; ++dimension) {
1494 steps(dimension) =
1495 steps(dimension - 1) * points_per_dimension_(dimension - 1);
1496 }
1497
1499 indices(number_parameters, Filler::ZEROS);
1500
1502 position_in_grid(number_parameters, Filler::ZEROS);
1503
1504 Vector<T, N> const
1505 span = upper_ - lower_;
1506
1507 for (LongIndex point = 1; point <= total_number_points; ++point) {
1508
1509 //std::cout << "Indices : ";
1510
1511 for (Index dimension = 0; dimension < number_parameters; ++dimension) {
1512
1513/* if ( points_per_dimension_(dimension) == 1 ) {
1514
1515 position_in_grid(dimension) = lower_(dimension);
1516
1517 } else {
1518
1519 position_in_grid(dimension) = indices(dimension) * span(dimension) /
1520 (points_per_dimension_(dimension) - 1) + lower_(dimension);
1521 }
1522*/
1523 position_in_grid(dimension) = indices(dimension) * span(dimension) /
1524 (points_per_dimension_(dimension) - 1) + lower_(dimension);
1525 visitor(position_in_grid);
1526
1527 //std::cout << indices(dimension) << " ";
1528
1529 // Check if index needs to be increased or rolled back
1530 if (point % steps(dimension) == 0) {
1531 ++indices(dimension);
1532 }
1533 if (indices(dimension) == points_per_dimension_(dimension)) {
1534 indices(dimension) = 0;
1535 }
1536
1537 }
1538 visitor(position_in_grid);
1539 //std::cout << std::endl;
1540 //std::cout << "Position : " << position_in_grid << std::endl;
1541
1542 }
1543
1544 return;
1545}
1546
1547} // namespace minitensor
1548#include <iterator>
1549
1550namespace minitensor {
1551
1552//
1553// Length of a segment
1554//
1555template<typename T, Index N>
1557T
1558length(Vector<T, N> const & p0, Vector<T, N> const & p1)
1559{
1560 return norm(p1 - p0);
1561}
1562
1563//
1564// Area of a triangle
1565//
1566template<typename T, Index N>
1568T
1569area(Vector<T, N> const & p0, Vector<T, N> const & p1,
1570 Vector<T, N> const & p2)
1571{
1572 Vector<T, N> const u = p1 - p0;
1573 Vector<T, N> const v = p2 - p0;
1574
1575 T const area = 0.5 * norm(cross(u, v));
1576
1577 return area;
1578}
1579
1580//
1581// Area of a quadrilateral.
1582// Taken from:
1583// Calculation of the volume of a general hexahedron for flow predictions
1584// Davies, D. E.; Salmond, D. J.
1585// AIAA Journal (ISSN 0001-1452), vol. 23, June 1985, p. 954-956.
1586// Their vertex naming convention: ABCD
1587// Our convention: 0123
1588//
1589template<typename T, Index N>
1591T
1593 Vector<T, N> const & p0, Vector<T, N> const & p1,
1594 Vector<T, N> const & p2, Vector<T, N> const & p3)
1595{
1596 Vector<T, N> const v31 = p1 - p3;
1597 Vector<T, N> const v02 = p2 - p0;
1598
1599 T const area = 0.5 * norm(cross(v31, v02));
1600
1601 return area;
1602}
1603
1604//
1605// Volume of tetrahedron
1606//
1607template<typename T, Index N>
1609T
1610volume(Vector<T, N> const & p0, Vector<T, N> const & p1,
1611 Vector<T, N> const & p2, Vector<T, N> const & p3)
1612{
1613 Vector<T, N> const u = p1 - p0;
1614 Vector<T, N> const v = p2 - p0;
1615 Vector<T, N> const w = p3 - p0;
1616
1617 T const volume = std::abs(dot(u, cross(v, w)))/ 6.0;
1618
1619 return volume;
1620}
1621
1622//
1623// Volume of hexahedron
1624// Taken from:
1625// Calculation of the volume of a general hexahedron for flow predictions
1626// Davies, D. E.; Salmond, D. J.
1627// AIAA Journal (ISSN 0001-1452), vol. 23, June 1985, p. 954-956.
1628// Their vertex naming convention: ABCDEFGH
1629// Our convention: 45670123
1630//
1631template<typename T, Index N>
1633T
1635 Vector<T, N> const & p0, Vector<T, N> const & p1,
1636 Vector<T, N> const & p2, Vector<T, N> const & p3,
1637 Vector<T, N> const & p4, Vector<T, N> const & p5,
1638 Vector<T, N> const & p6, Vector<T, N> const & p7)
1639{
1640 Vector<T, N> const v24 = p4 - p2;
1641 Vector<T, N> const v25 = p5 - p2;
1642 Vector<T, N> const v20 = p0 - p2;
1643 Vector<T, N> const v27 = p7 - p2;
1644
1645 Vector<T, N> const v75 = p5 - p7;
1646 Vector<T, N> const v46 = p6 - p4;
1647 Vector<T, N> const v50 = p0 - p5;
1648 Vector<T, N> const v41 = p1 - p4;
1649 Vector<T, N> const v07 = p7 - p0;
1650 Vector<T, N> const v43 = p3 - p4;
1651
1652 Vector<T, N> const v26 = p6 - p2;
1653 Vector<T, N> const v16 = p6 - p1;
1654 Vector<T, N> const v21 = p1 - p2;
1655 Vector<T, N> const v31 = p1 - p3;
1656 Vector<T, N> const v23 = p3 - p2;
1657 Vector<T, N> const v63 = p3 - p6;
1658
1659 Vector<T, N> const v7546 = cross(v75, v46);
1660 Vector<T, N> const v5041 = cross(v50, v41);
1661 Vector<T, N> const v0743 = cross(v07, v43);
1662 Vector<T, N> const v2616 = cross(v26, v16);
1663 Vector<T, N> const v2131 = cross(v21, v31);
1664 Vector<T, N> const v2363 = cross(v23, v63);
1665
1666 T const V1 = dot(v24, v7546 + v5041 + v0743);
1667 T const V2 = dot(v25, v7546 + v2616);
1668 T const V3 = dot(v20, v5041 + v2131);
1669 T const V4 = dot(v27, v0743 + v2363);
1670
1671 T const volume = (V1 + V2 + V3 + V4) / 12.0;
1672
1673 return volume;
1674}
1675
1676//
1677// Centroids of segment, triangle, tetrahedron, quadrilateral
1678// and hexahedron
1679// For these we can just take the average of the vertices.
1680// WARNING: This is not the center of mass.
1681//
1682template<typename T, Index N>
1683Vector<T, N>
1684centroid(std::vector<Vector<T, N>> const & points)
1685{
1686 Vector<T, N> C(points[0].get_dimension());
1687 C.clear();
1688 typedef typename std::vector<Vector<T, N>>::size_type sizeT;
1689 sizeT const n = points.size();
1690
1691 for (sizeT i = 0; i < n; ++i) {
1692 C += points[i];
1693 }
1694 return C / static_cast<T>(n);
1695}
1696
1697//
1698// The surface normal of a face
1699// Input: 3 independent nodes on the face
1700// Output: unit normal vector
1701//
1702template<typename T, Index N>
1704Vector<T, N>
1706 Vector<T, N> const & p1,
1707 Vector<T, N> const & p2)
1708{
1709 // Construct 2 independent vectors
1710 Vector<T, N> const v0 = p1 - p0;
1711 Vector<T, N> const v1 = p2 - p0;
1712
1713 Vector<T, N> const n = unit(cross(v0, v1));
1714
1715 return n;
1716}
1717
1718//
1719// Given 3 points p0, p1, p2 that define a plane
1720// determine if point p is in the same side of the normal
1721// to the plane as defined by the right hand rule.
1722//
1723template<typename T, Index N>
1725bool
1727 Vector<T, N> const & p,
1728 Vector<T, N> const & p0,
1729 Vector<T, N> const & p1,
1730 Vector<T, N> const & p2,
1731 T const tolerance)
1732{
1733 Vector<T, N> const v0 = p1 - p0;
1734 Vector<T, N> const v1 = p2 - p0;
1735 T const h = std::min(norm(v0), norm(v1));
1736 Vector<T, N> const n = unit(cross(v0, v1));
1737 Vector<T, N> const v = p - p0;
1738
1739 T const s = dot(v, n);
1740
1741 if (s < -tolerance * h) return false;
1742
1743 return true;
1744}
1745
1746//
1747// Given two iterators to a container of points,
1748// find the associated bounding box.
1749// \param start, end: define sequence of points
1750// \return vectors that define the bounding box
1751//
1752template<typename T, typename I, Index N>
1754std::pair<Vector<T, N>, Vector<T, N>>
1755bounding_box(I start, I end)
1756{
1757 I
1758 it = start;
1759
1761 min = (*it);
1762
1764 max = min;
1765
1766 Index const
1767 dimension = min.get_dimension();
1768
1769 ++it;
1770
1771 for (; it != end; ++it) {
1772
1773 Vector<T, N> const &
1774 point = (*it);
1775
1776 for (Index i = 0; i < dimension; ++i) {
1777 T const s = point(i);
1778 if (s < min(i)) min(i) = s;
1779 if (s > max(i)) max(i) = s;
1780 }
1781
1782 }
1783
1784 return std::make_pair(min, max);
1785}
1786
1790template<typename T, typename I>
1792std::pair<Vector<T, DYNAMIC>, Vector<T, DYNAMIC>>
1793bounding_box(I start, I end)
1794{
1795 return bounding_box<T, I, DYNAMIC>(start, end);
1796}
1797
1798//
1799// Determine if a given point is inside a bounding box.
1800// \param p the point
1801// \param min, max points defining the box
1802// \return whether the point is inside
1803//
1804template<typename T, Index N>
1806bool
1808 Vector<T, N> const & p,
1809 Vector<T, N> const & min,
1810 Vector<T, N> const & max)
1811{
1812 Index const
1813 dimension = p.get_dimension();
1814
1815 assert(min.get_dimension() == dimension);
1816 assert(max.get_dimension() == dimension);
1817
1818 for (Index i = 0; i < dimension; ++i) {
1819 T const & s = p(i);
1820 if (s < min(i)) return false;
1821 if (s > max(i)) return false;
1822 }
1823
1824 return true;
1825}
1826
1827//
1828// Generate random point inside bounding box
1829// \param min, max the bounding box
1830// \return p point inside box
1831//
1832template<typename T, Index N>
1834Vector<T, N>
1836{
1837 Index const
1838 dimension = min.get_dimension();
1839
1840 assert(max.get_dimension() == dimension);
1841
1842 Vector<T, N> p(dimension);
1843
1844 for (Index i = 0; i < dimension; ++i) {
1845 p(i) = (max(i) - min(i)) * T(std::rand())/T(RAND_MAX) + min(i);
1846 }
1847
1848 return p;
1849}
1850
1851//
1852// Given 4 points p0, p1, p2, p3 that define a tetrahedron
1853// determine if point p is inside it.
1854//
1855template<typename T, Index N>
1857bool
1859 Vector<T, N> const & p,
1860 Vector<T, N> const & p0,
1861 Vector<T, N> const & p1,
1862 Vector<T, N> const & p2,
1863 Vector<T, N> const & p3,
1864 T const tolerance)
1865{
1866 if (in_normal_side(p, p0, p1, p2, tolerance) == false) return false;
1867 if (in_normal_side(p, p0, p3, p1, tolerance) == false) return false;
1868 if (in_normal_side(p, p1, p3, p2, tolerance) == false) return false;
1869 if (in_normal_side(p, p2, p3, p0, tolerance) == false) return false;
1870
1871 return true;
1872}
1873
1874//
1875// Given 8 points that define a hexahedron
1876// determine if point p is inside it.
1877// Assumption: faces are planar
1878//
1879template<typename T, Index N>
1881bool
1883 Vector<T, N> const & p,
1884 Vector<T, N> const & p0,
1885 Vector<T, N> const & p1,
1886 Vector<T, N> const & p2,
1887 Vector<T, N> const & p3,
1888 Vector<T, N> const & p4,
1889 Vector<T, N> const & p5,
1890 Vector<T, N> const & p6,
1891 Vector<T, N> const & p7,
1892 T const tolerance)
1893{
1894 if (in_normal_side(p, p0, p1, p2, tolerance) == false) return false;
1895 if (in_normal_side(p, p0, p4, p5, tolerance) == false) return false;
1896 if (in_normal_side(p, p1, p5, p6, tolerance) == false) return false;
1897 if (in_normal_side(p, p2, p6, p7, tolerance) == false) return false;
1898 if (in_normal_side(p, p3, p7, p4, tolerance) == false) return false;
1899 if (in_normal_side(p, p4, p7, p6, tolerance) == false) return false;
1900
1901 return true;
1902}
1903
1904//
1905// Closest point
1906// \param p the point
1907// \param n vector of points to test
1908// \return index to closest point
1909//
1910template<typename T, Index N>
1911typename std::vector<Vector<T, N>>::size_type
1912closest_point(Vector<T, N> const & p, std::vector<Vector<T, N>> const & n)
1913{
1914 assert(n.size() > 0);
1915
1916 typename std::vector<Vector<T, N>>::size_type
1917 index = 0;
1918
1919 Vector<T, N> const
1920 v0 = p - n[0];
1921
1922 T
1923 min = norm_square(v0);
1924
1925 for (typename std::vector<Vector<T, N>>::size_type i = 1;
1926 i < n.size();
1927 ++i) {
1928
1929 Vector<T, N> const
1930 vi = p - n[i];
1931
1932 T const
1933 s = norm_square(vi);
1934
1935 if (s < min) {
1936 min = s;
1937 index = i;
1938 }
1939
1940 }
1941
1942 return index;
1943}
1944
1945// Median of a sequence defined by random
1946// access iterators. Undefined for empty set.
1947// \param begin, end Iterators that define the sequence
1948// \return median of sequence
1949//
1950template<typename T, typename Iterator>
1952T
1953median(Iterator begin, Iterator end)
1954{
1955 // Firewall
1956 if (begin == end) {
1957 MT_ERROR_EXIT("Median undefined for empty set.");
1958 }
1959
1960 Index const
1961 size = static_cast<Index>(std::distance(begin, end));
1962
1963 T
1964 median;
1965
1966 Index const
1967 mid_index = size / 2;
1968
1969 Iterator
1970 mid_iterator = begin + mid_index;
1971 std::partial_sort(begin, mid_iterator, end);
1972
1973 if (size % 2 == 0) {
1974
1975 // Even number of elements
1976 T const
1977 b = *mid_iterator;
1978
1979 Iterator
1980 previous = mid_iterator - 1;
1981
1982 T const
1983 a = *previous;
1984
1985 median = (a + b) / 2.0;
1986
1987 } else {
1988
1989 // Odd number of elements
1990 median = *mid_iterator;
1991
1992 }
1993
1994 return median;
1995}
1996
1997//
1998// Given quadrilateral nodes and a position
1999// in parametric coordinates, interpolate.
2000// \param xi position in parametric coordinates
2001// \param p0 ... corner nodes
2002// \return interpolated position
2003//
2004template<typename T, Index N>
2006Vector<T, N>
2009 Vector<T, N> const & p0,
2010 Vector<T, N> const & p1,
2011 Vector<T, N> const & p2,
2012 Vector<T, N> const & p3)
2013{
2014
2015 T const
2016 N0 = 0.25 * (1 - xi(0)) * (1 - xi(1));
2017
2018 T const
2019 N1 = 0.25 * (1 + xi(0)) * (1 - xi(1));
2020
2021 T const
2022 N2 = 0.25 * (1 + xi(0)) * (1 + xi(1));
2023
2024 T const
2025 N3 = 0.25 * (1 - xi(0)) * (1 + xi(1));
2026
2027 Vector<T, N> const
2028 p = N0 * p0 + N1 * p1 + N2 * p2 + N3 * p3;
2029
2030 return p;
2031}
2032
2033//
2034// Given triangle nodes and a position
2035// in parametric coordinates, interpolate.
2036// \param xi position in parametric coordinates
2037// \param p0 ... corner nodes
2038// \return interpolated position
2039//
2040template<typename T, Index N>
2042Vector<T, N>
2045 Vector<T, N> const & p0,
2046 Vector<T, N> const & p1,
2047 Vector<T, N> const & p2)
2048{
2049 xi(2) = 1.0 - xi(0) - xi(1);
2050
2051 Vector<T, N> const
2052 p = xi(0) * p0 + xi(1) * p1 + xi(2) * p2;
2053
2054 return p;
2055}
2056
2057//
2058// Given hexahedron nodes and a position
2059// in parametric coordinates, interpolate.
2060// \param xi position in parametric coordinates
2061// \param p0 ... corner nodes
2062// \return interpolated position
2063//
2064template<typename T, Index N>
2066Vector<T, N>
2069 Vector<T, N> const & p0,
2070 Vector<T, N> const & p1,
2071 Vector<T, N> const & p2,
2072 Vector<T, N> const & p3,
2073 Vector<T, N> const & p4,
2074 Vector<T, N> const & p5,
2075 Vector<T, N> const & p6,
2076 Vector<T, N> const & p7)
2077{
2078
2079 T const
2080 N0 = 0.125 * (1 - xi(0)) * (1 - xi(1)) * (1 - xi(2));
2081
2082 T const
2083 N1 = 0.125 * (1 + xi(0)) * (1 - xi(1)) * (1 - xi(2));
2084
2085 T const
2086 N2 = 0.125 * (1 + xi(0)) * (1 + xi(1)) * (1 - xi(2));
2087
2088 T const
2089 N3 = 0.125 * (1 - xi(0)) * (1 + xi(1)) * (1 - xi(2));
2090
2091 T const
2092 N4 = 0.125 * (1 - xi(0)) * (1 - xi(1)) * (1 + xi(2));
2093
2094 T const
2095 N5 = 0.125 * (1 + xi(0)) * (1 - xi(1)) * (1 + xi(2));
2096
2097 T const
2098 N6 = 0.125 * (1 + xi(0)) * (1 + xi(1)) * (1 + xi(2));
2099
2100 T const
2101 N7 = 0.125 * (1 - xi(0)) * (1 + xi(1)) * (1 + xi(2));
2102
2103 Vector<T, N> const
2104 p =
2105 N0 * p0 + N1 * p1 + N2 * p2 + N3 * p3 +
2106 N4 * p4 + N5 * p5 + N6 * p6 + N7 * p7;
2107
2108 return p;
2109}
2110
2111//
2112// Given tetrahedron nodes and a position
2113// in parametric coordinates, interpolate.
2114// \param xi position in parametric coordinates
2115// \param p0 ... corner nodes
2116// \return interpolated position
2117//
2118template<typename T, Index N>
2120Vector<T, N>
2123 Vector<T, N> const & p0,
2124 Vector<T, N> const & p1,
2125 Vector<T, N> const & p2,
2126 Vector<T, N> const & p3)
2127{
2128 xi(3) = 1.0 - xi(0) - xi(1) - xi(2);
2129
2130 Vector<T, N> const
2131 p = xi(0) * p0 + xi(1) * p1 + xi(2) * p2 + xi(3) * p3;
2132
2133 return p;
2134}
2135
2136//
2137// Given element type and nodes and a position
2138// in parametric coordinates, interpolate.
2139// \param type element type
2140// \param xi position in parametric coordinates
2141// \param v ... corner nodes
2142// \return interpolated position
2143//
2144template<typename T, Index M, Index N>
2146Vector<T, N>
2148 ELEMENT::Type element_type,
2149 Vector<T, M> &xi,
2150 std::vector<Vector<T, N>> const &v)
2151{
2152 Vector<T, N> p;
2153
2154 switch (element_type) {
2155
2157 p = interpolate_triangle(xi, v[0], v[1], v[2]);
2158 break;
2159
2161 p = interpolate_quadrilateral(xi, v[0], v[1], v[2], v[3]);
2162 break;
2163
2165 p = interpolate_tetrahedron(xi, v[0], v[1], v[2], v[3]);
2166 break;
2167
2170 xi, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
2171 break;
2172
2173 default:
2174 MT_ERROR_EXIT("Unknown element type in interpolation.");
2175 break;
2176
2177 }
2178
2179 return p;
2180}
2181
2182//
2183// Given a vector of points, determine
2184// distances between all of them.
2185// \param vector of points
2186// \return distance matrix
2187//
2188template<typename T, Index N>
2189std::vector< std::vector<T>>
2190distance_matrix(std::vector<Vector<T, N>> const & points)
2191{
2192 Index const
2193 number_points = points.size();
2194
2195 std::vector< std::vector<T>>
2196 distances(number_points);
2197
2198 for (Index i = 0; i < number_points; ++i) {
2199
2200 distances[i].resize(number_points);
2201
2202 distances[i][i] = 0.0;
2203
2204 for (Index j = i + 1; j < number_points; ++j) {
2205
2206 T const
2207 distance = norm(points[i] - points[j]);
2208
2209 distances[i][j] = distance;
2210 distances[j][i] = distance;
2211
2212 }
2213
2214 }
2215
2216 return distances;
2217}
2218
2219//
2220// Given a distance matrix, determine the minimum
2221// distance between two distinct points.
2222// \param distance matrix
2223// \return minimum distance
2224//
2225template<typename T>
2226std::vector<T>
2227minimum_distances(std::vector< std::vector<T>> const & distances)
2228{
2229 Index const
2230 number_points = distances.size();
2231
2232 std::vector<T>
2233 minima(number_points);
2234
2235 // First row
2236 T
2237 minimum = distances[0][1];
2238
2239 for (Index j = 2; j < number_points; ++j) {
2240 minimum = std::min(minimum, distances[0][j]);
2241 }
2242
2243 minima[0] = minimum;
2244
2245 // Remaining rows
2246 for (Index i = 1; i < number_points; ++i) {
2247
2248 minimum = distances[i][0];
2249
2250 for (Index j = 1; j < number_points; ++j) {
2251
2252 if (i == j) continue;
2253
2254 minimum = std::min(minimum, distances[i][j]);
2255
2256 }
2257
2258 minima[i] = minimum;
2259
2260 }
2261
2262 return minima;
2263}
2264
2266} // namespace minitensor
2267
2268#endif // MiniTensor_Geometry_h
#define KOKKOS_INLINE_FUNCTION
#define MT_ERROR_EXIT(...)
KOKKOS_INLINE_FUNCTION void clear()
static constexpr Index value
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, N > cross(Vector< S, N > const &u, Vector< T, N > const &v)
KOKKOS_INLINE_FUNCTION Index get_dimension() const
KOKKOS_INLINE_FUNCTION Vector< typename Promote< S, T >::type, M > dot(Matrix< T, M, N > const &A, Vector< S, N > const &u)
KOKKOS_INLINE_FUNCTION T norm_square(Vector< T, N > const &u)
KOKKOS_INLINE_FUNCTION Tensor< typename Promote< S, T >::type, N > dot2(Tensor3< T, N > const &A, Vector< S > const &u)
KOKKOS_INLINE_FUNCTION Vector< T, N > normal(Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, 3 > get_arg_minimum() const
KOKKOS_INLINE_FUNCTION void operator()(Vector< T, dimension_const< N, 3 >::value > const &parameters)
KOKKOS_INLINE_FUNCTION Vector< T, N > interpolate_element(ELEMENT::Type element_type, Vector< T, M > &xi, std::vector< Vector< T, N > > const &v)
std::vector< T > minimum_distances(std::vector< std::vector< T > > const &distances)
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_maximum() const
KOKKOS_INLINE_FUNCTION bool in_normal_side(Vector< T, N > const &p, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, T const tolerance=0)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal(Vector< T, dimension_const< N, 2 >::value > const &parameters) const
KOKKOS_INLINE_FUNCTION bool in_hexahedron(Vector< T, N > const &p, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3, Vector< T, N > const &p4, Vector< T, N > const &p5, Vector< T, N > const &p6, Vector< T, N > const &p7, T const tolerance=0)
KOKKOS_INLINE_FUNCTION T get_minimum() const
KOKKOS_INLINE_FUNCTION void operator()(Vector< T, dimension_const< N, 2 >::value > const &parameters)
KOKKOS_INLINE_FUNCTION bool in_box(Vector< T, N > const &p, Vector< T, N > const &min, Vector< T, N > const &max)
KOKKOS_INLINE_FUNCTION CartesianParametrization(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal(Vector< T, dimension_const< N, 2 >::value > const &parameters) const
KOKKOS_INLINE_FUNCTION T get_maximum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_maximum() const
Vector< Index, N > points_per_dimension_
KOKKOS_INLINE_FUNCTION std::pair< Vector< T, N >, Vector< T, N > > bounding_box(I start, I end)
KOKKOS_INLINE_FUNCTION void operator()(Vector< T, dimension_const< N, 2 >::value > const &parameters)
KOKKOS_INLINE_FUNCTION Vector< T, N > random_in_box(Vector< T, N > const &min, Vector< T, N > const &max)
KOKKOS_INLINE_FUNCTION bool in_tetrahedron(Vector< T, N > const &p, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3, T const tolerance=0)
std::vector< std::vector< T > > distance_matrix(std::vector< Vector< T, N > > const &points)
KOKKOS_INLINE_FUNCTION T get_maximum() const
KOKKOS_INLINE_FUNCTION SphericalParametrization(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION TangentParametrization(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION Vector< T, 3 > get_arg_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > interpolate_tetrahedron(Vector< T, dimension_const< N, 4 >::value > &xi, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3)
KOKKOS_INLINE_FUNCTION Vector< T, N > interpolate_quadrilateral(Vector< T, dimension_const< N, 2 >::value > &xi, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3)
KOKKOS_INLINE_FUNCTION StereographicParametrization(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION ParametricGrid()
KOKKOS_INLINE_FUNCTION void operator()(Vector< T, dimension_const< N, 2 >::value > const &parameters)
KOKKOS_INLINE_FUNCTION Vector< T, 3 > get_arg_maximum() const
std::vector< Vector< T, N > >::size_type closest_point(Vector< T, N > const &p, std::vector< Vector< T, N > > const &n)
KOKKOS_INLINE_FUNCTION T area(Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_maximum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_maximum() const
KOKKOS_INLINE_FUNCTION Vector< T, 3 > get_arg_maximum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_maximum() const
KOKKOS_INLINE_FUNCTION ProjectiveParametrization(Tensor4< T, N > const &A)
KOKKOS_INLINE_FUNCTION T get_maximum() const
KOKKOS_INLINE_FUNCTION T get_maximum() const
KOKKOS_INLINE_FUNCTION T length(Vector< T, N > const &p0, Vector< T, N > const &p1)
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_minimum() const
KOKKOS_INLINE_FUNCTION T get_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_maximum() const
KOKKOS_INLINE_FUNCTION T volume(Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal(Vector< T, dimension_const< N, 3 >::value > const &parameters) const
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal(Vector< T, dimension_const< N, 3 >::value > const &parameters) const
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_maximum() const
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_minimum() const
KOKKOS_INLINE_FUNCTION Vector< T, N > interpolate_hexahedron(Vector< T, dimension_const< N, 3 >::value > &xi, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2, Vector< T, N > const &p3, Vector< T, N > const &p4, Vector< T, N > const &p5, Vector< T, N > const &p6, Vector< T, N > const &p7)
Vector< T, N > centroid(std::vector< Vector< T, N > > const &points)
KOKKOS_INLINE_FUNCTION ELEMENT::Type find_type(Index const dimension, Index const number_nodes)
KOKKOS_INLINE_FUNCTION Vector< T, 2 > get_arg_maximum() const
KOKKOS_INLINE_FUNCTION T get_minimum() const
KOKKOS_INLINE_FUNCTION T get_minimum() const
KOKKOS_INLINE_FUNCTION void traverse(Visitor &visitor) const
KOKKOS_INLINE_FUNCTION Vector< T, N > interpolate_triangle(Vector< T, dimension_const< N, 3 >::value > &xi, Vector< T, N > const &p0, Vector< T, N > const &p1, Vector< T, N > const &p2)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal_minimum() const
KOKKOS_INLINE_FUNCTION T median(Iterator begin, Iterator end)
KOKKOS_INLINE_FUNCTION T get_maximum() const
KOKKOS_INLINE_FUNCTION T get_minimum() const
KOKKOS_INLINE_FUNCTION void operator()(Vector< T, dimension_const< N, 3 >::value > const &parameters)
KOKKOS_INLINE_FUNCTION Vector< T, N > get_normal(Vector< T, dimension_const< N, 2 >::value > const &parameters) const
KOKKOS_INLINE_FUNCTION T norm(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION T det(Tensor< T, N > const &A)
KOKKOS_INLINE_FUNCTION Quaternion< T > unit(Quaternion< T > const &q)
uint32_t Index
Indexing type.
KOKKOS_INLINE_FUNCTION T min(T const &a, T const &b)
KOKKOS_INLINE_FUNCTION T max(T const &a, T const &b)
uint64_t LongIndex
High count type.
Set to constant value if not dynamic.