48 #if defined (HAVE_CONFIG_H)
96 return (
id <
N - 1 ?
id :
N - 2);
116 bilinear (
const double u11,
const double u21,
const double u12,
117 const double u22,
const double x,
const double y)
119 return (u11 * (1-
x) * (1-y) +
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))) );
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);
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);
202 ((
V.x == 0) && (
V.y == 0)));
209 const double zeta,
const double xi,
213 Vector2 V0, V1, S0, X1, Xnxt, S1;
234 X1 =
add2d (X0f, S0);
249 Xnxt =
add2d (X0f, S);
256 buffer(i, 0) = Xnxt.
x;
257 buffer(i, 1) = Xnxt.
y;
259 if (i + 1 >= maxnverts)
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)
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 +
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))) );
323 double fcx, fcy, fcz;
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),
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),
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),
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);
376 ((
V.x == 0) && (
V.y == 0) && (
V.z == 0)));
383 const double zeta,
const double xi,
const double rho,
387 Vector3 V0, V1, S0, X1, Xnxt, S1;
408 X1 =
add3d (X0f, S0);
425 Xnxt =
add3d (X0f, S);
432 buffer(i, 0) = Xnxt.
x;
433 buffer(i, 1) = Xnxt.
y;
434 buffer(i, 2) = Xnxt.
z;
436 if (i + 1 >= maxnverts)
448 const int nargin = args.
length ();
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 ();
465 Matrix buffer (maxnverts, 2);
467 euler2d (cols, rows, U,
V, TX, TY, zeta,
xi, h, maxnverts,
479 const int nargin = args.
length ();
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 ();
496 const int ndims = dims.
ndims ();
498 error (
"%s: dimension must be 3", fcn);
501 Matrix buffer (maxnverts, 3);
503 euler3d (dims(1), dims(0), dims(2), U,
V, W, TX, TY, TZ, zeta,
xi, rho,
504 h, maxnverts, buffer, &nverts);
511 DEFUN (__streameuler2d__, args, ,
526 DEFUN (__streameuler3d__, args, ,
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type columns(void) const
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
OCTAVE_API Matrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Vector representing the dimensions (size) of an Array.
octave_idx_type ndims(void) const
Number of dimensions.
octave_idx_type length(void) const
Array< octave_value > array_value(void) const
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
OCTINTERP_API void print_usage(void)
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
void error(const char *fmt,...)
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
std::complex< T > floor(const std::complex< T > &x)
F77_RET_T const F77_DBLE * x
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]
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)
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)
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)