1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2019-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 /*
27
28 References:
29
30 @article{
31 title = {Particle Tracing Algorithms for 3D Curvilinear Grids},
32 year = {2000},
33 author = {Nielson, Gregory and Uller, H. and Sadarjoen, I. and Walsum, Theo and Hin, Andrea and Post, Frits}
34 }
35
36 @article{
37 title = {Sources of error in the graphical analysis of CFD results},
38 publisher = {Journal of Scientific Computing},
39 year = {1988},
40 volume = {3},
41 number = {2},
42 pages = {149--164},
43 author = {Buning, Pieter G.},
44 }
45
46 */
47
48 #if defined (HAVE_CONFIG_H)
49 # include "config.h"
50 #endif
51
52 #include "defun.h"
53 #include "error.h"
54 #include "ovl.h"
55
56 // Coordinates of a point in C-Space (unit square mesh)
57
58 typedef struct
59 {
60 double x, y;
61 } Vector2;
62
63 // The integer value and the fractional value from a point in C-Space.
64 // Equivalent to the cell index the point is located in and the local
65 // coordinates of the point in the cell.
66
67 typedef struct
68 {
69 double fcx, fcy;
70 signed long idx, idy;
71 } Cell2;
72
73 typedef struct
74 {
75 double x, y, z;
76 } Vector3;
77
78 typedef struct
79 {
80 double fcx, fcy, fcz;
81 signed long idx, idy, idz;
82 } Cell3;
83
84 static inline void
number_to_fractional(signed long * id,double * fc,const double u)85 number_to_fractional (signed long *id, double *fc, const double u)
86 {
87 *id = floor (u);
88 *fc = u - *id;
89 }
90
91 static inline octave_idx_type
handle_border_index(const octave_idx_type id,const octave_idx_type N)92 handle_border_index (const octave_idx_type id, const octave_idx_type N)
93 {
94 return (id < N - 1 ? id : N - 2);
95 }
96
97 static inline void
handle_border(octave_idx_type * id2,double * fc2,const octave_idx_type id1,const double fc1,const octave_idx_type N)98 handle_border (octave_idx_type *id2, double *fc2, const octave_idx_type id1,
99 const double fc1, const octave_idx_type N)
100 {
101 if (id1 < N - 1)
102 {
103 *id2 = id1;
104 *fc2 = fc1;
105 }
106 else
107 {
108 *id2 = N - 2;
109 *fc2 = 1.0;
110 }
111 }
112
113 static inline double
bilinear(const double u11,const double u21,const double u12,const double u22,const double x,const double y)114 bilinear (const double u11, const double u21, const double u12,
115 const double u22, const double x, const double y)
116 {
117 return (u11 * (1-x) * (1-y) +
118 u21 * x * (1-y) +
119 u12 * (1-x) * y +
120 u22 * x * y);
121 }
122
123 static inline Cell2
vector_to_cell2d(const Vector2 X)124 vector_to_cell2d (const Vector2 X)
125 {
126 Cell2 Z;
127
128 number_to_fractional (&Z.idx, &Z.fcx, X.x);
129 number_to_fractional (&Z.idy, &Z.fcy, X.y);
130
131 return (Z);
132 }
133
134 static inline bool
is_in_definition_set2d(const Cell2 X,const octave_idx_type cols,const octave_idx_type rows)135 is_in_definition_set2d (const Cell2 X, const octave_idx_type cols,
136 const octave_idx_type rows)
137 {
138 return ( (((X.idx >= 0) && (X.idx < cols-1)) ||
139 ((X.idx == cols-1) && (X.fcx == 0.0))) &&
140 (((X.idy >= 0) && (X.idy < rows-1)) ||
141 ((X.idy == rows-1) && (X.fcy == 0.0))) );
142 }
143
144 static inline Vector2
add2d(const Cell2 X,const Vector2 Y)145 add2d (const Cell2 X, const Vector2 Y)
146 {
147 Vector2 Z = {X.idx + X.fcx + Y.x,
148 X.idy + X.fcy + Y.y};
149
150 return (Z);
151 }
152
153 static inline Vector2
vector_interpolation2d(const Cell2 X,const Matrix & u,const Matrix & v,const octave_idx_type cols,const octave_idx_type rows)154 vector_interpolation2d (const Cell2 X, const Matrix& u, const Matrix& v,
155 const octave_idx_type cols, const octave_idx_type rows)
156 {
157 Vector2 V;
158 double fcx, fcy;
159 octave_idx_type idx, idy;
160
161 handle_border (&idx, &fcx, X.idx, X.fcx, cols);
162 handle_border (&idy, &fcy, X.idy, X.fcy, rows);
163
164 V.x = bilinear (u(idy, idx), u(idy, idx+1), u(idy+1, idx),
165 u(idy+1, idx+1), fcx, fcy);
166 V.y = bilinear (v(idy, idx), v(idy, idx+1), v(idy+1, idx),
167 v(idy+1, idx+1), fcx, fcy);
168
169 return (V);
170 }
171
172 // Apply the Jacobian matrix on the vector V.
173 // The step vector length is set to h.
174
175 static inline Vector2
calculate_step_vector2d(const Cell2 X,const Vector2 V,const RowVector & tx,const RowVector & ty,const octave_idx_type cols,const octave_idx_type rows,const double h)176 calculate_step_vector2d (const Cell2 X, const Vector2 V,
177 const RowVector& tx, const RowVector& ty,
178 const octave_idx_type cols, const octave_idx_type rows,
179 const double h)
180 {
181 Vector2 S;
182
183 const octave_idx_type idx = handle_border_index (X.idx, cols);
184 const octave_idx_type idy = handle_border_index (X.idy, rows);
185
186 const double x = V.x * tx(idx);
187 const double y = V.y * ty(idy);
188 const double n = 1.0 / sqrt (x*x + y*y);
189 S.x = h * n * x;
190 S.y = h * n * y;
191
192 return (S);
193 }
194
195 static inline bool
is_singular2d(const Vector2 V)196 is_singular2d (const Vector2 V)
197 {
198 return ((octave::math::isnan (V.x) || octave::math::isnan (V.y)) ||
199 ((V.x == 0) && (V.y == 0)));
200 }
201
202 static void
euler2d(const octave_idx_type cols,const octave_idx_type rows,const Matrix & u,const Matrix & v,const RowVector & tx,const RowVector & ty,const double zeta,const double xi,const double h,const octave_idx_type maxnverts,Matrix & buffer,octave_idx_type * nverts)203 euler2d (const octave_idx_type cols, const octave_idx_type rows,
204 const Matrix& u, const Matrix& v,
205 const RowVector& tx, const RowVector& ty,
206 const double zeta, const double xi,
207 const double h, const octave_idx_type maxnverts,
208 Matrix& buffer, octave_idx_type *nverts)
209 {
210 Vector2 V0, V1, S0, X1, Xnxt, S1;
211 const Vector2 X0 = {zeta, xi};
212 Cell2 X0f, X1f;
213
214 octave_idx_type i = 0;
215
216 buffer(i, 0) = X0.x;
217 buffer(i, 1) = X0.y;
218
219 X0f = vector_to_cell2d (X0);
220 while (true)
221 {
222 if (! is_in_definition_set2d (X0f, cols, rows))
223 break;
224
225 V0 = vector_interpolation2d (X0f, u, v, cols, rows);
226 if (is_singular2d (V0))
227 break;
228
229 S0 = calculate_step_vector2d (X0f, V0, tx, ty, cols, rows, h);
230
231 X1 = add2d (X0f, S0);
232 X1f = vector_to_cell2d (X1);
233 if (! is_in_definition_set2d (X1f, cols, rows))
234 break;
235
236 V1 = vector_interpolation2d (X1f, u, v, cols, rows);
237 if (is_singular2d (V1))
238 break;
239
240 S1 = calculate_step_vector2d (X1f, V1, tx, ty, cols, rows, h);
241
242 // Runge Kutta - Heun's Scheme
243 const Vector2 S = {0.5 * (S0.x + S1.x),
244 0.5 * (S0.y + S1.y)};
245 Xnxt = add2d (X0f, S);
246
247 X0f = vector_to_cell2d (Xnxt);
248 if (! is_in_definition_set2d (X0f, cols, rows))
249 break;
250
251 i++;
252 buffer(i, 0) = Xnxt.x;
253 buffer(i, 1) = Xnxt.y;
254
255 if (i + 1 >= maxnverts)
256 break;
257 }
258
259 *nverts = i + 1;
260 }
261
262 static inline double
trilinear(const double u111,const double u211,const double u121,const double u221,const double u112,const double u212,const double u122,const double u222,const double x,const double y,const double z)263 trilinear (const double u111, const double u211, const double u121,
264 const double u221, const double u112, const double u212,
265 const double u122, const double u222,
266 const double x, const double y, const double z)
267 {
268 return (u111 * (1-x) * (1-y) * (1-z) +
269 u211 * x * (1-y) * (1-z) +
270 u121 * (1-x) * y * (1-z) +
271 u221 * x * y * (1-z) +
272 u112 * (1-x) * (1-y) * z +
273 u212 * x * (1-y) * z +
274 u122 * (1-x) * y * z +
275 u222 * x * y * z);
276 }
277
278 static inline Cell3
vector_to_cell3d(const Vector3 X)279 vector_to_cell3d (const Vector3 X)
280 {
281 Cell3 Z;
282
283 number_to_fractional (&Z.idx, &Z.fcx, X.x);
284 number_to_fractional (&Z.idy, &Z.fcy, X.y);
285 number_to_fractional (&Z.idz, &Z.fcz, X.z);
286
287 return (Z);
288 }
289
290 static inline bool
is_in_definition_set3d(const Cell3 X,const octave_idx_type nx,const octave_idx_type ny,const octave_idx_type nz)291 is_in_definition_set3d (const Cell3 X, const octave_idx_type nx,
292 const octave_idx_type ny, const octave_idx_type nz)
293 {
294 return ( (((X.idx >= 0) && (X.idx < nx-1)) ||
295 ((X.idx == nx-1) && (X.fcx == 0.0))) &&
296 (((X.idy >= 0) && (X.idy < ny-1)) ||
297 ((X.idy == ny-1) && (X.fcy == 0.0))) &&
298 (((X.idz >= 0) && (X.idz < nz-1)) ||
299 ((X.idz == nz-1) && (X.fcz == 0.0))) );
300 }
301
302 static inline Vector3
add3d(const Cell3 X,const Vector3 Y)303 add3d (const Cell3 X, const Vector3 Y)
304 {
305 Vector3 Z = {X.idx + X.fcx + Y.x,
306 X.idy + X.fcy + Y.y,
307 X.idz + X.fcz + Y.z};
308
309 return (Z);
310 }
311
312 static inline Vector3
vector_interpolation3d(const Cell3 X,const NDArray & u,const NDArray & v,const NDArray & w,const octave_idx_type nx,const octave_idx_type ny,const octave_idx_type nz)313 vector_interpolation3d (const Cell3 X, const NDArray& u, const NDArray& v,
314 const NDArray& w, const octave_idx_type nx,
315 const octave_idx_type ny, const octave_idx_type nz)
316 {
317 Vector3 V;
318 double fcx, fcy, fcz;
319 octave_idx_type idx, idy, idz;
320
321 handle_border (&idx, &fcx, X.idx, X.fcx, nx);
322 handle_border (&idy, &fcy, X.idy, X.fcy, ny);
323 handle_border (&idz, &fcz, X.idz, X.fcz, nz);
324
325 V.x = trilinear (u(idy, idx, idz), u(idy, idx+1, idz),
326 u(idy+1, idx, idz), u(idy+1, idx+1, idz),
327 u(idy, idx, idz+1), u(idy, idx+1, idz+1),
328 u(idy+1, idx, idz+1), u(idy+1, idx+1, idz+1),
329 fcx, fcy, fcz);
330 V.y = trilinear (v(idy, idx, idz), v(idy, idx+1, idz),
331 v(idy+1, idx, idz), v(idy+1, idx+1, idz),
332 v(idy, idx, idz+1), v(idy, idx+1, idz+1),
333 v(idy+1, idx, idz+1), v(idy+1, idx+1, idz+1),
334 fcx, fcy, fcz);
335 V.z = trilinear (w(idy, idx, idz), w(idy, idx+1, idz),
336 w(idy+1, idx, idz), w(idy+1, idx+1, idz),
337 w(idy, idx, idz+1), w(idy, idx+1, idz+1),
338 w(idy+1, idx, idz+1), w(idy+1, idx+1, idz+1),
339 fcx, fcy, fcz);
340
341 return (V);
342 }
343
344 static inline Vector3
calculate_step_vector3d(const Cell3 X,const Vector3 V,const RowVector & tx,const RowVector & ty,const RowVector & tz,const octave_idx_type nx,const octave_idx_type ny,const octave_idx_type nz,const double h)345 calculate_step_vector3d (const Cell3 X, const Vector3 V,
346 const RowVector& tx, const RowVector& ty, const RowVector& tz,
347 const octave_idx_type nx, const octave_idx_type ny,
348 const octave_idx_type nz, const double h)
349 {
350 Vector3 S;
351
352 const octave_idx_type idx = handle_border_index (X.idx, nx);
353 const octave_idx_type idy = handle_border_index (X.idy, ny);
354 const octave_idx_type idz = handle_border_index (X.idz, nz);
355
356 const double x = V.x * tx(idx);
357 const double y = V.y * ty(idy);
358 const double z = V.z * tz(idz);
359 const double n = 1.0 / sqrt (x*x + y*y + z*z);
360 S.x = h * n * x;
361 S.y = h * n * y;
362 S.z = h * n * z;
363
364 return (S);
365 }
366
367 static inline bool
is_singular3d(const Vector3 V)368 is_singular3d (const Vector3 V)
369 {
370 return ((octave::math::isnan (V.x) || octave::math::isnan (V.y) ||
371 octave::math::isnan (V.z)) ||
372 ((V.x == 0) && (V.y == 0) && (V.z == 0)));
373 }
374
375 static void
euler3d(const octave_idx_type nx,const octave_idx_type ny,const octave_idx_type nz,const NDArray & u,const NDArray & v,const NDArray & w,const RowVector & tx,const RowVector & ty,const RowVector & tz,const double zeta,const double xi,const double rho,const double h,const octave_idx_type maxnverts,Matrix & buffer,octave_idx_type * nverts)376 euler3d (const octave_idx_type nx, const octave_idx_type ny, const octave_idx_type nz,
377 const NDArray& u, const NDArray& v, const NDArray& w,
378 const RowVector& tx, const RowVector& ty, const RowVector& tz,
379 const double zeta, const double xi, const double rho,
380 const double h, const octave_idx_type maxnverts,
381 Matrix& buffer, octave_idx_type *nverts)
382 {
383 Vector3 V0, V1, S0, X1, Xnxt, S1;
384 const Vector3 X0 = {zeta, xi, rho};
385 Cell3 X0f, X1f;
386
387 octave_idx_type i = 0;
388 buffer(i, 0) = X0.x;
389 buffer(i, 1) = X0.y;
390 buffer(i, 2) = X0.z;
391
392 X0f = vector_to_cell3d (X0);
393 while (true)
394 {
395 if (! is_in_definition_set3d (X0f, nx, ny, nz))
396 break;
397
398 V0 = vector_interpolation3d (X0f, u, v, w, nx, ny, nz);
399 if (is_singular3d (V0))
400 break;
401
402 S0 = calculate_step_vector3d (X0f, V0, tx, ty, tz, nx, ny, nz, h);
403
404 X1 = add3d (X0f, S0);
405
406 X1f = vector_to_cell3d (X1);
407 if (! is_in_definition_set3d (X1f, nx, ny, nz))
408 break;
409
410 V1 = vector_interpolation3d (X1f, u, v, w, nx, ny, nz);
411 if (is_singular3d (V1))
412 break;
413
414 S1 = calculate_step_vector3d (X1f, V1, tx, ty, tz, nx, ny, nz, h);
415
416 // Runge Kutta - Heun's Scheme
417 const Vector3 S = {0.5 * (S0.x + S1.x),
418 0.5 * (S0.y + S1.y),
419 0.5 * (S0.z + S1.z)};
420 Xnxt = add3d (X0f, S);
421
422 X0f = vector_to_cell3d (Xnxt);
423 if (! is_in_definition_set3d (X0f, nx, ny, nz))
424 break;
425
426 i++;
427 buffer(i, 0) = Xnxt.x;
428 buffer(i, 1) = Xnxt.y;
429 buffer(i, 2) = Xnxt.z;
430
431 if (i + 1 >= maxnverts)
432 break;
433
434 }
435
436 *nverts = i + 1;
437 }
438
439 static octave_value
streameuler2d_internal(const octave_value_list & args)440 streameuler2d_internal (const octave_value_list& args)
441 {
442
443 const int nargin = args.length ();
444 if (nargin != 8)
445 print_usage ();
446
447 const Matrix U = args(0).matrix_value ();
448 const Matrix V = args(1).matrix_value ();
449 const RowVector TX = args(2).row_vector_value ();
450 const RowVector TY = args(3).row_vector_value ();
451 const double zeta = args(4).double_value ();
452 const double xi = args(5).double_value ();
453 const double h = args(6).double_value ();
454 const octave_idx_type maxnverts = args(7).idx_type_value ();
455
456 const octave_idx_type rows = U.rows ();
457 const octave_idx_type cols = U.columns ();
458
459 octave_idx_type nverts;
460 Matrix buffer (maxnverts, 2);
461
462 euler2d (cols, rows, U, V, TX, TY, zeta, xi, h, maxnverts,
463 buffer, &nverts);
464
465 Matrix xy = buffer.extract (0, 0, nverts-1, 1);
466
467 return octave_value (xy);
468 }
469
470 static octave_value
streameuler3d_internal(const octave_value_list & args,const char * fcn)471 streameuler3d_internal (const octave_value_list& args, const char *fcn)
472 {
473
474 const int nargin = args.length ();
475 if (nargin != 11)
476 print_usage ();
477
478 const NDArray U = args(0).array_value ();
479 const NDArray V = args(1).array_value ();
480 const NDArray W = args(2).array_value ();
481 const RowVector TX = args(3).row_vector_value ();
482 const RowVector TY = args(4).row_vector_value ();
483 const RowVector TZ = args(5).row_vector_value ();
484 const double zeta = args(6).double_value ();
485 const double xi = args(7).double_value ();
486 const double rho = args(8).double_value ();
487 const double h = args(9).double_value ();
488 const octave_idx_type maxnverts = args(10).idx_type_value ();
489
490 const dim_vector dims = args(0).dims ();
491 const int ndims = dims.ndims ();
492 if (ndims != 3)
493 error ("%s: dimension must be 3", fcn);
494
495 octave_idx_type nverts;
496 Matrix buffer (maxnverts, 3);
497
498 euler3d (dims(1), dims(0), dims(2), U, V, W, TX, TY, TZ, zeta, xi, rho,
499 h, maxnverts, buffer, &nverts);
500
501 Matrix xyz = buffer.extract (0, 0, nverts-1, 2);
502
503 return octave_value (xyz);
504 }
505
506 DEFUN (__streameuler2d__, args, ,
507 doc: /* -*- texinfo -*-
508 @deftypefn {} {} __streameuler2d__ (@var{U}, @var{V}, @var{TX}, @var{TY}, @var{ZETA}, @var{XI}, @var{H}, @var{MAXNVERTS})
509 Calculates the streamline in a vector field @code{[@var{U}, @var{V}]} starting
510 from a seed point at position @code{[@var{ZETA}, @var{XI}]}. The integrator
511 used is Heun's Scheme. The step size can be controlled by @var{H}. The
512 Jacobian matrix can be defined for each grid cell by
513 @code{[@var{TX}, @var{TY}]}.
514
515 @seealso{streamline, stream2, stream3, __streameuler3d__}
516 @end deftypefn */)
517 {
518 return streameuler2d_internal (args);
519 }
520
521 DEFUN (__streameuler3d__, args, ,
522 doc: /* -*- texinfo -*-
523 @deftypefn {} {} __streameuler3d__ (@var{U}, @var{V}, @var{W}, @var{TX}, @var{TY}, @var{TZ}, @var{ZETA}, @var{XI}, @var{RHO}, @var{H}, @var{MAXNVERTS})
524 Calculates the streamline in a vector field @code{[@var{U}, @var{V}, @var{W}]}
525 starting from a seed point at position
526 @code{[@var{ZETA}, @var{XI}, @var{RHO}]}. The integrator used is Heun's
527 Scheme. The step size can be controlled by @var{H}. The Jacobian matrix can
528 be defined for each grid cell by @code{[@var{TX}, @var{TY}, @var{TZ}]}.
529
530 @seealso{streamline, stream2, stream3, __streameuler2d__}
531 @end deftypefn */)
532 {
533 return streameuler3d_internal (args, "__streameuler3d__");
534 }
535
536