From fc1d2b1c135231da9621a8d6c1d80692e0c59fe6 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Mon, 5 Oct 2026 23:40:27 +0200 Subject: [PATCH 1/2] Pipe sector with nodes on the rotation axis for the periodic3d_axis test 30 degree sector of a pipe with one cell in the circumferential direction, written by wedge.py (python3 wedge.py wedge.su2 30 1 1 8 10). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UkNcoCEH8nFNrHWzhJCar --- navierstokes/periodic3D_axis/wedge.py | 55 +++ navierstokes/periodic3D_axis/wedge.su2 | 467 +++++++++++++++++++++++++ 2 files changed, 522 insertions(+) create mode 100644 navierstokes/periodic3D_axis/wedge.py create mode 100644 navierstokes/periodic3D_axis/wedge.su2 diff --git a/navierstokes/periodic3D_axis/wedge.py b/navierstokes/periodic3D_axis/wedge.py new file mode 100644 index 00000000..e5cae344 --- /dev/null +++ b/navierstokes/periodic3D_axis/wedge.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python3 +"""wedge.py [r0] +3D pipe (radius 1, length 2, axis = z) made of nsect sectors of degrees with ncirc cells each in the +circumferential direction. If nsect * angle = 360 the mesh is the full pipe (no periodic markers), otherwise the +markers per1 (theta = 0) and per2 (theta = nsect * angle) are written. +r0 = 0 (default): the mesh contains the axis (prisms at the axis, the axis nodes lie on both periodic markers); +r0 > 0: annular mesh without axis nodes (marker hub at r = r0). +Markers: per1, per2, inlet (z = 0), outlet (z = 2), wall (r = 1) [, hub]. +Also writes .map with one line per node: node, j (circumferential index), i (radial index), k (axial +index); axis nodes have j = 0, i = 0. Nodes with the same (j, i, k) are the same node in a sector and in the full pipe.""" +import sys, math +out, deg, ncirc, nsect, nr, nz = sys.argv[1], float(sys.argv[2]), int(sys.argv[3]), int(sys.argv[4]), int(sys.argv[5]), int(sys.argv[6]) +r0 = float(sys.argv[7]) if len(sys.argv) > 7 else 0.0 +R, Lz = 1.0, 2.0 +M = ncirc * nsect # cells in the circumferential direction +full = abs(nsect * deg - 360.0) < 1e-9 +axis = r0 == 0.0 +X, ijk = [], [] +def add(x, y, z, j, i, k): X.append((x, y, z)); ijk.append((j, i, k)); return len(X) - 1 +A = [add(0.0, 0.0, Lz*k/nz, 0, 0, k) for k in range(nz+1)] if axis else None +i0 = 1 if axis else 0 +nj = M if full else M + 1 # node lines in the circumferential direction +P = [[[None]*(nz+1) for _ in range(nr+1)] for _ in range(nj)] +for j in range(nj): + th = math.radians(deg) * j / ncirc + for i in range(i0, nr+1): + r = r0 + (R-r0)*i/nr + for k in range(nz+1): P[j][i][k] = add(r*math.cos(th), r*math.sin(th), Lz*k/nz, j, i, k) +E, Mk = [], {"per1": [], "per2": [], "inlet": [], "outlet": [], "wall": []} +if full: del Mk["per1"], Mk["per2"] +if not axis: Mk["hub"] = [] +for k in range(nz): + for j in range(M): + a, b = P[j], P[(j+1) % nj] + if axis: E.append((13, A[k], a[1][k], b[1][k], A[k+1], a[1][k+1], b[1][k+1])) + else: Mk["hub"].append((9, a[0][k], b[0][k], b[0][k+1], a[0][k+1])) + for i in range(i0, nr): + E.append((12, a[i][k], a[i+1][k], b[i+1][k], b[i][k], a[i][k+1], a[i+1][k+1], b[i+1][k+1], b[i][k+1])) + Mk["wall"].append((9, a[nr][k], b[nr][k], b[nr][k+1], a[nr][k+1])) + if not full: + for j, m in ((0, "per1"), (M, "per2")): + if axis: Mk[m].append((9, A[k], P[j][1][k], P[j][1][k+1], A[k+1])) + for i in range(i0, nr): Mk[m].append((9, P[j][i][k], P[j][i+1][k], P[j][i+1][k+1], P[j][i][k+1])) +for k, m in ((0, "inlet"), (nz, "outlet")): + for j in range(M): + a, b = P[j], P[(j+1) % nj] + if axis: Mk[m].append((5, A[k], a[1][k], b[1][k])) + for i in range(i0, nr): Mk[m].append((9, a[i][k], a[i+1][k], b[i+1][k], b[i][k])) +L = ["NDIME= 3", "NELEM= %d" % len(E)] + [" ".join(map(str, e)) + " %d" % n for n, e in enumerate(E)] +L += ["NPOIN= %d" % len(X)] + ["%.17g %.17g %.17g %d" % (x[0], x[1], x[2], n) for n, x in enumerate(X)] +L += ["NMARK= %d" % len(Mk)] +for m, f in Mk.items(): L += ["MARKER_TAG= " + m, "MARKER_ELEMS= %d" % len(f)] + [" ".join(map(str, e)) for e in f] +open(out, "w").write("\n".join(L) + "\n") +open(out + ".map", "w").write("\n".join("%d %d %d %d" % ((n,) + t) for n, t in enumerate(ijk)) + "\n") +print(out, "points", len(X), "elements", len(E), "axis nodes", len(A) if axis else 0, "full" if full else "sector") diff --git a/navierstokes/periodic3D_axis/wedge.su2 b/navierstokes/periodic3D_axis/wedge.su2 new file mode 100644 index 00000000..d474b746 --- /dev/null +++ b/navierstokes/periodic3D_axis/wedge.su2 @@ -0,0 +1,467 @@ +NDIME= 3 +NELEM= 80 +13 0 11 99 1 12 100 0 +12 11 22 110 99 12 23 111 100 1 +12 22 33 121 110 23 34 122 111 2 +12 33 44 132 121 34 45 133 122 3 +12 44 55 143 132 45 56 144 133 4 +12 55 66 154 143 56 67 155 144 5 +12 66 77 165 154 67 78 166 155 6 +12 77 88 176 165 78 89 177 166 7 +13 1 12 100 2 13 101 8 +12 12 23 111 100 13 24 112 101 9 +12 23 34 122 111 24 35 123 112 10 +12 34 45 133 122 35 46 134 123 11 +12 45 56 144 133 46 57 145 134 12 +12 56 67 155 144 57 68 156 145 13 +12 67 78 166 155 68 79 167 156 14 +12 78 89 177 166 79 90 178 167 15 +13 2 13 101 3 14 102 16 +12 13 24 112 101 14 25 113 102 17 +12 24 35 123 112 25 36 124 113 18 +12 35 46 134 123 36 47 135 124 19 +12 46 57 145 134 47 58 146 135 20 +12 57 68 156 145 58 69 157 146 21 +12 68 79 167 156 69 80 168 157 22 +12 79 90 178 167 80 91 179 168 23 +13 3 14 102 4 15 103 24 +12 14 25 113 102 15 26 114 103 25 +12 25 36 124 113 26 37 125 114 26 +12 36 47 135 124 37 48 136 125 27 +12 47 58 146 135 48 59 147 136 28 +12 58 69 157 146 59 70 158 147 29 +12 69 80 168 157 70 81 169 158 30 +12 80 91 179 168 81 92 180 169 31 +13 4 15 103 5 16 104 32 +12 15 26 114 103 16 27 115 104 33 +12 26 37 125 114 27 38 126 115 34 +12 37 48 136 125 38 49 137 126 35 +12 48 59 147 136 49 60 148 137 36 +12 59 70 158 147 60 71 159 148 37 +12 70 81 169 158 71 82 170 159 38 +12 81 92 180 169 82 93 181 170 39 +13 5 16 104 6 17 105 40 +12 16 27 115 104 17 28 116 105 41 +12 27 38 126 115 28 39 127 116 42 +12 38 49 137 126 39 50 138 127 43 +12 49 60 148 137 50 61 149 138 44 +12 60 71 159 148 61 72 160 149 45 +12 71 82 170 159 72 83 171 160 46 +12 82 93 181 170 83 94 182 171 47 +13 6 17 105 7 18 106 48 +12 17 28 116 105 18 29 117 106 49 +12 28 39 127 116 29 40 128 117 50 +12 39 50 138 127 40 51 139 128 51 +12 50 61 149 138 51 62 150 139 52 +12 61 72 160 149 62 73 161 150 53 +12 72 83 171 160 73 84 172 161 54 +12 83 94 182 171 84 95 183 172 55 +13 7 18 106 8 19 107 56 +12 18 29 117 106 19 30 118 107 57 +12 29 40 128 117 30 41 129 118 58 +12 40 51 139 128 41 52 140 129 59 +12 51 62 150 139 52 63 151 140 60 +12 62 73 161 150 63 74 162 151 61 +12 73 84 172 161 74 85 173 162 62 +12 84 95 183 172 85 96 184 173 63 +13 8 19 107 9 20 108 64 +12 19 30 118 107 20 31 119 108 65 +12 30 41 129 118 31 42 130 119 66 +12 41 52 140 129 42 53 141 130 67 +12 52 63 151 140 53 64 152 141 68 +12 63 74 162 151 64 75 163 152 69 +12 74 85 173 162 75 86 174 163 70 +12 85 96 184 173 86 97 185 174 71 +13 9 20 108 10 21 109 72 +12 20 31 119 108 21 32 120 109 73 +12 31 42 130 119 32 43 131 120 74 +12 42 53 141 130 43 54 142 131 75 +12 53 64 152 141 54 65 153 142 76 +12 64 75 163 152 65 76 164 153 77 +12 75 86 174 163 76 87 175 164 78 +12 86 97 185 174 87 98 186 175 79 +NPOIN= 187 +0 0 0 0 +0 0 0.20000000000000001 1 +0 0 0.40000000000000002 2 +0 0 0.59999999999999998 3 +0 0 0.80000000000000004 4 +0 0 1 5 +0 0 1.2 6 +0 0 1.3999999999999999 7 +0 0 1.6000000000000001 8 +0 0 1.8 9 +0 0 2 10 +0.125 0 0 11 +0.125 0 0.20000000000000001 12 +0.125 0 0.40000000000000002 13 +0.125 0 0.59999999999999998 14 +0.125 0 0.80000000000000004 15 +0.125 0 1 16 +0.125 0 1.2 17 +0.125 0 1.3999999999999999 18 +0.125 0 1.6000000000000001 19 +0.125 0 1.8 20 +0.125 0 2 21 +0.25 0 0 22 +0.25 0 0.20000000000000001 23 +0.25 0 0.40000000000000002 24 +0.25 0 0.59999999999999998 25 +0.25 0 0.80000000000000004 26 +0.25 0 1 27 +0.25 0 1.2 28 +0.25 0 1.3999999999999999 29 +0.25 0 1.6000000000000001 30 +0.25 0 1.8 31 +0.25 0 2 32 +0.375 0 0 33 +0.375 0 0.20000000000000001 34 +0.375 0 0.40000000000000002 35 +0.375 0 0.59999999999999998 36 +0.375 0 0.80000000000000004 37 +0.375 0 1 38 +0.375 0 1.2 39 +0.375 0 1.3999999999999999 40 +0.375 0 1.6000000000000001 41 +0.375 0 1.8 42 +0.375 0 2 43 +0.5 0 0 44 +0.5 0 0.20000000000000001 45 +0.5 0 0.40000000000000002 46 +0.5 0 0.59999999999999998 47 +0.5 0 0.80000000000000004 48 +0.5 0 1 49 +0.5 0 1.2 50 +0.5 0 1.3999999999999999 51 +0.5 0 1.6000000000000001 52 +0.5 0 1.8 53 +0.5 0 2 54 +0.625 0 0 55 +0.625 0 0.20000000000000001 56 +0.625 0 0.40000000000000002 57 +0.625 0 0.59999999999999998 58 +0.625 0 0.80000000000000004 59 +0.625 0 1 60 +0.625 0 1.2 61 +0.625 0 1.3999999999999999 62 +0.625 0 1.6000000000000001 63 +0.625 0 1.8 64 +0.625 0 2 65 +0.75 0 0 66 +0.75 0 0.20000000000000001 67 +0.75 0 0.40000000000000002 68 +0.75 0 0.59999999999999998 69 +0.75 0 0.80000000000000004 70 +0.75 0 1 71 +0.75 0 1.2 72 +0.75 0 1.3999999999999999 73 +0.75 0 1.6000000000000001 74 +0.75 0 1.8 75 +0.75 0 2 76 +0.875 0 0 77 +0.875 0 0.20000000000000001 78 +0.875 0 0.40000000000000002 79 +0.875 0 0.59999999999999998 80 +0.875 0 0.80000000000000004 81 +0.875 0 1 82 +0.875 0 1.2 83 +0.875 0 1.3999999999999999 84 +0.875 0 1.6000000000000001 85 +0.875 0 1.8 86 +0.875 0 2 87 +1 0 0 88 +1 0 0.20000000000000001 89 +1 0 0.40000000000000002 90 +1 0 0.59999999999999998 91 +1 0 0.80000000000000004 92 +1 0 1 93 +1 0 1.2 94 +1 0 1.3999999999999999 95 +1 0 1.6000000000000001 96 +1 0 1.8 97 +1 0 2 98 +0.10825317547305484 0.062499999999999993 0 99 +0.10825317547305484 0.062499999999999993 0.20000000000000001 100 +0.10825317547305484 0.062499999999999993 0.40000000000000002 101 +0.10825317547305484 0.062499999999999993 0.59999999999999998 102 +0.10825317547305484 0.062499999999999993 0.80000000000000004 103 +0.10825317547305484 0.062499999999999993 1 104 +0.10825317547305484 0.062499999999999993 1.2 105 +0.10825317547305484 0.062499999999999993 1.3999999999999999 106 +0.10825317547305484 0.062499999999999993 1.6000000000000001 107 +0.10825317547305484 0.062499999999999993 1.8 108 +0.10825317547305484 0.062499999999999993 2 109 +0.21650635094610968 0.12499999999999999 0 110 +0.21650635094610968 0.12499999999999999 0.20000000000000001 111 +0.21650635094610968 0.12499999999999999 0.40000000000000002 112 +0.21650635094610968 0.12499999999999999 0.59999999999999998 113 +0.21650635094610968 0.12499999999999999 0.80000000000000004 114 +0.21650635094610968 0.12499999999999999 1 115 +0.21650635094610968 0.12499999999999999 1.2 116 +0.21650635094610968 0.12499999999999999 1.3999999999999999 117 +0.21650635094610968 0.12499999999999999 1.6000000000000001 118 +0.21650635094610968 0.12499999999999999 1.8 119 +0.21650635094610968 0.12499999999999999 2 120 +0.3247595264191645 0.18749999999999997 0 121 +0.3247595264191645 0.18749999999999997 0.20000000000000001 122 +0.3247595264191645 0.18749999999999997 0.40000000000000002 123 +0.3247595264191645 0.18749999999999997 0.59999999999999998 124 +0.3247595264191645 0.18749999999999997 0.80000000000000004 125 +0.3247595264191645 0.18749999999999997 1 126 +0.3247595264191645 0.18749999999999997 1.2 127 +0.3247595264191645 0.18749999999999997 1.3999999999999999 128 +0.3247595264191645 0.18749999999999997 1.6000000000000001 129 +0.3247595264191645 0.18749999999999997 1.8 130 +0.3247595264191645 0.18749999999999997 2 131 +0.43301270189221935 0.24999999999999997 0 132 +0.43301270189221935 0.24999999999999997 0.20000000000000001 133 +0.43301270189221935 0.24999999999999997 0.40000000000000002 134 +0.43301270189221935 0.24999999999999997 0.59999999999999998 135 +0.43301270189221935 0.24999999999999997 0.80000000000000004 136 +0.43301270189221935 0.24999999999999997 1 137 +0.43301270189221935 0.24999999999999997 1.2 138 +0.43301270189221935 0.24999999999999997 1.3999999999999999 139 +0.43301270189221935 0.24999999999999997 1.6000000000000001 140 +0.43301270189221935 0.24999999999999997 1.8 141 +0.43301270189221935 0.24999999999999997 2 142 +0.54126587736527421 0.31249999999999994 0 143 +0.54126587736527421 0.31249999999999994 0.20000000000000001 144 +0.54126587736527421 0.31249999999999994 0.40000000000000002 145 +0.54126587736527421 0.31249999999999994 0.59999999999999998 146 +0.54126587736527421 0.31249999999999994 0.80000000000000004 147 +0.54126587736527421 0.31249999999999994 1 148 +0.54126587736527421 0.31249999999999994 1.2 149 +0.54126587736527421 0.31249999999999994 1.3999999999999999 150 +0.54126587736527421 0.31249999999999994 1.6000000000000001 151 +0.54126587736527421 0.31249999999999994 1.8 152 +0.54126587736527421 0.31249999999999994 2 153 +0.649519052838329 0.37499999999999994 0 154 +0.649519052838329 0.37499999999999994 0.20000000000000001 155 +0.649519052838329 0.37499999999999994 0.40000000000000002 156 +0.649519052838329 0.37499999999999994 0.59999999999999998 157 +0.649519052838329 0.37499999999999994 0.80000000000000004 158 +0.649519052838329 0.37499999999999994 1 159 +0.649519052838329 0.37499999999999994 1.2 160 +0.649519052838329 0.37499999999999994 1.3999999999999999 161 +0.649519052838329 0.37499999999999994 1.6000000000000001 162 +0.649519052838329 0.37499999999999994 1.8 163 +0.649519052838329 0.37499999999999994 2 164 +0.75777222831138391 0.43749999999999994 0 165 +0.75777222831138391 0.43749999999999994 0.20000000000000001 166 +0.75777222831138391 0.43749999999999994 0.40000000000000002 167 +0.75777222831138391 0.43749999999999994 0.59999999999999998 168 +0.75777222831138391 0.43749999999999994 0.80000000000000004 169 +0.75777222831138391 0.43749999999999994 1 170 +0.75777222831138391 0.43749999999999994 1.2 171 +0.75777222831138391 0.43749999999999994 1.3999999999999999 172 +0.75777222831138391 0.43749999999999994 1.6000000000000001 173 +0.75777222831138391 0.43749999999999994 1.8 174 +0.75777222831138391 0.43749999999999994 2 175 +0.86602540378443871 0.49999999999999994 0 176 +0.86602540378443871 0.49999999999999994 0.20000000000000001 177 +0.86602540378443871 0.49999999999999994 0.40000000000000002 178 +0.86602540378443871 0.49999999999999994 0.59999999999999998 179 +0.86602540378443871 0.49999999999999994 0.80000000000000004 180 +0.86602540378443871 0.49999999999999994 1 181 +0.86602540378443871 0.49999999999999994 1.2 182 +0.86602540378443871 0.49999999999999994 1.3999999999999999 183 +0.86602540378443871 0.49999999999999994 1.6000000000000001 184 +0.86602540378443871 0.49999999999999994 1.8 185 +0.86602540378443871 0.49999999999999994 2 186 +NMARK= 5 +MARKER_TAG= per1 +MARKER_ELEMS= 80 +9 0 11 12 1 +9 11 22 23 12 +9 22 33 34 23 +9 33 44 45 34 +9 44 55 56 45 +9 55 66 67 56 +9 66 77 78 67 +9 77 88 89 78 +9 1 12 13 2 +9 12 23 24 13 +9 23 34 35 24 +9 34 45 46 35 +9 45 56 57 46 +9 56 67 68 57 +9 67 78 79 68 +9 78 89 90 79 +9 2 13 14 3 +9 13 24 25 14 +9 24 35 36 25 +9 35 46 47 36 +9 46 57 58 47 +9 57 68 69 58 +9 68 79 80 69 +9 79 90 91 80 +9 3 14 15 4 +9 14 25 26 15 +9 25 36 37 26 +9 36 47 48 37 +9 47 58 59 48 +9 58 69 70 59 +9 69 80 81 70 +9 80 91 92 81 +9 4 15 16 5 +9 15 26 27 16 +9 26 37 38 27 +9 37 48 49 38 +9 48 59 60 49 +9 59 70 71 60 +9 70 81 82 71 +9 81 92 93 82 +9 5 16 17 6 +9 16 27 28 17 +9 27 38 39 28 +9 38 49 50 39 +9 49 60 61 50 +9 60 71 72 61 +9 71 82 83 72 +9 82 93 94 83 +9 6 17 18 7 +9 17 28 29 18 +9 28 39 40 29 +9 39 50 51 40 +9 50 61 62 51 +9 61 72 73 62 +9 72 83 84 73 +9 83 94 95 84 +9 7 18 19 8 +9 18 29 30 19 +9 29 40 41 30 +9 40 51 52 41 +9 51 62 63 52 +9 62 73 74 63 +9 73 84 85 74 +9 84 95 96 85 +9 8 19 20 9 +9 19 30 31 20 +9 30 41 42 31 +9 41 52 53 42 +9 52 63 64 53 +9 63 74 75 64 +9 74 85 86 75 +9 85 96 97 86 +9 9 20 21 10 +9 20 31 32 21 +9 31 42 43 32 +9 42 53 54 43 +9 53 64 65 54 +9 64 75 76 65 +9 75 86 87 76 +9 86 97 98 87 +MARKER_TAG= per2 +MARKER_ELEMS= 80 +9 0 99 100 1 +9 99 110 111 100 +9 110 121 122 111 +9 121 132 133 122 +9 132 143 144 133 +9 143 154 155 144 +9 154 165 166 155 +9 165 176 177 166 +9 1 100 101 2 +9 100 111 112 101 +9 111 122 123 112 +9 122 133 134 123 +9 133 144 145 134 +9 144 155 156 145 +9 155 166 167 156 +9 166 177 178 167 +9 2 101 102 3 +9 101 112 113 102 +9 112 123 124 113 +9 123 134 135 124 +9 134 145 146 135 +9 145 156 157 146 +9 156 167 168 157 +9 167 178 179 168 +9 3 102 103 4 +9 102 113 114 103 +9 113 124 125 114 +9 124 135 136 125 +9 135 146 147 136 +9 146 157 158 147 +9 157 168 169 158 +9 168 179 180 169 +9 4 103 104 5 +9 103 114 115 104 +9 114 125 126 115 +9 125 136 137 126 +9 136 147 148 137 +9 147 158 159 148 +9 158 169 170 159 +9 169 180 181 170 +9 5 104 105 6 +9 104 115 116 105 +9 115 126 127 116 +9 126 137 138 127 +9 137 148 149 138 +9 148 159 160 149 +9 159 170 171 160 +9 170 181 182 171 +9 6 105 106 7 +9 105 116 117 106 +9 116 127 128 117 +9 127 138 139 128 +9 138 149 150 139 +9 149 160 161 150 +9 160 171 172 161 +9 171 182 183 172 +9 7 106 107 8 +9 106 117 118 107 +9 117 128 129 118 +9 128 139 140 129 +9 139 150 151 140 +9 150 161 162 151 +9 161 172 173 162 +9 172 183 184 173 +9 8 107 108 9 +9 107 118 119 108 +9 118 129 130 119 +9 129 140 141 130 +9 140 151 152 141 +9 151 162 163 152 +9 162 173 174 163 +9 173 184 185 174 +9 9 108 109 10 +9 108 119 120 109 +9 119 130 131 120 +9 130 141 142 131 +9 141 152 153 142 +9 152 163 164 153 +9 163 174 175 164 +9 174 185 186 175 +MARKER_TAG= inlet +MARKER_ELEMS= 8 +5 0 11 99 +9 11 22 110 99 +9 22 33 121 110 +9 33 44 132 121 +9 44 55 143 132 +9 55 66 154 143 +9 66 77 165 154 +9 77 88 176 165 +MARKER_TAG= outlet +MARKER_ELEMS= 8 +5 10 21 109 +9 21 32 120 109 +9 32 43 131 120 +9 43 54 142 131 +9 54 65 153 142 +9 65 76 164 153 +9 76 87 175 164 +9 87 98 186 175 +MARKER_TAG= wall +MARKER_ELEMS= 10 +9 88 176 177 89 +9 89 177 178 90 +9 90 178 179 91 +9 91 179 180 92 +9 92 180 181 93 +9 93 181 182 94 +9 94 182 183 95 +9 95 183 184 96 +9 96 184 185 97 +9 97 185 186 98 From f0b9b4a4862e3d3178b749af5739d7c8321e914e Mon Sep 17 00:00:00 2001 From: rois1995 Date: Sat, 10 Oct 2026 15:13:13 +0200 Subject: [PATCH 2/2] Clarify the periodic pipe mesh generator --- navierstokes/periodic3D_axis/wedge.py | 145 ++++++++++++++++++-------- 1 file changed, 100 insertions(+), 45 deletions(-) diff --git a/navierstokes/periodic3D_axis/wedge.py b/navierstokes/periodic3D_axis/wedge.py index e5cae344..baa30fae 100644 --- a/navierstokes/periodic3D_axis/wedge.py +++ b/navierstokes/periodic3D_axis/wedge.py @@ -8,48 +8,103 @@ Markers: per1, per2, inlet (z = 0), outlet (z = 2), wall (r = 1) [, hub]. Also writes .map with one line per node: node, j (circumferential index), i (radial index), k (axial index); axis nodes have j = 0, i = 0. Nodes with the same (j, i, k) are the same node in a sector and in the full pipe.""" -import sys, math -out, deg, ncirc, nsect, nr, nz = sys.argv[1], float(sys.argv[2]), int(sys.argv[3]), int(sys.argv[4]), int(sys.argv[5]), int(sys.argv[6]) -r0 = float(sys.argv[7]) if len(sys.argv) > 7 else 0.0 -R, Lz = 1.0, 2.0 -M = ncirc * nsect # cells in the circumferential direction -full = abs(nsect * deg - 360.0) < 1e-9 -axis = r0 == 0.0 -X, ijk = [], [] -def add(x, y, z, j, i, k): X.append((x, y, z)); ijk.append((j, i, k)); return len(X) - 1 -A = [add(0.0, 0.0, Lz*k/nz, 0, 0, k) for k in range(nz+1)] if axis else None -i0 = 1 if axis else 0 -nj = M if full else M + 1 # node lines in the circumferential direction -P = [[[None]*(nz+1) for _ in range(nr+1)] for _ in range(nj)] -for j in range(nj): - th = math.radians(deg) * j / ncirc - for i in range(i0, nr+1): - r = r0 + (R-r0)*i/nr - for k in range(nz+1): P[j][i][k] = add(r*math.cos(th), r*math.sin(th), Lz*k/nz, j, i, k) -E, Mk = [], {"per1": [], "per2": [], "inlet": [], "outlet": [], "wall": []} -if full: del Mk["per1"], Mk["per2"] -if not axis: Mk["hub"] = [] -for k in range(nz): - for j in range(M): - a, b = P[j], P[(j+1) % nj] - if axis: E.append((13, A[k], a[1][k], b[1][k], A[k+1], a[1][k+1], b[1][k+1])) - else: Mk["hub"].append((9, a[0][k], b[0][k], b[0][k+1], a[0][k+1])) - for i in range(i0, nr): - E.append((12, a[i][k], a[i+1][k], b[i+1][k], b[i][k], a[i][k+1], a[i+1][k+1], b[i+1][k+1], b[i][k+1])) - Mk["wall"].append((9, a[nr][k], b[nr][k], b[nr][k+1], a[nr][k+1])) - if not full: - for j, m in ((0, "per1"), (M, "per2")): - if axis: Mk[m].append((9, A[k], P[j][1][k], P[j][1][k+1], A[k+1])) - for i in range(i0, nr): Mk[m].append((9, P[j][i][k], P[j][i+1][k], P[j][i+1][k+1], P[j][i][k+1])) -for k, m in ((0, "inlet"), (nz, "outlet")): - for j in range(M): - a, b = P[j], P[(j+1) % nj] - if axis: Mk[m].append((5, A[k], a[1][k], b[1][k])) - for i in range(i0, nr): Mk[m].append((9, a[i][k], a[i+1][k], b[i+1][k], b[i][k])) -L = ["NDIME= 3", "NELEM= %d" % len(E)] + [" ".join(map(str, e)) + " %d" % n for n, e in enumerate(E)] -L += ["NPOIN= %d" % len(X)] + ["%.17g %.17g %.17g %d" % (x[0], x[1], x[2], n) for n, x in enumerate(X)] -L += ["NMARK= %d" % len(Mk)] -for m, f in Mk.items(): L += ["MARKER_TAG= " + m, "MARKER_ELEMS= %d" % len(f)] + [" ".join(map(str, e)) for e in f] -open(out, "w").write("\n".join(L) + "\n") -open(out + ".map", "w").write("\n".join("%d %d %d %d" % ((n,) + t) for n, t in enumerate(ijk)) + "\n") -print(out, "points", len(X), "elements", len(E), "axis nodes", len(A) if axis else 0, "full" if full else "sector") + +import math +import sys + +output = sys.argv[1] +sector_angle = float(sys.argv[2]) +circumferential_cells_per_sector = int(sys.argv[3]) +sectors = int(sys.argv[4]) +radial_cells = int(sys.argv[5]) +axial_cells = int(sys.argv[6]) +inner_radius = float(sys.argv[7]) if len(sys.argv) > 7 else 0.0 +radius, length = 1.0, 2.0 +circumferential_cells = circumferential_cells_per_sector * sectors +full_pipe = abs(sectors * sector_angle - 360.0) < 1e-9 +includes_axis = inner_radius == 0.0 +points, node_indices = [], [] + + +def add_point(x, y, z, j, i, k): + points.append((x, y, z)) + node_indices.append((j, i, k)) + return len(points) - 1 + + +# All circumferential node lines share a single axial line at r = 0. +axis_nodes = [] +if includes_axis: + axis_nodes = [add_point(0.0, 0.0, length * k / axial_cells, 0, 0, k) + for k in range(axial_cells + 1)] +first_radial_index = 1 if includes_axis else 0 +circumferential_node_lines = circumferential_cells if full_pipe else circumferential_cells + 1 +node_ids = [[[None] * (axial_cells + 1) for _ in range(radial_cells + 1)] + for _ in range(circumferential_node_lines)] +for j in range(circumferential_node_lines): + theta = math.radians(sector_angle) * j / circumferential_cells_per_sector + for i in range(first_radial_index, radial_cells + 1): + r = inner_radius + (radius - inner_radius) * i / radial_cells + for k in range(axial_cells + 1): + node_ids[j][i][k] = add_point(r * math.cos(theta), r * math.sin(theta), + length * k / axial_cells, j, i, k) + +elements = [] +markers = {"per1": [], "per2": [], "inlet": [], "outlet": [], "wall": []} +if full_pipe: + del markers["per1"], markers["per2"] +if not includes_axis: + markers["hub"] = [] + +# SU2 element types: triangle = 5, quadrilateral = 9, hexahedron = 12, prism = 13. +for k in range(axial_cells): + for j in range(circumferential_cells): + # Wrap the last node line onto the first one only for the full pipe. + left = node_ids[j] + right = node_ids[(j + 1) % circumferential_node_lines] + if includes_axis: + elements.append((13, axis_nodes[k], left[1][k], right[1][k], + axis_nodes[k + 1], left[1][k + 1], right[1][k + 1])) + else: + markers["hub"].append((9, left[0][k], right[0][k], right[0][k + 1], left[0][k + 1])) + for i in range(first_radial_index, radial_cells): + elements.append((12, left[i][k], left[i + 1][k], right[i + 1][k], right[i][k], + left[i][k + 1], left[i + 1][k + 1], right[i + 1][k + 1], right[i][k + 1])) + markers["wall"].append((9, left[radial_cells][k], right[radial_cells][k], + right[radial_cells][k + 1], left[radial_cells][k + 1])) + if not full_pipe: + for j, marker in ((0, "per1"), (circumferential_cells, "per2")): + nodes = node_ids[j] + if includes_axis: + markers[marker].append((9, axis_nodes[k], nodes[1][k], nodes[1][k + 1], axis_nodes[k + 1])) + for i in range(first_radial_index, radial_cells): + markers[marker].append((9, nodes[i][k], nodes[i + 1][k], nodes[i + 1][k + 1], nodes[i][k + 1])) + +for k, marker in ((0, "inlet"), (axial_cells, "outlet")): + for j in range(circumferential_cells): + left = node_ids[j] + right = node_ids[(j + 1) % circumferential_node_lines] + if includes_axis: + markers[marker].append((5, axis_nodes[k], left[1][k], right[1][k])) + for i in range(first_radial_index, radial_cells): + markers[marker].append((9, left[i][k], left[i + 1][k], right[i + 1][k], right[i][k])) + +with open(output, "w") as mesh_file: + mesh_file.write("NDIME= 3\nNELEM= %d\n" % len(elements)) + for element_id, element in enumerate(elements): + mesh_file.write(" ".join(map(str, element)) + " %d\n" % element_id) + mesh_file.write("NPOIN= %d\n" % len(points)) + for node_id, (x, y, z) in enumerate(points): + mesh_file.write("%.17g %.17g %.17g %d\n" % (x, y, z, node_id)) + mesh_file.write("NMARK= %d\n" % len(markers)) + for marker, faces in markers.items(): + mesh_file.write("MARKER_TAG= %s\nMARKER_ELEMS= %d\n" % (marker, len(faces))) + for face in faces: + mesh_file.write(" ".join(map(str, face)) + "\n") + +with open(output + ".map", "w") as map_file: + for node_id, indices in enumerate(node_indices): + map_file.write("%d %d %d %d\n" % ((node_id,) + indices)) + +print(output, "points", len(points), "elements", len(elements), "axis nodes", len(axis_nodes), + "full" if full_pipe else "sector")