|
| 1 | +# 1. reading a single VTU file |
| 2 | + |
| 3 | +A single VTU file can be accessed via: |
| 4 | + |
| 5 | + |
| 6 | +```python |
| 7 | +import vtuIO |
| 8 | +``` |
| 9 | + |
| 10 | + |
| 11 | +```python |
| 12 | +vtufile = vtuIO.VTUIO("examples/square_1e2_pcs_0_ts_1_t_1.000000.vtu", dim=2) |
| 13 | +``` |
| 14 | + |
| 15 | +The `dim` argument is needed for correct interpolation. By defualt `dim=3` is assumed. |
| 16 | +Basic VTU properties, like fieldnames, points and corresponding fielddata as provided by the unstructured grid VTK class can be simply accessed as follows: |
| 17 | + |
| 18 | + |
| 19 | +```python |
| 20 | +vtufile.getFieldnames() |
| 21 | +``` |
| 22 | + |
| 23 | + |
| 24 | + |
| 25 | + |
| 26 | + ['D1_left_bottom_N1_right', 'Linear_1_to_minus1', 'pressure', 'v'] |
| 27 | + |
| 28 | + |
| 29 | + |
| 30 | +points (in this example the first 3 points) can be simmply accessed with |
| 31 | + |
| 32 | + |
| 33 | +```python |
| 34 | +vtufile.points[0:3] |
| 35 | +``` |
| 36 | + |
| 37 | + |
| 38 | + |
| 39 | + |
| 40 | + array([[0. , 0. ], |
| 41 | + [0.1, 0. ], |
| 42 | + [0.2, 0. ]]) |
| 43 | + |
| 44 | + |
| 45 | + |
| 46 | + |
| 47 | +```python |
| 48 | +vtufile.getField("v")[0:3] |
| 49 | +``` |
| 50 | + |
| 51 | + |
| 52 | + |
| 53 | + |
| 54 | + array([[ 2.00000000e+00, 0.00000000e+00], |
| 55 | + [ 2.00000000e+00, 1.62547932e-16], |
| 56 | + [ 2.00000000e+00, -9.91229679e-16]]) |
| 57 | + |
| 58 | + |
| 59 | + |
| 60 | +Aside basic VTU properties, the field data at any given point, e.g, |
| 61 | + |
| 62 | + |
| 63 | +```python |
| 64 | +points={'pt0': (0.5,0.5,0.0), 'pt1': (0.2,0.2,0.0)} |
| 65 | +``` |
| 66 | + |
| 67 | + |
| 68 | +```python |
| 69 | +can be retrieved via |
| 70 | +``` |
| 71 | + |
| 72 | + |
| 73 | +```python |
| 74 | +vtufile.getPointData("pressure", pts=points) |
| 75 | +``` |
| 76 | + |
| 77 | + |
| 78 | + |
| 79 | + |
| 80 | + {'pt0': 3.413510714673346e-17, 'pt1': 0.6000000000000001} |
| 81 | + |
| 82 | + |
| 83 | + |
| 84 | +## 1.1 Creating contour plots |
| 85 | + |
| 86 | + |
| 87 | +```python |
| 88 | +import matplotlib.pyplot as plt |
| 89 | +import matplotlib.tri as tri |
| 90 | +``` |
| 91 | + |
| 92 | + |
| 93 | +```python |
| 94 | +vtufile = vtuIO.VTUIO("examples/square2d_random.vtu", dim=2) |
| 95 | +``` |
| 96 | + |
| 97 | + |
| 98 | +```python |
| 99 | +field = vtufile.getField("gaussian_field_2"); |
| 100 | +``` |
| 101 | + |
| 102 | + |
| 103 | +```python |
| 104 | +triang = tri.Triangulation(vtufile.points[:,0], vtufile.points[:,1]) |
| 105 | +``` |
| 106 | + |
| 107 | + |
| 108 | +```python |
| 109 | +plt.tricontourf(triang,field) |
| 110 | +``` |
| 111 | + |
| 112 | + |
| 113 | + |
| 114 | + |
| 115 | + <matplotlib.tri.tricontour.TriContourSet at 0x7fc9e6cc5850> |
| 116 | + |
| 117 | + |
| 118 | + |
| 119 | + |
| 120 | + |
| 121 | + |
| 122 | + |
| 123 | + |
| 124 | + |
| 125 | +### _This random field was created using the ranmedi package:_ https://github.com/joergbuchwald/ranmedi/ |
| 126 | + |
| 127 | +## 1.2 Extracting Pointsetdata |
| 128 | + |
| 129 | +There are basically three interpolation methods available for extracting data at arbitrary points (`cubic` is only available for 1D and 2D). The default is `linear`. |
| 130 | + |
| 131 | + |
| 132 | +```python |
| 133 | +methods = ["nearest", "linear", "cubic"] |
| 134 | +``` |
| 135 | + |
| 136 | + |
| 137 | +```python |
| 138 | +vtufile = vtuIO.VTUIO("examples/square2d_random.vtu", dim=2) |
| 139 | +data_diag = {} |
| 140 | +for method in methods: |
| 141 | + data_diag[method] = vtufile.getPointSetData("gaussian_field_2", pointsetarray=diagonal, interpolation_method=method) |
| 142 | +``` |
| 143 | + |
| 144 | + |
| 145 | +```python |
| 146 | +import numpy as np |
| 147 | +x = np.linspace(0.0,64,num=100); |
| 148 | +``` |
| 149 | + |
| 150 | + |
| 151 | +```python |
| 152 | +diagonal = [(i,i,0) for i in x]; |
| 153 | +``` |
| 154 | + |
| 155 | + |
| 156 | +```python |
| 157 | +r_diag = np.sqrt(2*x*x) |
| 158 | +``` |
| 159 | + |
| 160 | + |
| 161 | +```python |
| 162 | +for method in methods: |
| 163 | + plt.plot(r_diag, data_diag[method], label=method) |
| 164 | +plt.legend() |
| 165 | +``` |
| 166 | + |
| 167 | + |
| 168 | + |
| 169 | + |
| 170 | + <matplotlib.legend.Legend at 0x7fc9dd0106a0> |
| 171 | + |
| 172 | + |
| 173 | + |
| 174 | + |
| 175 | + |
| 176 | + |
| 177 | + |
| 178 | + |
| 179 | + |
| 180 | +# 2. Writing VTU files |
| 181 | +some simple methods also exist for adding new fields to an existing VTU file or save it separately: |
| 182 | + |
| 183 | + |
| 184 | +```python |
| 185 | +vtufile = vtuIO.VTUIO("examples/square_1e2_pcs_0_ts_1_t_1.000000.vtu", dim=2) |
| 186 | +``` |
| 187 | + |
| 188 | + |
| 189 | +```python |
| 190 | +p_size = len(vtufile.getField("pressure")) |
| 191 | +``` |
| 192 | + |
| 193 | + |
| 194 | +```python |
| 195 | +p0 = np.ones(p_size) * 1e6 |
| 196 | +``` |
| 197 | + |
| 198 | + |
| 199 | +```python |
| 200 | +vtufile.writeField(p0, "initialPressure", "mesh_initialpressure.vtu") |
| 201 | +``` |
| 202 | + |
| 203 | + |
| 204 | +```python |
| 205 | +def p_init(x,y,z): |
| 206 | + if x<0.5: |
| 207 | + return -0.5e6 |
| 208 | + else: |
| 209 | + return 0.5e6 |
| 210 | +``` |
| 211 | + |
| 212 | + |
| 213 | +```python |
| 214 | +vtufile.func2Field(p_init, "p_init", "mesh_initialpressure.vtu") |
| 215 | +``` |
| 216 | + |
| 217 | + |
| 218 | +```python |
| 219 | +def null(x,y,z): |
| 220 | + return 0.0 |
| 221 | +``` |
| 222 | + |
| 223 | + |
| 224 | +```python |
| 225 | +vtufile.func2MdimField([p_init,p_init,null,null], "sigma00","mesh_initialpressure.vtu") |
| 226 | +``` |
| 227 | + |
| 228 | +# 3. Reading time-series data from PVD files: |
| 229 | + |
| 230 | +Similar to reading VTU files, it is possible extract time series data from a list of vtufiles given as a PVD file. For extracting grid data at arbitrary points within the mesh, there are two methods available. The stadard method is linear interpolation between cell nodes and the other is the value of the closest node: |
| 231 | + |
| 232 | + |
| 233 | +```python |
| 234 | +pvdfile = vtuIO.PVDIO("examples", "square_1e2_pcs_0.pvd", dim=2) |
| 235 | +``` |
| 236 | + |
| 237 | + examples/square_1e2_pcs_0.pvd |
| 238 | + |
| 239 | + |
| 240 | + |
| 241 | +```python |
| 242 | +time = pvdfile.timesteps |
| 243 | +``` |
| 244 | + |
| 245 | + |
| 246 | +```python |
| 247 | +points={'pt0': (0.3,0.5,0.0), 'pt1': (0.24,0.21,0.0)} |
| 248 | +``` |
| 249 | + |
| 250 | + |
| 251 | +```python |
| 252 | +pressure_linear = pvdfile.readTimeSeries("pressure", points) |
| 253 | +``` |
| 254 | + |
| 255 | + |
| 256 | +```python |
| 257 | +pressure_nearest = pvdfile.readTimeSeries("pressure", points, interpolation_method="nearest") |
| 258 | +``` |
| 259 | + |
| 260 | +As point pt0 is a node in the mesh, both values at $t=1$ agree, whereas pt1 is not a mesh node point resulting in different values. |
| 261 | + |
| 262 | + |
| 263 | +```python |
| 264 | +plt.plot(time, pressure_linear["pt0"], "b-", label="pt0 linear interpolated") |
| 265 | +plt.plot(time, pressure_nearest["pt0"], "b--", label="pt0 closest point value") |
| 266 | +plt.plot(time, pressure_linear["pt1"], "r-", label="pt1 linear interpolated") |
| 267 | +plt.plot(time, pressure_nearest["pt1"], "r--", label="pt1 closest point value") |
| 268 | +plt.legend() |
| 269 | +plt.xlabel("t") |
| 270 | +plt.ylabel("p") |
| 271 | +``` |
| 272 | + |
| 273 | + |
| 274 | + |
| 275 | + |
| 276 | + Text(0, 0.5, 'p') |
| 277 | + |
| 278 | + |
| 279 | + |
| 280 | + |
| 281 | + |
| 282 | + |
| 283 | + |
| 284 | + |
| 285 | + |
| 286 | +# 4. Reading point set data from PVD files |
| 287 | + |
| 288 | +Define two discretized axes: |
| 289 | + |
| 290 | + |
| 291 | +```python |
| 292 | +x = np.linspace(0,1,101) |
| 293 | +``` |
| 294 | + |
| 295 | + |
| 296 | +```python |
| 297 | +xaxis = [(i,0,0) for i in x] |
| 298 | +diagonal = [(i,i,0) for i in x] |
| 299 | +r_diag = np.sqrt(2*x*x) |
| 300 | +``` |
| 301 | + |
| 302 | + |
| 303 | +```python |
| 304 | +t1 = 0.2543 |
| 305 | +t2 = 0.9 |
| 306 | +``` |
| 307 | + |
| 308 | + |
| 309 | +```python |
| 310 | +pressure_xaxis_t1 = pvdfile.readPointSetData(t1, "pressure", pointsetarray=xaxis) |
| 311 | +pressure_diagonal_t1 = pvdfile.readPointSetData(t1, "pressure", pointsetarray=diagonal) |
| 312 | +pressure_xaxis_t2 = pvdfile.readPointSetData(t2, "pressure", pointsetarray=xaxis) |
| 313 | +pressure_diagonal_t2 = pvdfile.readPointSetData(t2, "pressure", pointsetarray=diagonal) |
| 314 | +``` |
| 315 | + |
| 316 | + |
| 317 | +```python |
| 318 | +plt.plot(x, pressure_xaxis_t1, label="p_x t=t1") |
| 319 | +plt.plot(r_diag, pressure_diagonal_t1, label="p_diag t=t1") |
| 320 | +plt.plot(x, pressure_xaxis_t2, label="p_x t=t1") |
| 321 | +plt.plot(r_diag, pressure_diagonal_t2, label="p_diag t=t1") |
| 322 | +plt.xlabel("r") |
| 323 | +plt.ylabel("p") |
| 324 | +plt.legend() |
| 325 | +``` |
| 326 | + |
| 327 | + |
| 328 | + |
| 329 | + |
| 330 | + <matplotlib.legend.Legend at 0x7fc9dc7fadc0> |
| 331 | + |
| 332 | + |
| 333 | + |
| 334 | + |
| 335 | + |
| 336 | + |
| 337 | + |
| 338 | + |
| 339 | + |
| 340 | + |
| 341 | +```python |
| 342 | + |
| 343 | +``` |
0 commit comments