cloth simulation

cloth_gl.cpp

#include <stdio.h>
#include <math.h>
#include <windows.h>
#include <gl/GL.h>

// simulation
#define SIZE_X 129
#define SIZE_Y 129
#define RIGID_FACTOR   0.5f
#define DAMPING_FACTOR 0.998f
#define TIME_STEP      (1.0f / 60.0f)
#define RELAX_COUNT    3
#define GRAVITY        -9.81f
// sphere
#define SPHERE_RADIUS  16.0f
// camera
#define CAMERA_FOV    (3.14159f * 0.5f)
#define CAMERA_DIST   SIZE_Y
#define CAMERA_Z_NEAR 1.0f
#define CAMERA_Z_FAR  (CAMERA_DIST * 2)
#define CAMERA_Y      (-SIZE_Y * 0.5f)

int width, height;

struct vec3 {
    float x, y, z;

    vec3() {}
    vec3(float x, float y, float z) : x(x), y(y), z(z) {}

    vec3& operator += (const vec3 &v) { x += v.x; y += v.y; z += v.z; return *this; }
    vec3& operator -= (const vec3 &v) { x -= v.x; y -= v.y; z -= v.z; return *this; }
    vec3& operator *= (float s)       { x *= s;   y *= s;   z *= s;   return *this; }

    vec3 operator + (const vec3 &v) const { return vec3(x + v.x, y + v.y, z + v.z); }
    vec3 operator - (const vec3 &v) const { return vec3(x - v.x, y - v.y, z - v.z); }
    vec3 operator * (float s)       const { return vec3(x * s,   y * s,   z * s);   }

    float length2() const { return x * x + y * y + z * z; }
    float length()  const { return sqrtf(length2()); }
    vec3  cross(const vec3 &v) const { return vec3(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x); }
};

// cloth geometry
typedef unsigned int Index;

struct Vertex {
    vec3 normal;
    vec3 pos;
    vec3 old;
    int  links[2];
} *vertices;
Index *indices;

int vCount, iCount;

// camera
vec3 camRot;

// sphere
vec3 sphere;

// input
bool keys[256];
vec3 mousePos, mouseLast, mouseDelta;

void init() {
    vCount = SIZE_X * SIZE_Y;
    vertices = new Vertex[vCount];
    for (int y = 0; y < SIZE_Y; y++)
        for (int x = 0; x < SIZE_X; x++) {
            int i = y * SIZE_X + x;
            Vertex &v = vertices[i];
            v.pos = v.old = vec3(float(x) - SIZE_X * 0.5f, float(-y), float(y) * float(x) / SIZE_X);
            v.links[0] = (x > SIZE_X - 2) ? -1 : i + 1;
            v.links[1] = (y > SIZE_Y - 2) ? -1 : i + SIZE_X;
        }

    iCount = (SIZE_X - 1) * (SIZE_Y - 1) * 2 * 3; // cells_count * triangles_per_cell * indices_per_triangle
    indices = new Index[iCount];
    Index *index = indices;
    for (int y = 0; y < SIZE_Y - 1; y++)
        for (int x = 0; x < SIZE_X - 1; x++) {
            int i = y * SIZE_X + x;

            *index++ = i;
            *index++ = i + 1;
            *index++ = i + SIZE_X + 1;

            *index++ = i;
            *index++ = i + SIZE_X + 1;
            *index++ = i + SIZE_X;
         }

    sphere = vec3(0, -SIZE_Y * 0.5f, 0);

    camRot = vec3(0, 0, 0);

    memset(keys, 0, sizeof(keys));
}

void free() {
    delete[] vertices;
    delete[] indices;
}

inline bool isConstraint(int index) {
    return index < SIZE_X;
}

void integrate() {
    for (int i = 0; i < vCount; i++) {
        if (isConstraint(i)) continue; // skip constraint points
        Vertex &v = vertices[i];
        vec3 delta = v.pos - v.old;
        v.old = v.pos;
        v.pos += delta * DAMPING_FACTOR + vec3(0, GRAVITY * (TIME_STEP * TIME_STEP), 0);
    }
}

void collide() {
    for (int i = 0; i < vCount; i++) {
        if (isConstraint(i)) continue;
        Vertex &v = vertices[i];
    // sphere
        vec3 delta = sphere - v.pos;
        float dist = delta.length2();
        if (dist < (SPHERE_RADIUS * SPHERE_RADIUS)) {
            dist = sqrtf(dist);
            v.pos += delta * (1.0f - SPHERE_RADIUS / dist);
        }
    // floor
        if (v.pos.y < -SIZE_Y) v.pos.y = -SIZE_Y;
    }
}

bool resolve(int aIndex, int bIndex) {
    bool ca = isConstraint(aIndex);
    bool cb = isConstraint(bIndex);

    if (ca && cb) return true;
    
    Vertex &a = vertices[aIndex];
    Vertex &b = vertices[bIndex];

    vec3  delta = b.pos - a.pos;
    float k     = ((1.0f - 1.0f / delta.length()) * RIGID_FACTOR);
    if (k > RIGID_FACTOR * 0.7f) return false;
    vec3  dir   = delta * k;

    if (!ca && !cb) {
        dir *= 0.5f;
        a.pos += dir;
        b.pos -= dir;
    } else if (ca) {
        b.pos -= dir;
    } else if (cb) {
        a.pos += dir;
    }
    return true;
}

void collapse(int index) {
    int x = index % SIZE_X;
    int y = index / SIZE_X;
    if (x > SIZE_X - 1 || y > SIZE_Y - 1) return;

    int i = (y * (SIZE_X - 1) + x) * 2 * 3 - 6;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;

    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
    indices[i++] = 0;
}

void relax() {
    for (int y = 0; y < SIZE_Y; y++)
        for (int x = 0; x < SIZE_X; x++) {
            int i = y * SIZE_X + x;
            Vertex &v = vertices[i];
            if ((v.links[0] > -1 && !resolve(i, v.links[0])) || (v.links[1] > -1 && !resolve(i, v.links[1]))) {
                v.links[0] = v.links[1] = -1;
                collapse(i);
            }
        }
}

void update() {
    integrate();
    for (int i = 0; i < RELAX_COUNT; i++)
        relax();
    collide();
}

void updateCamera(float deltaTime) {
    if (keys[VK_LBUTTON]) {
        if (keys[VK_CONTROL]) {
            sphere += vec3(mouseDelta.x, 0, mouseDelta.y) * 0.1f;
        } else {
            camRot.x -= mouseDelta.y * 0.2f;
            camRot.y -= mouseDelta.x * 0.2f;
        }
        mouseDelta = vec3(0, 0, 0);
    }
}

void updateNormals() {
    for (int y = 0; y < SIZE_Y - 1; y++)
        for (int x = 0; x < SIZE_X - 1; x++) {
            int i = y * SIZE_X + x;
            vec3 a = vertices[i].pos - vertices[i + 1].pos;
            vec3 b = vertices[i].pos - vertices[i + SIZE_X].pos;
            vertices[i].normal = b.cross(a);
        }
}

void render() {
    glViewport(0, 0, width, height);

    glDisable(GL_CULL_FACE);
    glEnable(GL_DEPTH_TEST);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

// apply camera
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    float fy = CAMERA_Z_NEAR * tan(CAMERA_FOV * 0.5f);
    float fx = fy * float(width) / float(height);
    glFrustum(-fx, fx, -fy, fy, CAMERA_Z_NEAR, CAMERA_Z_FAR); // perspective projection
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(0, 0, -CAMERA_DIST);
    glRotatef(-camRot.z, 0, 0, 1);
    glRotatef(-camRot.x, 1, 0, 0);
    glRotatef(-camRot.y, 0, 1, 0);
    glTranslatef(0, -CAMERA_Y, 0);

// apply lighting
    glEnable(GL_NORMALIZE);
    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    GLfloat lightPos[]  = { 0, -SIZE_Y, SIZE_X, 1 };
    GLfloat lightDiff[] = { 1, 1, 1, 1 };
    GLfloat lightSpec[] = { 1, 1, 1, 1 };

    glLightfv(GL_LIGHT0, GL_POSITION, lightPos);
    glLightf(GL_LIGHT0, GL_LINEAR_ATTENUATION, 0.001f);

// render cloth
    updateNormals();
//    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_NORMAL_ARRAY);
    glVertexPointer(3, GL_FLOAT, sizeof(Vertex), &vertices->pos);
    glNormalPointer(GL_FLOAT, sizeof(Vertex), &vertices->normal);
    glDrawElements(GL_TRIANGLES, iCount, GL_UNSIGNED_INT, indices);
    glDisableClientState(GL_VERTEX_ARRAY);
    glDisableClientState(GL_NORMAL_ARRAY);
//    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
}


// =============================== WINDOWS THINGS =============================== 
int getTime() {
    LARGE_INTEGER Freq, Count;
    QueryPerformanceFrequency(&Freq);
    QueryPerformanceCounter(&Count);
    return int(Count.QuadPart * 1000L / Freq.QuadPart);
}

HGLRC initGL(HDC hDC) {
    PIXELFORMATDESCRIPTOR pfd;
    memset(&pfd, 0, sizeof(pfd));
    pfd.nSize = sizeof(pfd);
    pfd.nVersion = 1;
    pfd.dwFlags = PFD_DRAW_TO_WINDOW | PFD_SUPPORT_OPENGL | PFD_DOUBLEBUFFER;
    pfd.cColorBits = 32;
    pfd.cDepthBits = 24;

    int format = ChoosePixelFormat(hDC, &pfd);
    SetPixelFormat(hDC, format, &pfd);
    HGLRC hRC = wglCreateContext(hDC);
    wglMakeCurrent(hDC, hRC);
    return hRC;
}

void freeGL(HGLRC hRC) {
    wglMakeCurrent(0, 0);
    wglDeleteContext(hRC);
}

static LRESULT CALLBACK WndProc(HWND hWnd, UINT msg, WPARAM wParam, LPARAM lParam) {
    switch (msg) {
        case WM_SIZE :
            width = LOWORD(lParam);
            height = HIWORD(lParam);
            break;
        case WM_DESTROY :
            PostQuitMessage(0);
            break;

        case WM_KEYDOWN :
        case WM_KEYUP   :
            keys[wParam] = msg == WM_KEYDOWN;
            break;
        case WM_MOUSEMOVE     :
        case WM_LBUTTONDOWN   :
        case WM_LBUTTONUP     :
        case WM_LBUTTONDBLCLK : {
            mousePos = vec3(float(short(LOWORD(lParam))), float(short(HIWORD(lParam))), 0.0f);

            if (msg != WM_MOUSEMOVE) {
                bool down = msg != WM_LBUTTONUP && msg != WM_RBUTTONUP && msg != WM_MBUTTONUP;
                if (down)
                    SetCapture(hWnd);
                else
                    ReleaseCapture();
                keys[VK_LBUTTON] = down;
                mouseLast  = mousePos;
            }

            mouseDelta = mousePos - mouseLast;
            mouseLast  = mousePos;
            break;
        }
        default:
            return DefWindowProc(hWnd, msg, wParam, lParam);
    }
    return 0;
}

void main() {
    RECT r = {0, 0, 1280, 720};
    AdjustWindowRect(&r, WS_OVERLAPPEDWINDOW, false);
    HWND hWnd = CreateWindow("static", "Cloth Simulation", WS_OVERLAPPEDWINDOW, 0, 0, r.right - r.left, r.bottom - r.top, 0, 0, 0, 0);

    HDC hDC = GetDC(hWnd);
    HGLRC hRC = initGL(hDC);

    SetWindowLong(hWnd, GWL_WNDPROC, (LONG)&WndProc);
    ShowWindow(hWnd, SW_SHOWDEFAULT);

    init();

    float dt = 0.0f;
    int lastTime = getTime();
    MSG msg;

    do {
        if (PeekMessage(&msg, 0, 0, 0, PM_REMOVE)) {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        } else {
            int time = getTime();
            float delta = (time - lastTime) * 0.001f;
            lastTime = time;

            updateCamera(delta);
            dt += delta;
            while (dt >= TIME_STEP) {
                update();
                dt -= TIME_STEP;
            }

            render();
            SwapBuffers(hDC);
        }
    } while (msg.message != WM_QUIT);

    free();

    freeGL(hRC);
    ReleaseDC(hWnd, hDC);
    DestroyWindow(hWnd);
}