diff --git a/OpenSim/Region/Terrain.BasicTerrain/libTerrainBSD/Channel/Manipulators/NavierStokes.cs b/OpenSim/Region/Terrain.BasicTerrain/libTerrainBSD/Channel/Manipulators/NavierStokes.cs new file mode 100644 index 0000000000..3eae46b025 --- /dev/null +++ b/OpenSim/Region/Terrain.BasicTerrain/libTerrainBSD/Channel/Manipulators/NavierStokes.cs @@ -0,0 +1,214 @@ +/* +* Copyright (c) Contributors, http://www.openmetaverse.org/ +* See CONTRIBUTORS.TXT for a full list of copyright holders. +* +* Redistribution and use in source and binary forms, with or without +* modification, are permitted provided that the following conditions are met: +* * Redistributions of source code must retain the above copyright +* notice, this list of conditions and the following disclaimer. +* * Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions and the following disclaimer in the +* documentation and/or other materials provided with the distribution. +* * Neither the name of the OpenSim Project nor the +* names of its contributors may be used to endorse or promote products +* derived from this software without specific prior written permission. +* +* THIS SOFTWARE IS PROVIDED BY THE DEVELOPERS ``AS IS AND ANY +* EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +* DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY +* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +* +*/ + +using System; +using System.Collections.Generic; +using System.Text; + +namespace libTerrain +{ + partial class Channel + { + // Navier Stokes Algorithms ported from + // "Real-Time Fluid Dynamics for Games" by Jos Stam. + // presented at GDC 2003. + + // Poorly ported from C++. (I gave up making it properly native somewhere after nsSetBnd) + + private static int nsIX(int i, int j, int N) + { + return ((i) + (N + 2) * (j)); + } + + private static void nsSwap(ref double x0, ref double x) + { + double tmp = x0; + x0 = x; + x = tmp; + } + + private static void nsSwap(ref double[] x0, ref double[] x) + { + double[] tmp = x0; + x0 = x; + x = tmp; + } + + private void nsAddSource(int N, ref double[] x, ref double[] s, double dt) + { + int i; + int size = (N + 2) * (N + 2); + for (i = 0; i < size; i++) + { + x[i] += dt * s[i]; + } + } + + private void nsSetBnd(int N, int b, ref double[] x) + { + int i; + for (i = 0; i <= N; i++) + { + x[nsIX(0, i, N)] = b == 1 ? -x[nsIX(1, i, N)] : x[nsIX(1, i, N)]; + x[nsIX(0, N + 1, N)] = b == 1 ? -x[nsIX(N, i, N)] : x[nsIX(N, i, N)]; + x[nsIX(i, 0, N)] = b == 2 ? -x[nsIX(i, 1, N)] : x[nsIX(i, 1, N)]; + x[nsIX(i, N + 1, N)] = b == 2 ? -x[nsIX(i, N, N)] : x[nsIX(i, N, N)]; + } + x[nsIX(0, 0, N)] = 0.5f * (x[nsIX(1, 0, N)] + x[nsIX(0, 1, N)]); + x[nsIX(0, N + 1, N)] = 0.5f * (x[nsIX(1, N + 1, N)] + x[nsIX(0, N, N)]); + x[nsIX(N + 1, 0, N)] = 0.5f * (x[nsIX(N, 0, N)] + x[nsIX(N + 1, 1, N)]); + x[nsIX(N + 1, N + 1, N)] = 0.5f * (x[nsIX(N, N + 1, N)] + x[nsIX(N + 1, N, N)]); + } + + private void nsLinSolve(int N, int b, ref double[] x, ref double[] x0, double a, double c) + { + int i, j; + for (i = 1; i <= N; i++) + { + for (j = 1; j <= N; j++) + { + x[nsIX(i, j, N)] = (x0[nsIX(i, j, N)] + a * + (x[nsIX(i - 1, j, N)] + + x[nsIX(i + 1, j, N)] + + x[nsIX(i, j - 1, N)] + x[nsIX(i, j + 1, N)]) + ) / c; + } + } + + nsSetBnd(N, b, ref x); + } + + private void nsDiffuse(int N, int b, ref double[] x, ref double[] x0, double diff, double dt) + { + double a = dt * diff * N * N; + nsLinSolve(N, b, ref x, ref x0, a, 1 + 4 * a); + } + + private void nsAdvect(int N, int b, ref double[] d, ref double[] d0, ref double[] u, ref double[] v, double dt) + { + int i, j, i0, j0, i1, j1; + double x, y, s0, t0, s1, t1, dt0; + + dt0 = dt * N; + + for (i = 1; i <= N; i++) + { + for (j = 1; j <= N; j++) + { + x = i - dt0 * u[nsIX(i, j, N)]; + y = j - dt0 * v[nsIX(i, j, N)]; + + if (x < 0.5) + x = 0.5; + if (x > N + 0.5) + x = N + 0.5; + i0 = (int)x; + i1 = i0 + 1; + + if (y < 0.5) + y = 0.5; + if (y > N + 0.5) + y = N + 0.5; + j0 = (int)y; + j1 = j0 + 1; + + s1 = x - i0; + s0 = 1 - s1; + t1 = y - j0; + t0 = 1 - t1; + + d[nsIX(i, j, N)] = s0 * (t0 * d0[nsIX(i0, j0, N)] + t1 * d0[nsIX(i0, j1, N)]) + + s1 * (t0 * d0[nsIX(i1, j0, N)] + t1 * d0[nsIX(i1, j1, N)]); + } + } + + nsSetBnd(N, b, ref d); + } + + public void nsProject(int N, ref double[] u, ref double[] v, ref double[] p, ref double[] div) + { + int i, j; + + for (i = 1; i <= N; i++) + { + for (j = 1; j <= N; j++) + { + div[nsIX(i, j, N)] = -0.5 * (u[nsIX(i + 1, j, N)] - u[nsIX(i - 1, j, N)] + v[nsIX(i, j + 1, N)] - v[nsIX(i, j - 1, N)]) / N; + p[nsIX(i, j, N)] = 0; + } + } + + nsSetBnd(N, 0, ref div); + nsSetBnd(N, 0, ref p); + + nsLinSolve(N, 0, ref p, ref div, 1, 4); + + for (i = 1; i <= N; i++) + { + for (j = 1; j <= N; j++) + { + u[nsIX(i, j, N)] -= 0.5 * N * (p[nsIX(i + 1, j, N)] - p[nsIX(i - 1, j, N)]); + v[nsIX(i, j, N)] -= 0.5 * N * (p[nsIX(i, j + 1, N)] - p[nsIX(i, j - 1, N)]); + } + } + + nsSetBnd(N, 1, ref u); + nsSetBnd(N, 2, ref v); + } + + private void nsDensStep(int N, ref double[] x, ref double[] x0, ref double[] u, ref double[] v, double diff, double dt) + { + nsAddSource(N, ref x, ref x0, dt); + nsSwap(ref x0, ref x); + nsDiffuse(N, 0, ref x, ref x0, diff, dt); + nsSwap(ref x0, ref x); + nsAdvect(N, 0, ref x, ref x0, ref u, ref v, dt); + } + + private void nsVelStep(int N, ref double[] u, ref double[] v, ref double[] u0, ref double[] v0, double visc, double dt) + { + nsAddSource(N, ref u, ref u0, dt); + nsAddSource(N, ref v, ref v0, dt); + nsSwap(ref u0, ref u); + nsDiffuse(N, 1, ref u, ref u0, visc, dt); + nsSwap(ref v0, ref v); + nsDiffuse(N, 2, ref v, ref v0, visc, dt); + nsProject(N, ref u, ref v, ref u0, ref v0); + nsSwap(ref u0, ref u); + nsSwap(ref v0, ref v); + nsAdvect(N, 1, ref u, ref u0, ref u0, ref v0, dt); + nsAdvect(N, 2, ref v, ref v0, ref u0, ref v0, dt); + nsProject(N, ref u, ref v, ref u0, ref v0); + } + + public void navierSimulate() + { + + } + } +} \ No newline at end of file