diff --git a/cpp/dolfinx/io/cells.cpp b/cpp/dolfinx/io/cells.cpp index 8c142e5e7af..b83b0e1b06f 100644 --- a/cpp/dolfinx/io/cells.cpp +++ b/cpp/dolfinx/io/cells.cpp @@ -33,20 +33,17 @@ int cell_degree(mesh::CellType type, int num_nodes) return n - 1; } case mesh::CellType::tetrahedron: - // n * (n + 1) * (n + 2) = 6*A - // return n - 1; - switch (num_nodes) + { + int n = 0; + while (n * (n + 1) * (n + 2) < 6 * num_nodes) + ++n; + if (n * (n + 1) * (n + 2) != 6 * num_nodes) { - case 4: - return 1; - case 10: - return 2; - case 20: - return 3; - default: throw std::runtime_error("Unknown tetrahedron layout. Number of nodes: " + std::to_string(num_nodes)); } + return n - 1; + } case mesh::CellType::quadrilateral: { const int n = std::sqrt(num_nodes); @@ -102,10 +99,44 @@ std::uint16_t vec_pop(std::vector& v, int i) return value; } //----------------------------------------------------------------------------- +std::vector +vtk_triangle_remainders(std::vector remainders) +{ + std::vector map(remainders.size()); + std::size_t j = 0; + int degree; + while (remainders.size() > 0) + { + if (remainders.size() == 1) + { + map[j++] = vec_pop(remainders, 0); + break; + } + + degree = cell_degree(mesh::CellType::triangle, remainders.size()); + + map[j++] = vec_pop(remainders, 0); + map[j++] = vec_pop(remainders, degree - 1); + map[j++] = vec_pop(remainders, -1); + + for (int i = 0; i < degree - 1; ++i) + map[j++] = vec_pop(remainders, 0); + + for (int i = 1, k = degree * (degree - 1) / 2; i < degree; + k -= degree - i, ++i) + map[j++] = vec_pop(remainders, -k); + + for (int i = 1, k = 1; i < degree; k += i, ++i) + map[j++] = vec_pop(remainders, -k); + } + + return map; +} +//----------------------------------------------------------------------------- std::vector vtk_triangle(int num_nodes) { - // Vertices std::vector map(num_nodes); + // Vertices std::iota(map.begin(), map.begin() + 3, 0); int j = 3; @@ -123,51 +154,239 @@ std::vector vtk_triangle(int num_nodes) // Interior VTK is ordered as a lower order triangle, while FEniCS // orders them lexicographically. std::vector remainders(num_nodes - j); - std::iota(remainders.begin(), remainders.end(), 0); - const std::uint16_t base = 3 * degree; + std::iota(remainders.begin(), remainders.end(), 3 * degree); + + for (std::uint16_t r : vtk_triangle_remainders(remainders)) + { + map[j++] = r; + } + return map; +} +//----------------------------------------------------------------------------- +std::vector +vtk_tetrahedron_remainders(std::vector remainders) +{ + std::vector map(remainders.size()); + std::size_t j = 0; while (remainders.size() > 0) { if (remainders.size() == 1) { - map[j++] = base + vec_pop(remainders, 0); + map[j++] = vec_pop(remainders, 0); break; } - degree = cell_degree(mesh::CellType::triangle, remainders.size()); + const int deg + = cell_degree(mesh::CellType::tetrahedron, remainders.size()) + 1; - map[j++] = base + vec_pop(remainders, 0); - map[j++] = base + vec_pop(remainders, degree - 1); - map[j++] = base + vec_pop(remainders, -1); + map[j++] = vec_pop(remainders, 0); + map[j++] = vec_pop(remainders, deg - 2); + map[j++] = vec_pop(remainders, deg * (deg + 1) / 2 - 3); + map[j++] = vec_pop(remainders, -1); - for (int i = 0; i < degree - 1; ++i) - map[j++] = base + vec_pop(remainders, 0); + if (deg > 2) + { + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, 0); + } + int d = deg - 2; + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, d); + d += deg - 3 - i; + } + d = (deg - 2) * (deg - 1) / 2 - 1; + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, d); + d -= 2 + i; + } + d = (deg - 3) * (deg - 2) / 2; + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, d); + d += (deg - i) * (deg - i - 1) / 2 - 1; + } + d = (deg - 3) * (deg - 2) / 2 + deg - 3; + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, d); + d += (deg - 2 - i) * (deg - 1 - i) / 2 + deg - 4 - i; + } + d = (deg - 3) * (deg - 2) / 2 + deg - 3 + (deg - 2) * (deg - 1) / 2 - 1; + for (int i = 0; i < deg - 2; ++i) + { + map[j++] = vec_pop(remainders, d); + d += (deg - 3 - i) * (deg - 2 - i) / 2 + deg - i - 5; + } + } + if (deg > 3) + { + std::vector dofs((deg - 3) * (deg - 2) / 2); + int di = 0; + int d = (deg - 3) * (deg - 2) / 2; + for (int i = 0; i < deg - 3; ++i) + { + for (int ii = 0; ii < deg - 3 - i; ++ii) + dofs[di++] = vec_pop(remainders, d); + d += (deg - 2 - i) * (deg - 1 - i) / 2 - 1; + } + for (std::uint16_t r : vtk_triangle_remainders(dofs)) + map[j++] = r; - for (int i = 1, k = degree * (degree - 1) / 2; i < degree; - k -= degree - i, ++i) - map[j++] = base + vec_pop(remainders, -k); + di = 0; + int start = deg * deg - 4 * deg + 2; + int sub_i_start = deg - 3; + for (int i = 0; i < deg - 3; ++i) + { + d = start; + int sub_i = sub_i_start; + for (int ii = 0; ii < deg - 3 - i; ++ii) + { + dofs[di++] = vec_pop(remainders, d); + d += sub_i * (sub_i + 1) / 2 - 2 - i; + sub_i -= 1; + } + start -= 2 + i; + } + for (std::uint16_t r : vtk_triangle_remainders(dofs)) + map[j++] = r; - for (int i = 1, k = 1; i < degree; k += i, ++i) - map[j++] = base + vec_pop(remainders, -k); - } + di = 0; + start = (deg - 3) * (deg - 2) / 2; + sub_i_start = deg - 3; + for (int i = 0; i < deg - 3; ++i) + { + d = start; + int sub_i = sub_i_start; + for (int ii = 0; ii < deg - 3 - i; ++ii) + { + dofs[di++] = vec_pop(remainders, d); + d += sub_i * (sub_i + 1) / 2 - 1 - 2 * i; + sub_i -= 1; + } + start += deg - 4 - i; + } + for (std::uint16_t r : vtk_triangle_remainders(dofs)) + map[j++] = r; + di = 0; + int add_start = deg - 4; + for (int i = 0; i < deg - 3; ++i) + { + d = 0; + int add = add_start; + for (int ii = 0; ii < deg - 3 - i; ++ii) + { + dofs[di++] = vec_pop(remainders, d); + d += add; + add -= 1; + } + add_start -= 1; + } + for (std::uint16_t r : vtk_triangle_remainders(dofs)) + map[j++] = r; + } + } return map; } //----------------------------------------------------------------------------- std::vector vtk_tetrahedron(int num_nodes) { - switch (num_nodes) + const int degree = cell_degree(mesh::CellType::tetrahedron, num_nodes); + + std::vector map(num_nodes); + // Vertices + std::iota(map.begin(), map.begin() + 4, 0); + + if (degree < 2) + return map; + + int base = 4; + int j = 4; + const int edge_dofs = degree - 1; + for (int edge : std::vector({5, 2, 4, 3, 1, 0})) { - case 4: - return {0, 1, 2, 3}; - case 10: - return {0, 1, 2, 3, 9, 6, 8, 7, 5, 4}; - case 20: - return {0, 1, 2, 3, 14, 15, 8, 9, 13, 12, - 10, 11, 6, 7, 4, 5, 18, 16, 17, 19}; - default: - throw std::runtime_error("Unknown tetrahedron layout"); + if (edge == 4) + { + for (int i = 0; i < edge_dofs; ++i) + { + map[j++] = base + edge_dofs * (edge + 1) - 1 - i; + } + } + else + { + for (int i = 0; i < edge_dofs; ++i) + { + map[j++] = base + edge_dofs * edge + i; + } + } } + base += 6 * edge_dofs; + + if (degree < 3) + return map; + + const int n_face_dofs = (degree - 1) * (degree - 2) / 2; + + for (int face : std::vector({2, 0, 1, 3})) + { + std::vector face_dofs(n_face_dofs); + std::size_t fj = 0; + if (face == 2) + { + for (int i = 0; i < n_face_dofs; ++i) + { + face_dofs[fj++] = base + n_face_dofs * face + i; + } + } + else if (face == 0) + { + for (int i = degree - 3; i >= 0; --i) + { + int d = i; + for (int ii = 0; ii <= i; ++ii) + { + face_dofs[fj++] = base + n_face_dofs * face + d; + d += degree - 3 - ii; + } + } + } + else + { + for (int i = 0; i < degree - 2; ++i) + { + int d = i; + for (int ii = 0; ii < degree - 2 - i; ++ii) + { + face_dofs[fj++] = base + n_face_dofs * face + d; + d += degree - 2 - ii; + } + } + } + for (std::uint16_t r : vtk_triangle_remainders(face_dofs)) + { + map[j++] = r; + } + } + + base += 4 * n_face_dofs; + + if (degree < 4) + return map; + + std::vector remainders((degree - 1) * (degree - 2) + * (degree - 3) / 6); + std::iota(remainders.begin(), remainders.end(), base); + + for (std::uint16_t r : vtk_tetrahedron_remainders(remainders)) + { + map[j++] = r; + } + + return map; } //----------------------------------------------------------------------------- std::vector vtk_wedge(int num_nodes) @@ -198,13 +417,7 @@ std::vector vtk_pyramid(int num_nodes) //----------------------------------------------------------------------------- std::vector vtk_quadrilateral(int num_nodes) { - // Check that num_nodes is a square integer (since quadrilaterals are - // tensor products of intervals, the number of nodes for each interval - // should be an integer) - assert((std::sqrt(num_nodes) - std::floor(std::sqrt(num_nodes))) == 0); - - // Number of nodes in each direction - const int n = sqrt(num_nodes); + const int n = cell_degree(mesh::CellType::quadrilateral, num_nodes); std::vector map(num_nodes); // Vertices @@ -215,7 +428,7 @@ std::vector vtk_quadrilateral(int num_nodes) int j = 4; - const int edge_nodes = n - 2; + const int edge_nodes = n - 1; // Edges for (int k = 0; k < edge_nodes; ++k) @@ -235,8 +448,7 @@ std::vector vtk_quadrilateral(int num_nodes) //----------------------------------------------------------------------------- std::vector vtk_hexahedron(int num_nodes) { - const int n = std::cbrt(num_nodes); - assert(n * n * n == num_nodes); + std::uint16_t n = cell_degree(mesh::CellType::hexahedron, num_nodes); std::vector map(num_nodes); @@ -253,7 +465,7 @@ std::vector vtk_hexahedron(int num_nodes) // Edges int j = 8; int base = 8; - const int edge_nodes = n - 2; + const int edge_nodes = n - 1; const std::vector edges = {0, 3, 5, 1, 8, 10, 11, 9, 2, 4, 7, 6}; for (int e : edges) { diff --git a/python/test/unit/mesh/test_higher_order_mesh.py b/python/test/unit/mesh/test_higher_order_mesh.py index efe2032ea74..a153596ac19 100644 --- a/python/test/unit/mesh/test_higher_order_mesh.py +++ b/python/test/unit/mesh/test_higher_order_mesh.py @@ -738,3 +738,152 @@ def test_quadrilateral_cell_order_3(dtype): "Q", "quadrilateral", 3, gdim=2, lagrange_variant=basix.LagrangeVariant.equispaced, shape=(2,), dtype=dtype)) check_cell_volume(points, cell, domain, 5 / 6, dtype=dtype) + + +@pytest.mark.parametrize('order', range(1, 11)) +def test_vtk_perm_tetrahedron(order): + size = (order + 1) * (order + 2) * (order + 3) // 6 + p = perm_vtk(CellType.tetrahedron, size) + + if order == 1: + q = [0, 1, 2, 3] + if order == 2: + q = [0, 1, 2, 3, 9, 6, 8, 7, 5, 4] + if order == 3: + q = [0, 1, 2, 3, 14, 15, 8, 9, 13, 12, 10, 11, 6, 7, 4, 5, 18, 16, 17, 19] + if order == 4: + q = [0, 1, 2, 3, 19, 20, 21, 10, 11, 12, 18, 17, 16, 13, 14, 15, 7, 8, 9, 4, 5, 6, + 28, 29, 30, 23, 24, 22, 25, 27, 26, 31, 33, 32, 34] + if order == 5: + q = [0, 1, 2, 3, 24, 25, 26, 27, 12, 13, 14, 15, 23, 22, 21, 20, 16, 17, 18, 19, 8, + 9, 10, 11, 4, 5, 6, 7, 40, 42, 45, 41, 44, 43, 30, 33, 28, 32, 31, 29, 34, 39, + 36, 37, 38, 35, 46, 51, 48, 49, 50, 47, 52, 53, 54, 55] + if order == 6: + q = [0, 1, 2, 3, 29, 30, 31, 32, 33, 14, 15, 16, 17, 18, 28, 27, 26, 25, 24, 19, 20, + 21, 22, 23, 9, 10, 11, 12, 13, 4, 5, 6, 7, 8, 54, 57, 63, 55, 56, 60, 62, 61, + 58, 59, 37, 43, 34, 40, 42, 41, 38, 35, 36, 39, 44, 53, 47, 48, 51, 52, 50, 46, + 45, 49, 64, 73, 67, 68, 71, 72, 70, 66, 65, 69, 74, 76, 79, 83, 75, 78, 77, 80, + 81, 82] + if order == 7: + q = [0, 1, 2, 3, 34, 35, 36, 37, 38, 39, 16, 17, 18, 19, 20, 21, 33, 32, 31, 30, 29, + 28, 22, 23, 24, 25, 26, 27, 10, 11, 12, 13, 14, 15, 4, 5, 6, 7, 8, 9, 70, 74, + 84, 71, 72, 73, 78, 81, 83, 82, 79, 75, 76, 77, 80, 44, 54, 40, 48, 51, 53, 52, + 49, 45, 41, 42, 43, 47, 50, 46, 55, 69, 59, 60, 64, 67, 68, 66, 63, 58, 57, 56, + 61, 65, 62, 85, 99, 89, 90, 94, 97, 98, 96, 93, 88, 87, 86, 91, 95, 92, 100, + 103, 109, 119, 101, 102, 106, 108, 107, 104, 110, 116, 112, 117, 115, 118, 111, + 114, 113, 105] + if order == 8: + q = [0, 1, 2, 3, 39, 40, 41, 42, 43, 44, 45, 18, 19, 20, 21, 22, 23, 24, 38, 37, 36, + 35, 34, 33, 32, 25, 26, 27, 28, 29, 30, 31, 11, 12, 13, 14, 15, 16, 17, 4, 5, 6, + 7, 8, 9, 10, 88, 93, 108, 89, 90, 91, 92, 98, 102, 105, 107, 106, 103, 99, 94, + 95, 97, 104, 96, 101, 100, 51, 66, 46, 56, 60, 63, 65, 64, 61, 57, 52, 47, 48, + 49, 50, 55, 62, 53, 59, 58, 54, 67, 87, 72, 73, 78, 82, 85, 86, 84, 81, 77, 71, + 70, 69, 68, 74, 83, 76, 79, 80, 75, 109, 129, 114, 115, 120, 124, 127, 128, 126, + 123, 119, 113, 112, 111, 110, 116, 125, 118, 121, 122, 117, 130, 134, 144, 164, + 131, 132, 133, 138, 141, 143, 142, 139, 135, 145, 155, 161, 148, 157, 162, 154, + 160, 163, 146, 147, 156, 153, 159, 151, 149, 158, 152, 136, 140, 137, 150] + if order == 9: + q = [0, 1, 2, 3, 44, 45, 46, 47, 48, 49, 50, 51, 20, 21, 22, 23, 24, 25, 26, 27, 43, + 42, 41, 40, 39, 38, 37, 36, 28, 29, 30, 31, 32, 33, 34, 35, 12, 13, 14, 15, 16, + 17, 18, 19, 4, 5, 6, 7, 8, 9, 10, 11, 108, 114, 135, 109, 110, 111, 112, 113, + 120, 125, 129, 132, 134, 133, 130, 126, 121, 115, 116, 119, 131, 117, 118, 124, + 128, 127, 122, 123, 58, 79, 52, 64, 69, 73, 76, 78, 77, 74, 70, 65, 59, 53, 54, + 55, 56, 57, 63, 75, 60, 68, 72, 71, 66, 61, 62, 67, 80, 107, 86, 87, 93, 98, + 102, 105, 106, 104, 101, 97, 92, 85, 84, 83, 82, 81, 88, 103, 91, 94, 99, 100, + 96, 90, 89, 95, 136, 163, 142, 143, 149, 154, 158, 161, 162, 160, 157, 153, 148, + 141, 140, 139, 138, 137, 144, 159, 147, 150, 155, 156, 152, 146, 145, 151, 164, + 169, 184, 219, 165, 166, 167, 168, 174, 178, 181, 183, 182, 179, 175, 170, 185, + 200, 210, 216, 189, 203, 212, 217, 199, 209, 215, 218, 186, 188, 211, 187, 202, + 201, 198, 214, 193, 208, 206, 196, 190, 213, 197, 204, 207, 194, 171, 180, 173, + 176, 177, 172, 191, 192, 195, 205] + if order == 10: + q = [0, 1, 2, 3, 49, 50, 51, 52, 53, 54, 55, 56, 57, 22, 23, 24, 25, 26, 27, 28, 29, + 30, 48, 47, 46, 45, 44, 43, 42, 41, 40, 31, 32, 33, 34, 35, 36, 37, 38, 39, 13, + 14, 15, 16, 17, 18, 19, 20, 21, 4, 5, 6, 7, 8, 9, 10, 11, 12, 130, 137, 165, + 131, 132, 133, 134, 135, 136, 144, 150, 155, 159, 162, 164, 163, 160, 156, 151, + 145, 138, 139, 143, 161, 140, 141, 142, 149, 154, 158, 157, 152, 146, 147, 148, + 153, 65, 93, 58, 72, 78, 83, 87, 90, 92, 91, 88, 84, 79, 73, 66, 59, 60, 61, 62, + 63, 64, 71, 89, 67, 77, 82, 86, 85, 80, 74, 68, 69, 70, 76, 81, 75, 94, 129, + 101, 102, 109, 115, 120, 124, 127, 128, 126, 123, 119, 114, 108, 100, 99, 98, + 97, 96, 95, 103, 125, 107, 110, 116, 121, 122, 118, 113, 106, 105, 104, 111, + 117, 112, 166, 201, 173, 174, 181, 187, 192, 196, 199, 200, 198, 195, 191, 186, + 180, 172, 171, 170, 169, 168, 167, 175, 197, 179, 182, 188, 193, 194, 190, 185, + 178, 177, 176, 183, 189, 184, 202, 208, 229, 285, 203, 204, 205, 206, 207, 214, + 219, 223, 226, 228, 227, 224, 220, 215, 209, 230, 251, 266, 276, 282, 235, 255, + 269, 278, 283, 250, 265, 275, 281, 284, 231, 234, 277, 232, 233, 254, 268, 267, + 252, 253, 249, 280, 240, 264, 274, 272, 259, 244, 247, 262, 236, 279, 248, 256, + 270, 273, 263, 245, 241, 260, 210, 225, 213, 216, 221, 222, 218, 212, 211, 217, + 237, 239, 246, 271, 238, 243, 242, 257, 258, 261] + + pt = [0 for i in p] + for i, j in enumerate(p): + pt[j] = i + print(" ".join([f"{i}" for i in pt])) + print(" ".join([f"{i}" for i in q])) + + for i, j in enumerate(p): + assert q[j] == i + + +@pytest.mark.parametrize('order', range(1, 7)) +def test_vtk_perm_hexahedron(order): + size = (order + 1) ** 3 + p = perm_vtk(CellType.hexahedron, size) + + if order == 1: + q = [0, 1, 3, 2, 4, 5, 7, 6] + if order == 2: + q = [0, 1, 3, 2, 4, 5, 7, 6, 8, 11, 13, 9, 16, 18, 19, 17, 10, 12, 15, 14, 22, 23, + 21, 24, 20, 25, 26] + if order == 3: + q = [0, 1, 3, 2, 4, 5, 7, 6, 8, 9, 14, 15, 18, 19, 10, 11, 24, 25, 28, 29, 30, 31, + 26, 27, 12, 13, 16, 17, 22, 23, 20, 21, 40, 41, 42, 43, 44, 45, 46, 47, 36, + 37, 38, 39, 48, 49, 50, 51, 32, 33, 34, 35, 52, 53, 54, 55, 56, 57, 58, 59, + 60, 61, 62, 63] + if order == 4: + q = [0, 1, 3, 2, 4, 5, 7, 6, 8, 9, 10, 17, 18, 19, 23, 24, 25, 11, 12, 13, 32, 33, + 34, 38, 39, 40, 41, 42, 43, 35, 36, 37, 14, 15, 16, 20, 21, 22, 29, 30, 31, + 26, 27, 28, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, + 78, 79, 53, 54, 55, 56, 57, 58, 59, 60, 61, 80, 81, 82, 83, 84, 85, 86, 87, + 88, 44, 45, 46, 47, 48, 49, 50, 51, 52, 89, 90, 91, 92, 93, 94, 95, 96, 97, + 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, + 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124] + if order == 5: + q = [0, 1, 3, 2, 4, 5, 7, 6, 8, 9, 10, 11, 20, 21, 22, 23, 28, 29, 30, 31, 12, 13, + 14, 15, 40, 41, 42, 43, 48, 49, 50, 51, 52, 53, 54, 55, 44, 45, 46, 47, 16, + 17, 18, 19, 24, 25, 26, 27, 36, 37, 38, 39, 32, 33, 34, 35, 88, 89, 90, 91, + 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, + 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 72, 73, 74, 75, 76, 77, + 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 120, 121, 122, 123, 124, 125, 126, + 127, 128, 129, 130, 131, 132, 133, 134, 135, 56, 57, 58, 59, 60, 61, 62, 63, + 64, 65, 66, 67, 68, 69, 70, 71, 136, 137, 138, 139, 140, 141, 142, 143, 144, + 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, + 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, + 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, + 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, + 209, 210, 211, 212, 213, 214, 215] + if order == 6: + q = [0, 1, 3, 2, 4, 5, 7, 6, 8, 9, 10, 11, 12, 23, 24, 25, 26, 27, 33, 34, 35, 36, + 37, 13, 14, 15, 16, 17, 48, 49, 50, 51, 52, 58, 59, 60, 61, 62, 63, 64, 65, 66, + 67, 53, 54, 55, 56, 57, 18, 19, 20, 21, 22, 28, 29, 30, 31, 32, 43, 44, 45, 46, + 47, 38, 39, 40, 41, 42, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, + 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, + 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, + 161, 162, 163, 164, 165, 166, 167, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, + 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 168, + 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, + 185, 186, 187, 188, 189, 190, 191, 192, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, + 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 193, 194, 195, 196, + 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, + 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, + 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, + 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, + 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, + 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, + 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, + 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, + 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, + 341, 342] + + for i, j in enumerate(p): + assert q[j] == i