GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
stream-euler.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2019-2023 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 
57 
58 // Coordinates of a point in C-Space (unit square mesh)
59 
60 struct Vector2
61 {
62  double x, y;
63 };
64 
65 // The integer value and the fractional value from a point in C-Space.
66 // Equivalent to the cell index the point is located in and the local
67 // coordinates of the point in the cell.
68 
69 struct Cell2
70 {
71  double fcx, fcy;
72  signed long idx, idy;
73 };
74 
75 struct Vector3
76 {
77  double x, y, z;
78 };
79 
80 struct Cell3
81 {
82  double fcx, fcy, fcz;
83  signed long idx, idy, idz;
84 };
85 
86 static inline void
87 number_to_fractional (signed long *id, double *fc, const double u)
88 {
89  *id = floor (u);
90  *fc = u - *id;
91 }
92 
93 static inline octave_idx_type
95 {
96  return (id < N - 1 ? id : N - 2);
97 }
98 
99 static inline void
100 handle_border (octave_idx_type *id2, double *fc2, const octave_idx_type id1,
101  const double fc1, const octave_idx_type N)
102 {
103  if (id1 < N - 1)
104  {
105  *id2 = id1;
106  *fc2 = fc1;
107  }
108  else
109  {
110  *id2 = N - 2;
111  *fc2 = 1.0;
112  }
113 }
114 
115 static inline double
116 bilinear (const double u11, const double u21, const double u12,
117  const double u22, const double x, const double y)
118 {
119  return (u11 * (1-x) * (1-y) +
120  u21 * x * (1-y) +
121  u12 * (1-x) * y +
122  u22 * x * y);
123 }
124 
125 static inline Cell2
127 {
128  Cell2 Z;
129 
130  number_to_fractional (&Z.idx, &Z.fcx, X.x);
131  number_to_fractional (&Z.idy, &Z.fcy, X.y);
132 
133  return (Z);
134 }
135 
136 static inline bool
138  const octave_idx_type rows)
139 {
140  return ( (((X.idx >= 0) && (X.idx < cols-1)) ||
141  ((X.idx == cols-1) && (X.fcx == 0.0))) &&
142  (((X.idy >= 0) && (X.idy < rows-1)) ||
143  ((X.idy == rows-1) && (X.fcy == 0.0))) );
144 }
145 
146 static inline Vector2
147 add2d (const Cell2 X, const Vector2 Y)
148 {
149  Vector2 Z = {X.idx + X.fcx + Y.x,
150  X.idy + X.fcy + Y.y
151  };
152 
153  return (Z);
154 }
155 
156 static inline Vector2
157 vector_interpolation2d (const Cell2 X, const Matrix& u, const Matrix& v,
158  const octave_idx_type cols, const octave_idx_type rows)
159 {
160  Vector2 V;
161  double fcx, fcy;
162  octave_idx_type idx, idy;
163 
164  handle_border (&idx, &fcx, X.idx, X.fcx, cols);
165  handle_border (&idy, &fcy, X.idy, X.fcy, rows);
166 
167  V.x = bilinear (u(idy, idx), u(idy, idx+1), u(idy+1, idx),
168  u(idy+1, idx+1), fcx, fcy);
169  V.y = bilinear (v(idy, idx), v(idy, idx+1), v(idy+1, idx),
170  v(idy+1, idx+1), fcx, fcy);
171 
172  return (V);
173 }
174 
175 // Apply the Jacobian matrix on the vector V.
176 // The step vector length is set to h.
177 
178 static inline Vector2
180  const RowVector& tx, const RowVector& ty,
181  const octave_idx_type cols, const octave_idx_type rows,
182  const double h)
183 {
184  Vector2 S;
185 
186  const octave_idx_type idx = handle_border_index (X.idx, cols);
187  const octave_idx_type idy = handle_border_index (X.idy, rows);
188 
189  const double x = V.x * tx(idx);
190  const double y = V.y * ty(idy);
191  const double n = 1.0 / sqrt (x*x + y*y);
192  S.x = h * n * x;
193  S.y = h * n * y;
194 
195  return (S);
196 }
197 
198 static inline bool
200 {
201  return ((math::isnan (V.x) || math::isnan (V.y)) ||
202  ((V.x == 0) && (V.y == 0)));
203 }
204 
205 static void
206 euler2d (const octave_idx_type cols, const octave_idx_type rows,
207  const Matrix& u, const Matrix& v,
208  const RowVector& tx, const RowVector& ty,
209  const double zeta, const double xi,
210  const double h, const octave_idx_type maxnverts,
211  Matrix& buffer, octave_idx_type *nverts)
212 {
213  Vector2 V0, V1, S0, X1, Xnxt, S1;
214  const Vector2 X0 = {zeta, xi};
215  Cell2 X0f, X1f;
216 
217  octave_idx_type i = 0;
218 
219  buffer(i, 0) = X0.x;
220  buffer(i, 1) = X0.y;
221 
222  X0f = vector_to_cell2d (X0);
223  while (true)
224  {
225  if (! is_in_definition_set2d (X0f, cols, rows))
226  break;
227 
228  V0 = vector_interpolation2d (X0f, u, v, cols, rows);
229  if (is_singular2d (V0))
230  break;
231 
232  S0 = calculate_step_vector2d (X0f, V0, tx, ty, cols, rows, h);
233 
234  X1 = add2d (X0f, S0);
235  X1f = vector_to_cell2d (X1);
236  if (! is_in_definition_set2d (X1f, cols, rows))
237  break;
238 
239  V1 = vector_interpolation2d (X1f, u, v, cols, rows);
240  if (is_singular2d (V1))
241  break;
242 
243  S1 = calculate_step_vector2d (X1f, V1, tx, ty, cols, rows, h);
244 
245  // Runge Kutta - Heun's Scheme
246  const Vector2 S = {0.5 * (S0.x + S1.x),
247  0.5 * (S0.y + S1.y)
248  };
249  Xnxt = add2d (X0f, S);
250 
251  X0f = vector_to_cell2d (Xnxt);
252  if (! is_in_definition_set2d (X0f, cols, rows))
253  break;
254 
255  i++;
256  buffer(i, 0) = Xnxt.x;
257  buffer(i, 1) = Xnxt.y;
258 
259  if (i + 1 >= maxnverts)
260  break;
261  }
262 
263  *nverts = i + 1;
264 }
265 
266 static inline double
267 trilinear (const double u111, const double u211, const double u121,
268  const double u221, const double u112, const double u212,
269  const double u122, const double u222,
270  const double x, const double y, const double z)
271 {
272  return (u111 * (1-x) * (1-y) * (1-z) +
273  u211 * x * (1-y) * (1-z) +
274  u121 * (1-x) * y * (1-z) +
275  u221 * x * y * (1-z) +
276  u112 * (1-x) * (1-y) * z +
277  u212 * x * (1-y) * z +
278  u122 * (1-x) * y * z +
279  u222 * x * y * z);
280 }
281 
282 static inline Cell3
284 {
285  Cell3 Z;
286 
287  number_to_fractional (&Z.idx, &Z.fcx, X.x);
288  number_to_fractional (&Z.idy, &Z.fcy, X.y);
289  number_to_fractional (&Z.idz, &Z.fcz, X.z);
290 
291  return (Z);
292 }
293 
294 static inline bool
296  const octave_idx_type ny, const octave_idx_type nz)
297 {
298  return ( (((X.idx >= 0) && (X.idx < nx-1)) ||
299  ((X.idx == nx-1) && (X.fcx == 0.0))) &&
300  (((X.idy >= 0) && (X.idy < ny-1)) ||
301  ((X.idy == ny-1) && (X.fcy == 0.0))) &&
302  (((X.idz >= 0) && (X.idz < nz-1)) ||
303  ((X.idz == nz-1) && (X.fcz == 0.0))) );
304 }
305 
306 static inline Vector3
307 add3d (const Cell3 X, const Vector3 Y)
308 {
309  Vector3 Z = {X.idx + X.fcx + Y.x,
310  X.idy + X.fcy + Y.y,
311  X.idz + X.fcz + Y.z
312  };
313 
314  return (Z);
315 }
316 
317 static inline Vector3
318 vector_interpolation3d (const Cell3 X, const NDArray& u, const NDArray& v,
319  const NDArray& w, const octave_idx_type nx,
320  const octave_idx_type ny, const octave_idx_type nz)
321 {
322  Vector3 V;
323  double fcx, fcy, fcz;
324  octave_idx_type idx, idy, idz;
325 
326  handle_border (&idx, &fcx, X.idx, X.fcx, nx);
327  handle_border (&idy, &fcy, X.idy, X.fcy, ny);
328  handle_border (&idz, &fcz, X.idz, X.fcz, nz);
329 
330  V.x = trilinear (u(idy, idx, idz), u(idy, idx+1, idz),
331  u(idy+1, idx, idz), u(idy+1, idx+1, idz),
332  u(idy, idx, idz+1), u(idy, idx+1, idz+1),
333  u(idy+1, idx, idz+1), u(idy+1, idx+1, idz+1),
334  fcx, fcy, fcz);
335  V.y = trilinear (v(idy, idx, idz), v(idy, idx+1, idz),
336  v(idy+1, idx, idz), v(idy+1, idx+1, idz),
337  v(idy, idx, idz+1), v(idy, idx+1, idz+1),
338  v(idy+1, idx, idz+1), v(idy+1, idx+1, idz+1),
339  fcx, fcy, fcz);
340  V.z = trilinear (w(idy, idx, idz), w(idy, idx+1, idz),
341  w(idy+1, idx, idz), w(idy+1, idx+1, idz),
342  w(idy, idx, idz+1), w(idy, idx+1, idz+1),
343  w(idy+1, idx, idz+1), w(idy+1, idx+1, idz+1),
344  fcx, fcy, fcz);
345 
346  return (V);
347 }
348 
349 static inline Vector3
351  const RowVector& tx, const RowVector& ty, const RowVector& tz,
352  const octave_idx_type nx, const octave_idx_type ny,
353  const octave_idx_type nz, const double h)
354 {
355  Vector3 S;
356 
357  const octave_idx_type idx = handle_border_index (X.idx, nx);
358  const octave_idx_type idy = handle_border_index (X.idy, ny);
359  const octave_idx_type idz = handle_border_index (X.idz, nz);
360 
361  const double x = V.x * tx(idx);
362  const double y = V.y * ty(idy);
363  const double z = V.z * tz(idz);
364  const double n = 1.0 / sqrt (x*x + y*y + z*z);
365  S.x = h * n * x;
366  S.y = h * n * y;
367  S.z = h * n * z;
368 
369  return (S);
370 }
371 
372 static inline bool
374 {
375  return ((math::isnan (V.x) || math::isnan (V.y) || math::isnan (V.z)) ||
376  ((V.x == 0) && (V.y == 0) && (V.z == 0)));
377 }
378 
379 static void
381  const NDArray& u, const NDArray& v, const NDArray& w,
382  const RowVector& tx, const RowVector& ty, const RowVector& tz,
383  const double zeta, const double xi, const double rho,
384  const double h, const octave_idx_type maxnverts,
385  Matrix& buffer, octave_idx_type *nverts)
386 {
387  Vector3 V0, V1, S0, X1, Xnxt, S1;
388  const Vector3 X0 = {zeta, xi, rho};
389  Cell3 X0f, X1f;
390 
391  octave_idx_type i = 0;
392  buffer(i, 0) = X0.x;
393  buffer(i, 1) = X0.y;
394  buffer(i, 2) = X0.z;
395 
396  X0f = vector_to_cell3d (X0);
397  while (true)
398  {
399  if (! is_in_definition_set3d (X0f, nx, ny, nz))
400  break;
401 
402  V0 = vector_interpolation3d (X0f, u, v, w, nx, ny, nz);
403  if (is_singular3d (V0))
404  break;
405 
406  S0 = calculate_step_vector3d (X0f, V0, tx, ty, tz, nx, ny, nz, h);
407 
408  X1 = add3d (X0f, S0);
409 
410  X1f = vector_to_cell3d (X1);
411  if (! is_in_definition_set3d (X1f, nx, ny, nz))
412  break;
413 
414  V1 = vector_interpolation3d (X1f, u, v, w, nx, ny, nz);
415  if (is_singular3d (V1))
416  break;
417 
418  S1 = calculate_step_vector3d (X1f, V1, tx, ty, tz, nx, ny, nz, h);
419 
420  // Runge Kutta - Heun's Scheme
421  const Vector3 S = {0.5 * (S0.x + S1.x),
422  0.5 * (S0.y + S1.y),
423  0.5 * (S0.z + S1.z)
424  };
425  Xnxt = add3d (X0f, S);
426 
427  X0f = vector_to_cell3d (Xnxt);
428  if (! is_in_definition_set3d (X0f, nx, ny, nz))
429  break;
430 
431  i++;
432  buffer(i, 0) = Xnxt.x;
433  buffer(i, 1) = Xnxt.y;
434  buffer(i, 2) = Xnxt.z;
435 
436  if (i + 1 >= maxnverts)
437  break;
438 
439  }
440 
441  *nverts = i + 1;
442 }
443 
444 static octave_value
446 {
447 
448  const int nargin = args.length ();
449  if (nargin != 8)
450  print_usage ();
451 
452  const Matrix U = args(0).matrix_value ();
453  const Matrix V = args(1).matrix_value ();
454  const RowVector TX = args(2).row_vector_value ();
455  const RowVector TY = args(3).row_vector_value ();
456  const double zeta = args(4).double_value ();
457  const double xi = args(5).double_value ();
458  const double h = args(6).double_value ();
459  const octave_idx_type maxnverts = args(7).idx_type_value ();
460 
461  const octave_idx_type rows = U.rows ();
462  const octave_idx_type cols = U.columns ();
463 
464  octave_idx_type nverts;
465  Matrix buffer (maxnverts, 2);
466 
467  euler2d (cols, rows, U, V, TX, TY, zeta, xi, h, maxnverts,
468  buffer, &nverts);
469 
470  Matrix xy = buffer.extract (0, 0, nverts-1, 1);
471 
472  return octave_value (xy);
473 }
474 
475 static octave_value
476 streameuler3d_internal (const octave_value_list& args, const char *fcn)
477 {
478 
479  const int nargin = args.length ();
480  if (nargin != 11)
481  print_usage ();
482 
483  const NDArray U = args(0).array_value ();
484  const NDArray V = args(1).array_value ();
485  const NDArray W = args(2).array_value ();
486  const RowVector TX = args(3).row_vector_value ();
487  const RowVector TY = args(4).row_vector_value ();
488  const RowVector TZ = args(5).row_vector_value ();
489  const double zeta = args(6).double_value ();
490  const double xi = args(7).double_value ();
491  const double rho = args(8).double_value ();
492  const double h = args(9).double_value ();
493  const octave_idx_type maxnverts = args(10).idx_type_value ();
494 
495  const dim_vector dims = args(0).dims ();
496  const int ndims = dims.ndims ();
497  if (ndims != 3)
498  error ("%s: dimension must be 3", fcn);
499 
500  octave_idx_type nverts;
501  Matrix buffer (maxnverts, 3);
502 
503  euler3d (dims(1), dims(0), dims(2), U, V, W, TX, TY, TZ, zeta, xi, rho,
504  h, maxnverts, buffer, &nverts);
505 
506  Matrix xyz = buffer.extract (0, 0, nverts-1, 2);
507 
508  return octave_value (xyz);
509 }
510 
511 DEFUN (__streameuler2d__, args, ,
512  doc: /* -*- texinfo -*-
513 @deftypefn {} {@var{output} =} __streameuler2d__ (@var{U}, @var{V}, @var{TX}, @var{TY}, @var{ZETA}, @var{XI}, @var{H}, @var{MAXNVERTS})
514 Calculate the streamline in a vector field @code{[@var{U}, @var{V}]} starting
515 from a seed point at position @code{[@var{ZETA}, @var{XI}]}. The integrator
516 used is Heun's Scheme. The step size can be controlled by @var{H}. The
517 Jacobian matrix can be defined for each grid cell by
518 @code{[@var{TX}, @var{TY}]}.
519 
520 @seealso{streamline, stream2, stream3, __streameuler3d__}
521 @end deftypefn */)
522 {
523  return streameuler2d_internal (args);
524 }
525 
526 DEFUN (__streameuler3d__, args, ,
527  doc: /* -*- texinfo -*-
528 @deftypefn {} {@var{output} =} __streameuler3d__ (@var{U}, @var{V}, @var{W}, @var{TX}, @var{TY}, @var{TZ}, @var{ZETA}, @var{XI}, @var{RHO}, @var{H}, @var{MAXNVERTS})
529 Calculate the streamline in a vector field @code{[@var{U}, @var{V}, @var{W}]}
530 starting from a seed point at position
531 @code{[@var{ZETA}, @var{XI}, @var{RHO}]}. The integrator used is Heun's
532 Scheme. The step size can be controlled by @var{H}. The Jacobian matrix can
533 be defined for each grid cell by @code{[@var{TX}, @var{TY}, @var{TZ}]}.
534 
535 @seealso{streamline, stream2, stream3, __streameuler2d__}
536 @end deftypefn */)
537 {
538  return streameuler3d_internal (args, "__streameuler3d__");
539 }
540 
OCTAVE_END_NAMESPACE(octave)
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type columns(void) const
Definition: Array.h:471
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
Definition: Array.h:459
Definition: dMatrix.h:42
OCTAVE_API Matrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: dMatrix.cc:397
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:257
octave_idx_type length(void) const
Definition: ovl.h:113
Array< octave_value > array_value(void) const
Definition: ovl.h:90
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
void error(const char *fmt,...)
Definition: error.cc:979
F77_RET_T const F77_INT const F77_INT const F77_INT const F77_DBLE const F77_DBLE F77_INT F77_DBLE * V
F77_RET_T const F77_INT & N
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Z
bool isnan(bool)
Definition: lo-mappers.h:178
std::complex< T > floor(const std::complex< T > &x)
Definition: lo-mappers.h:130
F77_RET_T const F77_DBLE * x
octave_idx_type n
Definition: mx-inlines.cc:753
std::complex< double > w(std::complex< double > z, double relerr=0)
return octave_value(v1.char_array_value() . concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string()) ? '\'' :'"'))
static const double xi[33]
Definition: quadcc.cc:68
static void handle_border(octave_idx_type *id2, double *fc2, const octave_idx_type id1, const double fc1, const octave_idx_type N)
static bool is_in_definition_set2d(const Cell2 X, const octave_idx_type cols, const octave_idx_type rows)
static Cell2 vector_to_cell2d(const Vector2 X)
static octave_value streameuler3d_internal(const octave_value_list &args, const char *fcn)
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)
static Cell3 vector_to_cell3d(const Vector3 X)
static void number_to_fractional(signed long *id, double *fc, const double u)
Definition: stream-euler.cc:87
static bool is_in_definition_set3d(const Cell3 X, const octave_idx_type nx, const octave_idx_type ny, const octave_idx_type nz)
static octave_value streameuler2d_internal(const octave_value_list &args)
static octave_idx_type handle_border_index(const octave_idx_type id, const octave_idx_type N)
Definition: stream-euler.cc:94
static bool is_singular2d(const Vector2 V)
static Vector3 add3d(const Cell3 X, const Vector3 Y)
static Vector2 add2d(const Cell2 X, const Vector2 Y)
static Vector2 vector_interpolation2d(const Cell2 X, const Matrix &u, const Matrix &v, const octave_idx_type cols, const octave_idx_type rows)
static 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)
static 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)
static 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)
static 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)
static bool is_singular3d(const Vector3 V)
static double bilinear(const double u11, const double u21, const double u12, const double u22, const double x, const double y)
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)
signed long idy
Definition: stream-euler.cc:72
double fcx
Definition: stream-euler.cc:71
signed long idx
Definition: stream-euler.cc:72
double fcy
Definition: stream-euler.cc:71
double fcx
Definition: stream-euler.cc:82
signed long idz
Definition: stream-euler.cc:83
signed long idx
Definition: stream-euler.cc:83
double fcz
Definition: stream-euler.cc:82
signed long idy
Definition: stream-euler.cc:83
double fcy
Definition: stream-euler.cc:82
double y
Definition: stream-euler.cc:62
double x
Definition: stream-euler.cc:62
double x
Definition: stream-euler.cc:77
double z
Definition: stream-euler.cc:77
double y
Definition: stream-euler.cc:77