Cube and the flow

Some time  ago I've implemented an interesting technique that extends the standard Lattice Boltzmann model to handle moving boundaries. The coupling is one way (the fluid does not use any force on the object whereas opposite is the case), anyway the result is visually interesting. Because I had some real-time application in mind I didn't release it in any form but.. after staying for years on my hard disc I decided, yup.. let's give it to people :-)

In short: the method uses interpolated positions of the object and it's velocity to find a force acting on the fluid (or velocity corrections to be strict with the algorihm given in [1]). It wasn't very hard to implement (taking into account I already had the LBM implemented), you must take care about interpolation details to make it consistent.


Here we go with result (click to see videos):

 

More about LBM and other simulations on my website.


Literature
[1] Chang Shu, Ningyu Liu, Yongtian Chew, Zhiliang Lu, Numerical Simulation of Fish Motion by Using Lattice Boltzmann- Immersed Boundary Velocity Correction Method, Journal of Mechanical Science and Technology 21 (2007) 1139~114, link to pdf.

Lifeforce Simulation

IFS fractals

My first try in IFS fractals iterated from simple equations.

They are created from the following equation, You start with any x0,y0 and iterate ti get the pixels to draw x1,y1 and then x2,y2 and so on. Use the formula:
         x = a11 * x0 + a12 * y0 + c1;
        y = a21 * x0 + a22 * y0 + c2;

where a's and c's take from the tables.
The tree:
// a11 a12 a21 a22 c1 c2 s J


double a1[4][8]=    {    {-0.67, -0.02, -0.18, 0.81, 0.00, 1.02, 0.8530, 0.546},
                        { 0.40, 0.40, -0.10, 0.40, -0.04, 0.06, 0.6217, 0.200},
                        {-0.40, -0.40, -0.10, 0.40, 0.04, 0.06, 0.6217, 0.200},
                        {-0.10, 0.00, 0.44, 0.44, 0.00, -0.14, 0.6263, 0.044}
                    };
The leave:
double a1[4][8]=    {    {0.240, -0.007, 0.007, 0.015, 0.000, 0, 0.004},
                        {0.220, -0.330, 0.360, 0.100, 0.540, 0, 0.149},
                        {0.140, -0.360, -0.380, -0.100, 1.400, 0, 0.160},
                        {0.800, 0.100, -0.100, 0.800, 1.600, 0, 0.687}
                    };
The code written here was C++ with PPM output. The number of iterations was
1e11 for leave and 1e9 for tree.
by Maciej Matyka (motivated by K. Graczyk :-)
Literature: J. Kudrewicz, Fraktale i chaos, WNT 2015
Update: see animated version here!.

p.s. The full C++ code that produces these pictures is given below:

---- code starts here

#include <cstdlib>
#include <iostream>
#include <fstream>
#include <ctime>
#include <random>

using namespace std;

const double xmin = -2.25;
const double xmax = 2.25;//15;
const double ymin = 0;
const double ymax = 5;
const double dx = (xmax-xmin);
const double dy = (ymax-ymin);

const int W=int(dx*280);
const int H=int(dy*280);

int pixels[W][H]={0};


//vector<float> X;
//vector<float> Y;

void plot(double x, double y)
{
    int mx = int((x+dx/2)*W/dx);        // -1.25 : 1.25
    int my = int((y)*H/dy);                // 0 : 6
    if(mx<W&&mx>0&&my>0&&my<H)
    {
        pixels[ mx ][ my ] = 1;
    }

}

void saveppm()
{
    ofstream file("pix.ppm");
    int c;
    file << "P3" << endl;
    file << "# dkfjsdljf" << endl;
    file << W << " " << H << endl;
    file << 255 << endl;

    for(int j=0; j<H; j++)
    {
        for(int i=0; i<W; i++)

        {
            c = 255-pixels[i][H-j] * 255;
            file << c << " " << c << " " << c << " ";
        }
        file << endl;
    }

    file.close();
}

Result of the code (in ppm format).

// two stressed pilots
int main()
{
    double x,y;
    double x0,y0;


    mt19937 mt_rand(time(0));


    srand(time(NULL));

    // J. Kudrewicz, p. 24
    // a11 a12 a21 a22 c1 c2 s J
    // choinka
    double a1[4][8]=    {    {-0.67, -0.02, -0.18, 0.81, 0.00, 1.02, 0.8530, 0.546},
                        { 0.40, 0.40, -0.10, 0.40, -0.04, 0.06, 0.6217, 0.200},
                        {-0.40, -0.40, -0.10, 0.40, 0.04, 0.06, 0.6217, 0.200},
                        {-0.10, 0.00, 0.44, 0.44, 0.00, -0.14, 0.6263, 0.044}
                    };
    // paproc
/*    double a1[4][8]=    {    {0.240, -0.007, 0.007, 0.015, 0.000, 0, 0.004},
                        {0.220, -0.330, 0.360, 0.100, 0.540, 0, 0.149},
                        {0.140, -0.360, -0.380, -0.100, 1.400, 0, 0.160},
                        {0.800, 0.100, -0.100, 0.800, 1.600, 0, 0.687}
                    };*/

    // iterate
    x0=0.1;
    y0=0.6;

    long int step=0;
    int i;
    const long int ITERATIONS = 1000000000;

    //cout << mt_rand() << endl;
    //return 0;

    while(step++ < ITERATIONS)
    {
        //i = rand()%4;       // 0-3
        i = mt_rand()%4;
        x = a1[i][0] * x0 + a1[i][1] * y0 + a1[i][4];
        y = a1[i][2] * x0 + a1[i][3] * y0 + a1[i][5];
        plot(x,y);
        x0=x; y0=y;

        if(!(step % (ITERATIONS/50)))
            cout << "*";

    }
    cout << step << endl;
    saveppm();

    return 0;
}
---- code ends here