initial

20 files changed, 2161 insertions, 0 deletions

diff --git a/wave-sim/src/glwidget.cpp b/wave-sim/src/glwidget.cpp
new file mode 100755
index 0000000..ce8af46
--- /dev/null
+++ b/wave-sim/src/glwidget.cpp
@@ -0,0 +1,194 @@
+#include "glwidget.h"
+
+#include <QApplication>
+#include <QKeyEvent>
+#include <iostream>
+
+#define SPEED 1.5
+#define ROTATE_SPEED 0.0025
+
+using namespace std;
+
+GLWidget::GLWidget(QWidget *parent) :
+    QOpenGLWidget(parent),
+    m_deltaTimeProvider(),
+    m_intervalTimer(),
+    m_sim(),
+    m_camera(),
+    m_shader(),
+    m_forward(),
+    m_sideways(),
+    m_vertical(),
+    m_lastX(),
+    m_lastY(),
+    m_capture(false)
+{
+    // GLWidget needs all mouse move events, not just mouse drag events
+    setMouseTracking(true);
+
+    // Hide the cursor since this is a fullscreen app
+    QApplication::setOverrideCursor(Qt::ArrowCursor);
+
+    // GLWidget needs keyboard focus
+    setFocusPolicy(Qt::StrongFocus);
+
+    // Function tick() will be called once per interva
+    connect(&m_intervalTimer, SIGNAL(timeout()), this, SLOT(tick()));
+}
+
+GLWidget::~GLWidget()
+{
+    if (m_shader != nullptr) delete m_shader;
+}
+
+// ================== Basic OpenGL Overrides
+
+void GLWidget::initializeGL()
+{
+    // Initialize GL extension wrangler
+    glewExperimental = GL_TRUE;
+    GLenum err = glewInit();
+    if (err != GLEW_OK) fprintf(stderr, "Error while initializing GLEW: %s\n", glewGetErrorString(err));
+    fprintf(stdout, "Successfully initialized GLEW %s\n", glewGetString(GLEW_VERSION));
+
+    // Set clear color to white
+    glClearColor(1, 1, 1, 1);
+
+    // Enable depth-testing and backface culling
+    glEnable(GL_DEPTH_TEST);
+    glEnable(GL_CULL_FACE);
+    glCullFace(GL_BACK);
+
+    // Initialize the shader and simulation
+    m_shader = new Shader(":/resources/shaders/shader.vert", ":/resources/shaders/shader.frag");
+    m_sim.init();
+
+    // Initialize camera with a reasonable transform
+    Eigen::Vector3f eye    = {0, 2, -5};
+    Eigen::Vector3f target = {0, 1,  0};
+    m_camera.lookAt(eye, target);
+    m_camera.setOrbitPoint(target);
+    m_camera.setPerspective(120, width() / static_cast<float>(height()), 0.1, 50);
+
+    m_deltaTimeProvider.start();
+//    m_intervalTimer.start(1);
+//    m_intervalTimer.start(1000 / 60);
+    m_intervalTimer.start(1000 / 120);
+}
+
+void GLWidget::paintGL()
+{
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+    m_shader->bind();
+    m_shader->setUniform("proj", m_camera.getProjection());
+    m_shader->setUniform("view", m_camera.getView());
+    m_sim.draw(m_shader);
+    m_shader->unbind();
+}
+
+void GLWidget::resizeGL(int w, int h)
+{
+    glViewport(0, 0, w, h);
+    m_camera.setAspect(static_cast<float>(w) / h);
+}
+
+// ================== Event Listeners
+
+void GLWidget::mousePressEvent(QMouseEvent *event)
+{
+    m_capture = true;
+    m_lastX = event->position().x();
+    m_lastY = event->position().y();
+}
+
+void GLWidget::mouseMoveEvent(QMouseEvent *event)
+{
+    if (!m_capture) return;
+
+    int currX  = event->position().x();
+    int currY  = event->position().y();
+
+    int deltaX = currX - m_lastX;
+    int deltaY = currY - m_lastY;
+
+    if (deltaX == 0 && deltaY == 0) return;
+
+    m_camera.rotate(deltaY * ROTATE_SPEED,
+                    -deltaX * ROTATE_SPEED);
+
+    m_lastX = currX;
+    m_lastY = currY;
+}
+
+void GLWidget::mouseReleaseEvent(QMouseEvent *event)
+{
+    m_capture = false;
+}
+
+void GLWidget::wheelEvent(QWheelEvent *event)
+{
+    float zoom = 1 - event->pixelDelta().y() * 0.1f / 120.f;
+    m_camera.zoom(zoom);
+}
+
+void GLWidget::keyPressEvent(QKeyEvent *event)
+{
+    if (event->isAutoRepeat()) return;
+
+    switch (event->key())
+    {
+    case Qt::Key_W: m_forward  += SPEED; break;
+    case Qt::Key_S: m_forward  -= SPEED; break;
+    case Qt::Key_A: m_sideways -= SPEED; break;
+    case Qt::Key_D: m_sideways += SPEED; break;
+    case Qt::Key_F: m_vertical -= SPEED; break;
+    case Qt::Key_R: m_vertical += SPEED; break;
+    case Qt::Key_C: m_camera.toggleIsOrbiting(); break;
+    case Qt::Key_T: m_sim.toggleWire(); break;
+    case Qt::Key_O: m_sim.toggleForceRender(); break;
+    case Qt::Key_Escape: QApplication::quit();
+    }
+}
+
+void GLWidget::keyReleaseEvent(QKeyEvent *event)
+{
+    if (event->isAutoRepeat()) return;
+
+    switch (event->key())
+    {
+    case Qt::Key_W: m_forward  -= SPEED; break;
+    case Qt::Key_S: m_forward  += SPEED; break;
+    case Qt::Key_A: m_sideways += SPEED; break;
+    case Qt::Key_D: m_sideways -= SPEED; break;
+    case Qt::Key_F: m_vertical += SPEED; break;
+    case Qt::Key_R: m_vertical -= SPEED; break;
+    }
+}
+
+// ================== Physics Tick
+
+
+void GLWidget::tick()
+{
+    float deltaSeconds = m_deltaTimeProvider.restart() / 1000.f;
+//    deltaSeconds = 0.5f;
+//    std::cout << deltaSeconds << std::endl;
+    // .02 is optimal, .021 disappears after 4 for sphere
+    // .015 okay for ellipsoid
+    // .001 okay for cone
+    deltaSeconds = .01;
+    deltaSeconds = .014;
+    m_sim.update(deltaSeconds);
+
+    // Move camera
+    auto look = m_camera.getLook();
+    look.y() = 0;
+    look.normalize();
+    Eigen::Vector3f perp(-look.z(), 0, look.x());
+    Eigen::Vector3f moveVec = m_forward * look.normalized() + m_sideways * perp.normalized() + m_vertical * Eigen::Vector3f::UnitY();
+    moveVec *= deltaSeconds;
+    m_camera.move(moveVec);
+
+    // Flag this view for repainting (Qt will call paintGL() soon after)
+    update();
+}
diff --git a/wave-sim/src/glwidget.h b/wave-sim/src/glwidget.h
new file mode 100755
index 0000000..d108ccc
--- /dev/null
+++ b/wave-sim/src/glwidget.h
@@ -0,0 +1,62 @@
+#pragma once
+
+#ifdef __APPLE__
+#define GL_SILENCE_DEPRECATION
+#endif
+
+#include "simulation.h"
+#include "graphics/camera.h"
+#include "graphics/shader.h"
+
+#include <QOpenGLWidget>
+#include <QElapsedTimer>
+#include <QTimer>
+#include <memory>
+
+class GLWidget : public QOpenGLWidget
+{
+    Q_OBJECT
+
+public:
+    GLWidget(QWidget *parent = nullptr);
+    ~GLWidget();
+
+private:
+    static const int FRAMES_TO_AVERAGE = 30;
+
+private:
+    // Basic OpenGL Overrides
+    void initializeGL()         override;
+    void paintGL()              override;
+    void resizeGL(int w, int h) override;
+
+    // Event Listeners
+    void mousePressEvent  (QMouseEvent *event) override;
+    void mouseMoveEvent   (QMouseEvent *event) override;
+    void mouseReleaseEvent(QMouseEvent *event) override;
+    void wheelEvent       (QWheelEvent *event) override;
+    void keyPressEvent    (QKeyEvent   *event) override;
+    void keyReleaseEvent  (QKeyEvent   *event) override;
+
+private:
+    QElapsedTimer m_deltaTimeProvider; // For measuring elapsed time
+    QTimer        m_intervalTimer;     // For triggering timed events
+
+    Simulation m_sim;
+    Camera     m_camera;
+    Shader    *m_shader;
+
+    int m_forward;
+    int m_sideways;
+    int m_vertical;
+
+    int m_lastX;
+    int m_lastY;
+
+    bool m_capture;
+
+private slots:
+
+    // Physics Tick
+    void tick();
+};
diff --git a/wave-sim/src/graphics/camera.cpp b/wave-sim/src/graphics/camera.cpp
new file mode 100644
index 0000000..85fc7d9
--- /dev/null
+++ b/wave-sim/src/graphics/camera.cpp
@@ -0,0 +1,188 @@
+#include "graphics/camera.h"
+
+#include <iostream>
+
+Camera::Camera()
+    : m_position(0,0,0),
+      m_pitch(0), m_yaw(0),
+      m_look(0, 0, 1),
+      m_orbitPoint(0, 0, 0),
+      m_isOrbiting(false),
+      m_view(Eigen::Matrix4f::Identity()),
+      m_proj(Eigen::Matrix4f::Identity()),
+      m_viewDirty(true),
+      m_projDirty(true),
+      m_fovY(90), m_aspect(1), m_near(0.1f), m_far(50.f),
+      m_zoom(1)
+{}
+
+// ================== Position
+
+void Camera::setPosition(const Eigen::Vector3f &position)
+{
+    m_position  = position;
+    m_viewDirty = true;
+}
+
+void Camera::move(const Eigen::Vector3f &deltaPosition)
+{
+    if (deltaPosition.squaredNorm() == 0) return;
+
+    m_position += deltaPosition;
+
+    if (m_isOrbiting) {
+        m_orbitPoint += deltaPosition;
+    }
+
+    m_viewDirty = true;
+
+}
+
+// ================== Rotation
+
+void Camera::setRotation(float pitch, float yaw)
+{
+    m_pitch = pitch;
+    m_yaw   = yaw;
+    m_viewDirty = true;
+    updateLook();
+}
+
+void Camera::rotate(float deltaPitch, float deltaYaw)
+{
+    m_pitch += deltaPitch;
+    m_yaw   += deltaYaw;
+    m_pitch  = std::clamp(m_pitch, (float) -M_PI_2 + 0.01f, (float) M_PI_2 - 0.01f);
+    m_viewDirty = true;
+    updateLook();
+
+    if (m_isOrbiting) {
+        m_position = m_orbitPoint - m_look * m_zoom;
+    }
+}
+
+// ================== Position and Rotation
+
+void Camera::lookAt(const Eigen::Vector3f &eye, const Eigen::Vector3f &target)
+{
+    m_position  = eye;
+    m_look      = (target - eye).normalized();
+    m_viewDirty = true;
+    updatePitchAndYaw();
+}
+
+// ================== Orbiting
+
+void Camera::setOrbitPoint(const Eigen::Vector3f &orbitPoint)
+{
+    m_orbitPoint = orbitPoint;
+    m_viewDirty = true;
+}
+
+bool Camera::getIsOrbiting()
+{
+    return m_isOrbiting;
+}
+
+void Camera::setIsOrbiting(bool isOrbiting)
+{
+    m_isOrbiting = isOrbiting;
+    m_viewDirty = true;
+}
+
+void Camera::toggleIsOrbiting()
+{
+    m_isOrbiting = !m_isOrbiting;
+    m_viewDirty = true;
+
+    if (m_isOrbiting) {
+        m_zoom = (m_orbitPoint - m_position).norm();
+        m_look = (m_orbitPoint - m_position).normalized();
+        updatePitchAndYaw();
+    }
+}
+
+void Camera::zoom(float zoomMultiplier)
+{
+    if (!m_isOrbiting) return;
+
+    m_zoom *= zoomMultiplier;
+    m_position = m_orbitPoint - m_look * m_zoom;
+    m_viewDirty = true;
+}
+
+// ================== Intrinsics
+
+void Camera::setPerspective(float fovY, float aspect, float near, float far)
+{
+    m_fovY = fovY;
+    m_aspect = aspect;
+    m_near = near;
+    m_far = far;
+    m_projDirty = true;
+}
+
+void Camera::setAspect(float aspect)
+{
+    m_aspect = aspect;
+    m_projDirty = true;
+}
+
+// ================== Important Getters
+
+const Eigen::Matrix4f &Camera::getView()
+{
+    if (m_viewDirty) {
+        Eigen::Matrix3f R;
+        Eigen::Vector3f f = m_look.normalized();
+        Eigen::Vector3f u = Eigen::Vector3f::UnitY();
+        Eigen::Vector3f s = f.cross(u).normalized();
+        u = s.cross(f);
+        R.col(0) = s;
+        R.col(1) = u;
+        R.col(2) = -f;
+        m_view.topLeftCorner<3, 3>() = R.transpose();
+        m_view.topRightCorner<3, 1>() = -R.transpose() * m_position;
+        m_view(3, 3) = 1.f;
+        m_viewDirty = false;
+    }
+    return m_view;
+}
+
+const Eigen::Matrix4f &Camera::getProjection()
+{
+    if(m_projDirty) {
+        float theta = m_fovY * 0.5f;
+        float invRange = 1.f / (m_far - m_near);
+        float invtan = 1.f / tanf(theta);
+        m_proj(0, 0) = invtan / m_aspect;
+        m_proj(1, 1) = invtan;
+        m_proj(2, 2) = -(m_near + m_far) * invRange;
+        m_proj(3, 2) = -1;
+        m_proj(2, 3) = -2 * m_near * m_far * invRange;
+        m_proj(3, 3) = 0;
+        m_projDirty = false;
+    }
+    return m_proj;
+}
+
+const Eigen::Vector3f &Camera::getLook()
+{
+    return m_look;
+}
+
+// ================== Private Helpers
+
+void Camera::updateLook()
+{
+    m_look = Eigen::Vector3f(0, 0, 1);
+    m_look = Eigen::AngleAxis<float>(m_pitch, Eigen::Vector3f::UnitX()) * m_look;
+    m_look = Eigen::AngleAxis<float>(m_yaw,   Eigen::Vector3f::UnitY()) * m_look;
+    m_look = m_look.normalized();
+}
+
+void Camera::updatePitchAndYaw()
+{
+    m_pitch = asinf(-m_look.y());
+    m_yaw   = atan2f(m_look.x(), m_look.z());
+}
diff --git a/wave-sim/src/graphics/camera.h b/wave-sim/src/graphics/camera.h
new file mode 100644
index 0000000..04586af
--- /dev/null
+++ b/wave-sim/src/graphics/camera.h
@@ -0,0 +1,51 @@
+#pragma once
+
+#include "Eigen/Dense"
+
+class Camera
+{
+public:
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+    Camera();
+
+    void setPosition(const Eigen::Vector3f      &position);
+    void move       (const Eigen::Vector3f &deltaPosition);
+
+    void setRotation(float      pitch, float      yaw);
+    void rotate     (float deltaPitch, float deltaYaw);
+
+    void lookAt(const Eigen::Vector3f &eye, const Eigen::Vector3f &target);
+
+    void setOrbitPoint(const Eigen::Vector3f &target);
+    bool getIsOrbiting();
+    void setIsOrbiting(bool orbit);
+    void toggleIsOrbiting();
+    void zoom(float zoomMultiplier);
+
+    const Eigen::Matrix4f &getView();
+    const Eigen::Matrix4f &getProjection();
+    const Eigen::Vector3f &getLook();
+
+    void setPerspective(float fovY, float aspect, float near, float far);
+    void setAspect(float aspect);
+
+private:
+    void updateLook();
+    void updatePitchAndYaw();
+
+    // Do not mess with the order of these variables. Some Eigen voodoo will cause an inexplicable crash.
+
+    Eigen::Vector3f m_position;
+
+    float m_pitch, m_yaw;
+    Eigen::Vector3f m_look;
+
+    Eigen::Vector3f m_orbitPoint;
+    bool m_isOrbiting;
+
+    Eigen::Matrix4f m_view, m_proj;
+    bool m_viewDirty, m_projDirty;
+
+    float m_fovY, m_aspect, m_near, m_far;
+    float m_zoom;
+};
diff --git a/wave-sim/src/graphics/graphicsdebug.cpp b/wave-sim/src/graphics/graphicsdebug.cpp
new file mode 100644
index 0000000..b9d831c
--- /dev/null
+++ b/wave-sim/src/graphics/graphicsdebug.cpp
@@ -0,0 +1,126 @@
+#include <GL/glew.h>
+#include "graphics/graphicsdebug.h"
+
+#include <iostream>
+#include <vector>
+
+
+void checkError(std::string prefix) {
+    GLenum err = glGetError();
+    if (err != GL_NO_ERROR) {
+        std::cerr << prefix << (prefix == std::string("") ? "" : ": ") << "GL is in an error state before painting." << std::endl;
+        printGLErrorCodeInEnglish(err);
+    }
+}
+
+void printGLErrorCodeInEnglish(GLenum err) {
+    std::cerr << "GL error code " << err << ":" << std::endl;
+    switch(err) {
+        case GL_INVALID_ENUM:
+            std::cerr << "GL_INVALID_ENUM" << std::endl;
+            std::cerr << "An unacceptable value is specified for an enumerated argument. The offending command is ignored and has no other side effect than to set the error flag." << std::endl;
+            break;
+        case GL_INVALID_VALUE:
+            std::cerr << "GL_INVALID_VALUE" << std::endl;
+            std::cerr << "A numeric argument is out of range. The offending command is ignored and has no other side effect than to set the error flag." << std::endl;
+            break;
+        case GL_INVALID_OPERATION:
+            std::cerr << "GL_INVALID_OPERATION" << std::endl;
+            std::cerr << "The specified operation is not allowed in the current state. The offending command is ignored and has no other side effect than to set the error flag." << std::endl;
+            break;
+        case GL_INVALID_FRAMEBUFFER_OPERATION:
+            std::cerr << "GL_INVALID_FRAMEBUFFER_OPERATION" << std::endl;
+            std::cerr << "The framebuffer object is not complete. The offending command is ignored and has no other side effect than to set the error flag." << std::endl;
+            break;
+        case GL_OUT_OF_MEMORY:
+            std::cerr << "GL_OUT_OF_MEMORY" << std::endl;
+            std::cerr << "There is not enough memory left to execute the command. The state of the GL is undefined, except for the state of the error flags, after this error is recorded." << std::endl;
+            break;
+        case GL_STACK_UNDERFLOW:
+            std::cerr << "GL_STACK_UNDERFLOW" << std::endl;
+            std::cerr << "An attempt has been made to perform an operation that would cause an internal stack to underflow." << std::endl;
+            break;
+        case GL_STACK_OVERFLOW:
+            std::cerr << "GL_STACK_OVERFLOW" << std::endl;
+            std::cerr << "An attempt has been made to perform an operation that would cause an internal stack to overflow." << std::endl;
+            break;
+        default:
+            std::cerr << "Unknown GL error code" << std::endl;
+    }
+}
+
+void checkFramebufferStatus() {
+    GLenum status = glCheckFramebufferStatus(GL_FRAMEBUFFER);
+    if (status != GL_FRAMEBUFFER_COMPLETE) {
+        std::cerr << "Framebuffer is incomplete." << std::endl;
+        printFramebufferErrorCodeInEnglish(status);
+    }
+}
+
+void printFramebufferErrorCodeInEnglish(GLenum err) {
+    switch(err) {
+        case GL_FRAMEBUFFER_UNDEFINED:
+            std:: cerr << "GL_FRAMEBUFFER_UNDEFINED is returned if the specified framebuffer is the default read or draw framebuffer, but the default framebuffer does not exist." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_ATTACHMENT is returned if any of the framebuffer attachment points are framebuffer incomplete." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_MISSING_ATTACHMENT is returned if the framebuffer does not have at least one image attached to it." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_DRAW_BUFFER:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_DRAW_BUFFER is returned if the value of GL_FRAMEBUFFER_ATTACHMENT_OBJECT_TYPE is GL_NONE for any color attachment point(s) named by GL_DRAW_BUFFERi." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_READ_BUFFER:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_READ_BUFFER is returned if GL_READ_BUFFER is not GL_NONE and the value of GL_FRAMEBUFFER_ATTACHMENT_OBJECT_TYPE is GL_NONE for the color attachment point named by GL_READ_BUFFER." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_UNSUPPORTED:
+            std::cerr << "GL_FRAMEBUFFER_UNSUPPORTED is returned if the combination of internal formats of the attached images violates an implementation-dependent set of restrictions." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE is returned if the value of GL_RENDERBUFFER_SAMPLES is not the same for all attached renderbuffers; if the value of GL_TEXTURE_SAMPLES is the not same for all attached textures; or, if the attached images are a mix of renderbuffers and textures, the value of GL_RENDERBUFFER_SAMPLES does not match the value of GL_TEXTURE_SAMPLES." << std::endl;
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_MULTISAMPLE is also returned if the value of GL_TEXTURE_FIXED_SAMPLE_LOCATIONS is not the same for all attached textures; or, if the attached images are a mix of renderbuffers and textures, the value of GL_TEXTURE_FIXED_SAMPLE_LOCATIONS is not GL_TRUE for all attached textures." << std::endl;
+            break;
+        case GL_FRAMEBUFFER_INCOMPLETE_LAYER_TARGETS:
+            std::cerr << "GL_FRAMEBUFFER_INCOMPLETE_LAYER_TARGETS is returned if any framebuffer attachment is layered, and any populated attachment is not layered, or if all populated color attachments are not from textures of the same target." << std::endl;
+            break;
+    }
+}
+
+void checkShaderCompilationStatus(GLuint shaderID) {
+    GLint status;
+    glGetShaderiv(shaderID, GL_COMPILE_STATUS, &status);
+    if (status == GL_FALSE) {
+        std::cerr << "Error: Could not compile shader."  << std::endl;
+
+        GLint maxLength = 0;
+        glGetShaderiv(shaderID, GL_INFO_LOG_LENGTH, &maxLength);
+
+        // The maxLength includes the null character
+        std::vector<GLchar> errorLog(maxLength);
+        glGetShaderInfoLog(shaderID, maxLength, &maxLength, &errorLog[0]);
+
+        std::cerr << &errorLog[0] << std::endl;
+    } else {
+        std::cerr << "Shader compiled." << std::endl;
+    }
+}
+
+void checkShaderLinkStatus(GLuint shaderProgramID) {
+    GLint linked;
+    glGetProgramiv(shaderProgramID, GL_LINK_STATUS, &linked);
+    if (linked == GL_FALSE) {
+        std::cerr << "Shader failed to link" << std::endl;
+
+        GLint maxLength = 0;
+        glGetProgramiv(shaderProgramID, GL_INFO_LOG_LENGTH, &maxLength);
+
+        // The maxLength includes the null character
+        std::vector<GLchar> errorLog(maxLength);
+        glGetProgramInfoLog(shaderProgramID, maxLength, &maxLength, &errorLog[0]);
+
+        std::cerr << &errorLog[0] << std::endl;
+    } else {
+        std::cerr << "Shader linked successfully." << std::endl;
+    }
+}
diff --git a/wave-sim/src/graphics/graphicsdebug.h b/wave-sim/src/graphics/graphicsdebug.h
new file mode 100644
index 0000000..9be33b4
--- /dev/null
+++ b/wave-sim/src/graphics/graphicsdebug.h
@@ -0,0 +1,15 @@
+#pragma once
+
+#include <GL/glew.h>
+#include <string>
+
+#define GRAPHICS_DEBUG_LEVEL 0
+
+void checkError(std::string prefix = "");
+void printGLErrorCodeInEnglish(GLenum err);
+
+void checkFramebufferStatus();
+void printFramebufferErrorCodeInEnglish(GLenum err);
+
+void checkShaderCompilationStatus(GLuint shaderID);
+void checkShaderLinkStatus(GLuint shaderProgramID);
diff --git a/wave-sim/src/graphics/meshloader.cpp b/wave-sim/src/graphics/meshloader.cpp
new file mode 100644
index 0000000..84e78c7
--- /dev/null
+++ b/wave-sim/src/graphics/meshloader.cpp
@@ -0,0 +1,53 @@
+#include "graphics/meshloader.h"
+
+#define TINYOBJLOADER_IMPLEMENTATION
+#include "util/tiny_obj_loader.h"
+
+#include <iostream>
+
+#include <QString>
+#include <QFile>
+#include <QTextStream>
+#include <QRegularExpression>
+
+using namespace Eigen;
+
+bool MeshLoader::loadTetMesh(const std::string &filepath, std::vector<Eigen::Vector3d> &vertices, std::vector<Eigen::Vector4i> &tets)
+{
+    QString qpath = QString::fromStdString(filepath);
+    QFile file(qpath);
+
+    if(!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
+        std::cout << "Error opening file: " << filepath << std::endl;
+        return false;
+    }
+    QTextStream in(&file);
+
+    QRegularExpression vrxp("v (-?\\d*\\.?\\d+) +(-?\\d*\\.?\\d+) +(-?\\d*\\.?\\d+)");
+    QRegularExpression trxp("t (\\d+) +(\\d+) +(\\d+) +(\\d+)");
+
+    while(!in.atEnd()) {
+        QString line = in.readLine();
+        auto match = vrxp.match(line);
+        if(match.hasMatch()) {
+            vertices.emplace_back(match.captured(1).toDouble(),
+                                  match.captured(2).toDouble(),
+                                  match.captured(3).toDouble());
+            continue;
+        }
+        match = trxp.match(line);
+        if(match.hasMatch()) {
+            tets.emplace_back(match.captured(1).toInt(),
+                              match.captured(2).toInt(),
+                              match.captured(3).toInt(),
+                              match.captured(4).toInt());
+        }
+    }
+    file.close();
+    return true;
+}
+
+MeshLoader::MeshLoader()
+{
+
+}
diff --git a/wave-sim/src/graphics/meshloader.h b/wave-sim/src/graphics/meshloader.h
new file mode 100644
index 0000000..e6b87fd
--- /dev/null
+++ b/wave-sim/src/graphics/meshloader.h
@@ -0,0 +1,15 @@
+#pragma once
+
+#include <vector>
+#include "Eigen/Dense"
+#include "Eigen/StdVector"
+
+EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(Eigen::Matrix4i)
+
+class MeshLoader
+{
+public:
+    static bool loadTetMesh(const std::string &filepath, std::vector<Eigen::Vector3d> &vertices, std::vector<Eigen::Vector4i> &tets);
+private:
+    MeshLoader();
+};
diff --git a/wave-sim/src/graphics/shader.cpp b/wave-sim/src/graphics/shader.cpp
new file mode 100644
index 0000000..3789be0
--- /dev/null
+++ b/wave-sim/src/graphics/shader.cpp
@@ -0,0 +1,234 @@
+#include "graphics/shader.h"
+#include "graphics/graphicsdebug.h"
+
+#include <QFile>
+#include <QString>
+#include <QTextStream>
+#include <iostream>
+
+Shader::Shader(const std::string &vertexPath,
+               const std::string &fragmentPath)
+{
+    createProgramID();
+
+    std::vector<GLuint> shaders = {
+        createShaderFromString(getFileContents(vertexPath),   GL_VERTEX_SHADER),
+        createShaderFromString(getFileContents(fragmentPath), GL_FRAGMENT_SHADER)
+    };
+
+    buildShaderProgramFromShaders(shaders);
+    discoverShaderData();
+}
+
+Shader::~Shader()
+{
+    glDeleteProgram(m_programID);
+}
+
+Shader::Shader(Shader &&that)
+    : m_programID(that.m_programID),
+      m_attributes(std::move(that.m_attributes)),
+      m_uniforms(std::move(that.m_uniforms))
+{
+    that.m_programID = 0;
+}
+
+Shader& Shader::operator=(Shader &&that)
+{
+    this->~Shader();
+
+    m_programID = that.m_programID;
+    m_attributes = std::move(that.m_attributes);
+    m_uniforms = std::move(that.m_uniforms);
+
+    that.m_programID = 0;
+
+    return *this;
+}
+
+// ================== Regular Use
+
+void Shader::bind()   const { glUseProgram(m_programID); }
+
+void Shader::unbind() const { glUseProgram(0); }
+
+GLuint Shader::getUniformLocation(std::string name)
+{
+    return glGetUniformLocation(m_programID, name.c_str());
+}
+
+GLuint Shader::getEnumeratedUniformLocation(std::string name, int index)
+{
+    std::string n = name + "[" + std::to_string(index) + "]";
+    return glGetUniformLocation(m_programID, n.c_str());
+}
+
+// ================== Setting Uniforms
+
+// Note: the overload to set matrix uniforms is in the .h file
+
+void Shader::setUniform(const std::string &name, float f)
+{
+    glUniform1f(m_uniforms[name], f);
+}
+
+void Shader::setUniform(const std::string &name, int i)
+{
+    glUniform1i(m_uniforms[name], i);
+}
+
+void Shader::setUniform(const std::string &name, bool b)
+{
+    glUniform1i(m_uniforms[name], static_cast<GLint>(b));
+}
+
+// ================== Creating the Program
+
+void Shader::createProgramID()
+{
+    m_programID = glCreateProgram();
+}
+
+void Shader::buildShaderProgramFromShaders(const std::vector<GLuint> &shaderHandles)
+{
+    // Attach shaders
+    for (const GLuint &shaderHandle : shaderHandles) {
+        glAttachShader(m_programID, shaderHandle);
+    }
+
+    // Link program
+    glLinkProgram(m_programID);
+    checkShaderLinkStatus(m_programID);
+
+    // Detach and delete shaders
+    for (const GLuint &shaderHandle : shaderHandles) {
+        glDetachShader(m_programID, shaderHandle);
+        glDeleteShader(shaderHandle);
+    }
+}
+
+// ================== Creating Shaders From Filepaths
+
+std::string Shader::getFileContents(std::string filepath)
+{
+    QString filepathStr = QString::fromStdString(filepath);
+    QFile file(filepathStr);
+
+    if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
+        throw std::runtime_error(std::string("Failed to open shader: ") + filepath);
+    }
+
+    QTextStream stream(&file);
+    QString contents = stream.readAll();
+    file.close();
+
+    return contents.toStdString();
+}
+
+GLuint Shader::createShaderFromString(const std::string &str, GLenum shaderType)
+{
+    GLuint shaderHandle = glCreateShader(shaderType);
+
+    // Compile shader code
+    const char *codePtr = str.c_str();
+    glShaderSource(shaderHandle, 1, &codePtr, nullptr); // Assumes code is null terminated
+    glCompileShader(shaderHandle);
+
+    checkShaderCompilationStatus(shaderHandle);
+
+    return shaderHandle;
+}
+
+// ================== Discovering Attributes/Uniforms/Textures
+
+void Shader::discoverShaderData() {
+    discoverAttributes();
+    discoverUniforms();
+}
+
+void Shader::discoverAttributes() {
+    bind();
+    GLint attribCount;
+    glGetProgramiv(m_programID, GL_ACTIVE_ATTRIBUTES, &attribCount);
+    for (int i = 0; i < attribCount; i++) {
+        const GLsizei bufSize = 256;
+        GLsizei nameLength = 0;
+        GLint arraySize = 0;
+        GLenum type;
+        GLchar name[bufSize];
+        glGetActiveAttrib(m_programID, i, bufSize, &nameLength, &arraySize, &type, name);
+        name[std::min(nameLength, bufSize - 1)] = 0;
+        m_attributes[std::string(name)] = glGetAttribLocation(m_programID, name);
+    }
+    unbind();
+}
+
+void Shader::discoverUniforms() {
+    bind();
+    GLint uniformCount;
+    glGetProgramiv(m_programID, GL_ACTIVE_UNIFORMS, &uniformCount);
+    for (int i = 0; i < uniformCount; i++) {
+        const GLsizei bufSize = 256;
+        GLsizei nameLength = 0;
+        GLint arraySize = 0;
+        GLenum type;
+        GLchar name[bufSize];
+        glGetActiveUniform(m_programID, i, bufSize, &nameLength, &arraySize, &type, name);
+        name[std::min(nameLength, bufSize - 1)] = 0;
+
+        std::string strname(name);
+        if (isUniformArray(name, nameLength)) {
+            addUniformArray(strname, arraySize);
+        } else if (isTexture(type)) {
+            addTexture(strname);
+        } else {
+            addUniform(strname);
+        }
+    }
+    unbind();
+}
+
+bool Shader::isUniformArray(const GLchar *name, GLsizei nameLength) {
+    // Check if the last 3 characters are '[0]'
+    return (name[nameLength - 3] == '[') &&
+            (name[nameLength - 2] == '0') &&
+            (name[nameLength - 1] == ']');
+}
+
+bool Shader::isTexture(GLenum type) {
+    return (type == GL_SAMPLER_2D) ||
+            (type == GL_SAMPLER_CUBE);
+}
+
+void Shader::addUniformArray(const std::string &name, size_t size) {
+    std::string cleanName = name.substr(0, name.length() - 3);
+    for (auto i = static_cast<size_t>(0); i < size; i++) {
+        std::string enumeratedName = name;
+        enumeratedName[enumeratedName.length() - 2] = static_cast<char>('0' + i);
+        std::tuple< std::string, size_t > nameIndexTuple = std::make_tuple(cleanName, i);
+        m_uniformArrays[nameIndexTuple] = glGetUniformLocation(m_programID, enumeratedName.c_str());
+    }
+
+#if GRAPHICS_DEBUG_LEVEL > 0
+    m_trackedUniformArrays[name] = false;
+#endif
+}
+
+void Shader::addTexture(const std::string &name) {
+    GLint location = glGetUniformLocation(m_programID, name.c_str());
+    m_textureLocations[name] = location;
+    GLint slot = m_textureSlots.size();
+    m_textureSlots[location] = slot; // Assign slots in increasing order.
+
+#if GRAPHICS_DEBUG_LEVEL > 0
+    m_trackedTextures[name] = false;
+#endif
+}
+
+void Shader::addUniform(const std::string &name) {
+    m_uniforms[name] = glGetUniformLocation(m_programID, name.c_str());
+
+#if GRAPHICS_DEBUG_LEVEL > 0
+    m_trackedUniforms[name] = false;
+#endif
+}
diff --git a/wave-sim/src/graphics/shader.h b/wave-sim/src/graphics/shader.h
new file mode 100644
index 0000000..08f6654
--- /dev/null
+++ b/wave-sim/src/graphics/shader.h
@@ -0,0 +1,71 @@
+#pragma once
+
+#include <map>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include <GL/glew.h>
+#include "Eigen/Dense"
+#include "util/unsupportedeigenthing/OpenGLSupport"
+
+
+class Shader {
+public:
+    Shader(const std::string &vertexPath,
+           const std::string &fragmentPath);
+
+    virtual ~Shader();
+
+    Shader(Shader &that)            = delete;
+    Shader& operator=(Shader &that) = delete;
+    Shader(Shader &&that);
+    Shader& operator=(Shader &&that);
+
+    // Basic Usage
+    GLuint getProgramID() const { return m_programID; }
+    void   bind()   const;
+    void   unbind() const;
+    GLuint getUniformLocation(std::string name);
+    GLuint getEnumeratedUniformLocation(std::string name, int index);
+
+    // Setting Uniforms
+    void setUniform(const std::string &name, float f);
+    void setUniform(const std::string &name, int   i);
+    void setUniform(const std::string &name, bool  b);
+    template<typename type, int n, int m>
+    void setUniform(const std::string &name, const Eigen::Matrix<type, n, m> &mat)
+    {
+        glUniform(m_uniforms[name], mat);
+    }
+
+
+private:
+    // Creating the Program
+    void createProgramID();
+    void buildShaderProgramFromShaders(const std::vector<GLuint> &shaders);
+
+    // Creating Shaders From Filepaths
+    std::string getFileContents(std::string path);
+    GLuint createShaderFromString(const std::string &str, GLenum shaderType);
+
+    // Discovering attributes/uniforms/textures
+    void discoverShaderData();
+    void discoverAttributes();
+    void discoverUniforms();
+    bool isUniformArray(const GLchar *name , GLsizei nameLength);
+    bool isTexture(GLenum type);
+    void addUniform(const std::string &name);
+    void addUniformArray(const std::string &name, size_t size);
+    void addTexture(const std::string &name);
+
+    // Identifies the shader program associated with this shader
+    GLuint m_programID;
+
+    // Collections of known attributes/uniforms/textures
+    std::map<std::string, GLuint>                     m_attributes;
+    std::map<std::string, GLuint>                     m_uniforms;
+    std::map<std::tuple<std::string, size_t>, GLuint> m_uniformArrays;
+    std::map<std::string, GLuint>                     m_textureLocations; // name to uniform location
+    std::map<GLuint, GLuint>                          m_textureSlots;     // uniform location to texture slot
+};
diff --git a/wave-sim/src/graphics/shape.cpp b/wave-sim/src/graphics/shape.cpp
new file mode 100644
index 0000000..e59e009
--- /dev/null
+++ b/wave-sim/src/graphics/shape.cpp
@@ -0,0 +1,337 @@
+#include "shape.h"
+
+#include <iostream>
+
+#include "graphics/shader.h"
+
+using namespace Eigen;
+
+Shape::Shape()
+    : m_tetVao(-1),
+      m_numSurfaceVertices(),
+      m_verticesSize(),
+      m_modelMatrix(Eigen::Matrix4f::Identity()),
+      m_wireframe(false)
+{
+}
+
+//void Shape::init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &normals, const std::vector<Eigen::Vector3i> &triangles)
+//{
+//    if(vertices.size() != normals.size()) {
+//        std::cerr << "Vertices and normals are not the same size" << std::endl;
+//        return;
+//    }
+//    glGenBuffers(1, &m_surfaceVbo);
+//    glGenBuffers(1, &m_surfaceIbo);
+//    glGenVertexArrays(1, &m_surfaceVao);
+
+//    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+//    glBufferData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3 * 2, nullptr, GL_DYNAMIC_DRAW);
+//    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * vertices.size() * 3, static_cast<const void *>(vertices.data()));
+//    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3, sizeof(double) * vertices.size() * 3, static_cast<const void *>(normals.data()));
+//    glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+//    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_surfaceIbo);
+//    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int) * 3 * triangles.size(), static_cast<const void *>(triangles.data()), GL_STATIC_DRAW);
+//    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+//    glBindVertexArray(m_surfaceVao);
+//    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+//    glEnableVertexAttribArray(0);
+//    glVertexAttribPointer(0, 3, GL_DOUBLE, GL_FALSE, 0, static_cast<GLvoid *>(0));
+//    glEnableVertexAttribArray(1);
+//    glVertexAttribPointer(1, 3, GL_DOUBLE, GL_FALSE, 0, reinterpret_cast<GLvoid *>(sizeof(double) * vertices.size() * 3));
+//    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_surfaceIbo);
+//    glBindVertexArray(0);
+//    glBindBuffer(GL_ARRAY_BUFFER, 0);
+//    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+//    m_numSurfaceVertices = triangles.size() * 3;
+//    m_verticesSize = vertices.size();
+//    m_faces = triangles;
+//}
+
+void Shape::init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3i> &triangles)
+{
+    std::vector<Eigen::Vector3d> verts;
+    std::vector<Eigen::Vector3d> normals;
+    std::vector<Eigen::Vector3i> faces;
+    std::vector<Eigen::Vector3d> forces;
+    verts.reserve(triangles.size() * 3);
+    normals.reserve(triangles.size() * 3);
+    for(auto& f : triangles) {
+        auto& v1 = vertices[f[0]];
+        auto& v2 = vertices[f[1]];
+        auto& v3 = vertices[f[2]];
+        auto& e1 = v2 - v1;
+        auto& e2 = v3 - v1;
+        auto n = e1.cross(e2);
+        int s = verts.size();
+        faces.push_back(Eigen::Vector3i(s, s + 1, s + 2));
+        normals.push_back(n);
+        normals.push_back(n);
+        normals.push_back(n);
+        verts.push_back(v1);
+        verts.push_back(v2);
+        verts.push_back(v3);
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+    }
+    glGenBuffers(1, &m_surfaceVbo);
+    glGenBuffers(1, &m_surfaceIbo);
+    glGenVertexArrays(1, &m_surfaceVao);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3 * 3, nullptr, GL_DYNAMIC_DRAW);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3, sizeof(double) * verts.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3 * 2, sizeof(double) * verts.size() * 3, static_cast<const void *>(forces.data()));
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_surfaceIbo);
+    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int) * 3 * faces.size(), static_cast<const void *>(faces.data()), GL_STATIC_DRAW);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+    glBindVertexArray(m_surfaceVao);
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 3, GL_DOUBLE, GL_FALSE, 0, static_cast<GLvoid *>(0));
+    glEnableVertexAttribArray(1);
+    glVertexAttribPointer(1, 3, GL_DOUBLE, GL_FALSE, 0, reinterpret_cast<GLvoid *>(sizeof(double) * verts.size() * 3));
+    glEnableVertexAttribArray(2);
+    glVertexAttribPointer(2, 3, GL_DOUBLE, GL_FALSE, 0, reinterpret_cast<GLvoid *>(sizeof(double) * verts.size() * 3 * 2));
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_surfaceIbo);
+    glBindVertexArray(0);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+    m_numSurfaceVertices = faces.size() * 3;
+    m_verticesSize = vertices.size();
+    m_faces = triangles;
+
+    if (vertices.size() > 4) { //shape
+        m_red = 0.93;
+        m_green = 0.8;
+        m_blue = 1.f;
+        m_alpha = 1.f;
+    } else { //ground
+        m_red = 1;
+        m_green = 1;
+        m_blue = 1;
+        m_alpha = 1.f;
+    }
+    m_force = 0;
+//    m_red = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
+//    m_blue = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
+//    m_green = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
+//    m_alpha = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
+}
+
+void Shape::setColor(float r, float g, float b) {
+    m_red = r;
+    m_green = g;
+    m_blue = b;
+}
+
+void Shape::init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3i> &triangles, const std::vector<Eigen::Vector4i> &tetIndices)
+{
+    init(vertices, triangles);
+
+    std::vector<Eigen::Vector2i> lines;
+    for(Vector4i tet : tetIndices) {
+        lines.emplace_back(tet[0], tet[1]);
+        lines.emplace_back(tet[0], tet[2]);
+        lines.emplace_back(tet[0], tet[3]);
+        lines.emplace_back(tet[1], tet[2]);
+        lines.emplace_back(tet[1], tet[3]);
+        lines.emplace_back(tet[2], tet[3]);
+    }
+    glGenBuffers(1, &m_tetVbo);
+    glGenBuffers(1, &m_tetIbo);
+    glGenVertexArrays(1, &m_tetVao);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_tetVbo);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3, vertices.data(), GL_DYNAMIC_DRAW);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_tetIbo);
+    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(int) * 2 * lines.size(), static_cast<const void *>(lines.data()), GL_STATIC_DRAW);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+    glBindVertexArray(m_tetVao);
+    glBindBuffer(GL_ARRAY_BUFFER, m_tetVbo);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 3, GL_DOUBLE, GL_FALSE, 0, static_cast<GLvoid *>(0));
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, m_tetIbo);
+    glBindVertexArray(0);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+
+    m_numTetVertices = lines.size() * 2;
+}
+
+void Shape::setVertices(const std::vector<Eigen::Vector3d> &vertices)
+{
+    if(vertices.size() != m_verticesSize) {
+        std::cerr << "You can't set vertices to a vector that is a different length that what shape was inited with" << std::endl;
+        return;
+    }
+    std::vector<Eigen::Vector3d> verts;
+    std::vector<Eigen::Vector3d> normals;
+    std::vector<Eigen::Vector3d> forces;
+    verts.reserve(m_faces.size() * 3);
+    normals.reserve(m_faces.size() * 3);
+    for(auto& f : m_faces) {
+        auto& v1 = vertices[f[0]];
+        auto& v2 = vertices[f[1]];
+        auto& v3 = vertices[f[2]];
+        auto& e1 = v2 - v1;
+        auto& e2 = v3 - v1;
+        auto n = e1.cross(e2);
+        normals.push_back(n);
+        normals.push_back(n);
+        normals.push_back(n);
+        verts.push_back(v1);
+        verts.push_back(v2);
+        verts.push_back(v3);
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+    }
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3, sizeof(double) * verts.size() * 3, static_cast<const void *>(normals.data()));
+    if(m_tetVao != static_cast<GLuint>(-1)) {
+        glBindBuffer(GL_ARRAY_BUFFER, m_tetVbo);
+        glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * vertices.size() * 3, static_cast<const void *>(vertices.data()));
+    }
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3 * 2, sizeof(double) * vertices.size() * 3 * 2, static_cast<const void *>(forces.data()));
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+void Shape::setVerticesF(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &forces) {
+    if(vertices.size() != m_verticesSize) {
+        std::cerr << "You can't set vertices to a vector that is a different length that what shape was inited with" << std::endl;
+        return;
+    }
+    if(vertices.size() != forces.size()) {
+        std::cerr << "Vertices and forces are not the same size" << std::endl;
+        return;
+    }
+    std::vector<Eigen::Vector3d> verts;
+    std::vector<Eigen::Vector3d> normals;
+    std::vector<Eigen::Vector3d> glForces;
+    verts.reserve(m_faces.size() * 3);
+    normals.reserve(m_faces.size() * 3);
+    glForces.reserve(m_faces.size() * 3);
+
+    double maxForceNorm = 500;
+    for(auto& f : m_faces) {
+        auto& v1 = vertices[f[0]];
+        auto& v2 = vertices[f[1]];
+        auto& v3 = vertices[f[2]];
+//        auto& f1 = forces[f[0]].normalized();
+//        auto& f2 = forces[f[1]].normalized();
+//        auto& f3 = forces[f[2]].normalized();
+        auto& f1 = forces[f[0]];
+        auto& f2 = forces[f[1]];
+        auto& f3 = forces[f[2]];
+        maxForceNorm = std::max(f1.norm(), maxForceNorm);
+        maxForceNorm = std::max(f2.norm(), maxForceNorm);
+        maxForceNorm = std::max(f3.norm(), maxForceNorm);
+        auto& e1 = v2 - v1;
+        auto& e2 = v3 - v1;
+        auto n = e1.cross(e2);
+        normals.push_back(n);
+        normals.push_back(n);
+        normals.push_back(n);
+        verts.push_back(v1);
+        verts.push_back(v2);
+        verts.push_back(v3);
+        glForces.push_back(f1);
+        glForces.push_back(f2);
+        glForces.push_back(f3); // Cool effect if it's v1, v2, v3 instead
+    }
+//    std::cout << maxForceNorm << std::endl;
+    for(Eigen::Vector3d &f : glForces) {
+        f /= maxForceNorm;
+        f = f.cwiseAbs();
+    }
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3, sizeof(double) * verts.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * verts.size() * 3 * 2, sizeof(double) * verts.size() * 3, static_cast<const void *>(glForces.data()));
+    if(m_tetVao != static_cast<GLuint>(-1)) {
+        glBindBuffer(GL_ARRAY_BUFFER, m_tetVbo);
+        glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * vertices.size() * 3, static_cast<const void *>(vertices.data()));
+    }
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+void Shape::setModelMatrix(const Eigen::Affine3f &model)
+{
+    m_modelMatrix = model.matrix();
+}
+
+void Shape::toggleWireframe()
+{
+    m_wireframe = !m_wireframe;
+}
+
+void Shape::toggleForce() {
+    m_force = std::abs(m_force - 1);
+}
+
+void Shape::setVertices(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &normals)
+{
+    std::vector<Eigen::Vector3d> forces;
+    if(vertices.size() != normals.size()) {
+        std::cerr << "Vertices and normals are not the same size" << std::endl;
+        return;
+    }
+    if(vertices.size() != m_verticesSize) {
+        std::cerr << "You can't set vertices to a vector that is a different length that what shape was inited with" << std::endl;
+        return;
+    }
+    for(Eigen::Vector3d v : vertices) {
+        forces.push_back(Eigen::Vector3d(1.0, 1.0, 0.0));
+    }
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(double) * vertices.size() * 3, static_cast<const void *>(vertices.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3, sizeof(double) * vertices.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(double) * vertices.size() * 3 * 2, sizeof(double) * vertices.size() * 3, static_cast<const void *>(forces.data()));
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+void Shape::draw(Shader *shader)
+{
+    Eigen::Matrix3f m3 = m_modelMatrix.topLeftCorner(3, 3);
+    Eigen::Matrix3f inverseTransposeModel = m3.inverse().transpose();
+
+    if(m_wireframe && m_tetVao != static_cast<GLuint>(-1)) {
+        shader->setUniform("wire", 1);
+        shader->setUniform("model", m_modelMatrix);
+        shader->setUniform("inverseTransposeModel", inverseTransposeModel);
+        shader->setUniform("red",   1);
+        shader->setUniform("green", 1);
+        shader->setUniform("blue",  1);
+        shader->setUniform("alpha", 1);
+        shader->setUniform("displayForce", m_force);
+        glBindVertexArray(m_tetVao);
+        glDrawElements(GL_LINES, m_numTetVertices, GL_UNSIGNED_INT, reinterpret_cast<GLvoid *>(0));
+        glBindVertexArray(0);
+    } else {
+        shader->setUniform("wire", 0);
+        shader->setUniform("model", m_modelMatrix);
+        shader->setUniform("inverseTransposeModel", inverseTransposeModel);
+        shader->setUniform("red",   m_red);
+        shader->setUniform("green", m_green);
+        shader->setUniform("blue",  m_blue);
+        shader->setUniform("alpha", m_alpha);
+        shader->setUniform("displayForce", m_force);
+        glBindVertexArray(m_surfaceVao);
+        glDrawElements(GL_TRIANGLES, m_numSurfaceVertices, GL_UNSIGNED_INT, reinterpret_cast<GLvoid *>(0));
+        glBindVertexArray(0);
+    }
+}
diff --git a/wave-sim/src/graphics/shape.h b/wave-sim/src/graphics/shape.h
new file mode 100644
index 0000000..6540ef6
--- /dev/null
+++ b/wave-sim/src/graphics/shape.h
@@ -0,0 +1,59 @@
+#ifndef SHAPE_H
+#define SHAPE_H
+
+#include <GL/glew.h>
+#include <vector>
+
+#include <Eigen/Dense>
+
+class Shader;
+
+class Shape
+{
+public:
+    Shape();
+
+//    void init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &normals, const std::vector<Eigen::Vector3i> &triangles);
+    void init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3i> &triangles);
+    void init(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3i> &triangles, const std::vector<Eigen::Vector4i> &tetIndices);
+
+    void setVertices(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &normals);
+    void setVertices(const std::vector<Eigen::Vector3d> &vertices);
+    void setVerticesF(const std::vector<Eigen::Vector3d> &vertices, const std::vector<Eigen::Vector3d> &forces);
+
+    void setModelMatrix(const Eigen::Affine3f &model);
+
+    void toggleWireframe();
+
+    void draw(Shader *shader);
+
+    void setColor(float r, float g, float b);
+
+    void toggleForce();
+
+private:
+    GLuint m_surfaceVao;
+    GLuint m_tetVao;
+    GLuint m_surfaceVbo;
+    GLuint m_tetVbo;
+    GLuint m_surfaceIbo;
+    GLuint m_tetIbo;
+
+    unsigned int m_numSurfaceVertices;
+    unsigned int m_numTetVertices;
+    unsigned int m_verticesSize;
+    float m_red;
+    float m_blue;
+    float m_green;
+    float m_alpha;
+
+    std::vector<Eigen::Vector3i> m_faces;
+
+    Eigen::Matrix4f m_modelMatrix;
+
+    bool m_wireframe;
+
+    int m_force;
+};
+
+#endif // SHAPE_H
diff --git a/wave-sim/src/main.cpp b/wave-sim/src/main.cpp
new file mode 100755
index 0000000..112c77b
--- /dev/null
+++ b/wave-sim/src/main.cpp
@@ -0,0 +1,38 @@
+#include "mainwindow.h"
+#include <cstdlib>
+#include <ctime>
+
+#include <QApplication>
+#include <QSurfaceFormat>
+#include <QScreen>
+
+int main(int argc, char *argv[])
+{
+    srand(static_cast<unsigned>(time(0)));
+
+    // Create a Qt application
+    QApplication a(argc, argv);
+    QCoreApplication::setApplicationName("Simulation");
+    QCoreApplication::setOrganizationName("CS 2240");
+    QCoreApplication::setApplicationVersion(QT_VERSION_STR);
+
+    // Set OpenGL version to 4.1 and context to Core
+    QSurfaceFormat fmt;
+    fmt.setVersion(4, 1);
+    fmt.setProfile(QSurfaceFormat::CoreProfile);
+    QSurfaceFormat::setDefaultFormat(fmt);
+
+    // Create a GUI window
+    MainWindow w;
+    w.resize(600, 500);
+    int desktopArea = QGuiApplication::primaryScreen()->size().width() *
+                      QGuiApplication::primaryScreen()->size().height();
+    int widgetArea = w.width() * w.height();
+    if (((float)widgetArea / (float)desktopArea) < 0.75f)
+        w.show();
+    else
+        w.showMaximized();
+
+
+    return a.exec();
+}
diff --git a/wave-sim/src/mainwindow.cpp b/wave-sim/src/mainwindow.cpp
new file mode 100755
index 0000000..72560dd
--- /dev/null
+++ b/wave-sim/src/mainwindow.cpp
@@ -0,0 +1,16 @@
+#include "mainwindow.h"
+#include <QHBoxLayout>
+
+MainWindow::MainWindow()
+{
+    glWidget = new GLWidget();
+
+    QHBoxLayout *container = new QHBoxLayout;
+    container->addWidget(glWidget);
+    this->setLayout(container);
+}
+
+MainWindow::~MainWindow()
+{
+    delete glWidget;
+}
diff --git a/wave-sim/src/mainwindow.h b/wave-sim/src/mainwindow.h
new file mode 100755
index 0000000..721c8f8
--- /dev/null
+++ b/wave-sim/src/mainwindow.h
@@ -0,0 +1,17 @@
+#pragma once
+
+#include <QMainWindow>
+#include "glwidget.h"
+
+class MainWindow : public QWidget
+{
+    Q_OBJECT
+
+public:
+    MainWindow();
+    ~MainWindow();
+
+private:
+
+    GLWidget *glWidget;
+};
diff --git a/wave-sim/src/mainwindow.ui b/wave-sim/src/mainwindow.ui
new file mode 100755
index 0000000..45c61e5
--- /dev/null
+++ b/wave-sim/src/mainwindow.ui
@@ -0,0 +1,46 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>MainWindow</class>
+ <widget class="QMainWindow" name="MainWindow">
+  <property name="geometry">
+   <rect>
+    <x>0</x>
+    <y>0</y>
+    <width>800</width>
+    <height>600</height>
+   </rect>
+  </property>
+  <property name="windowTitle">
+   <string>simulation</string>
+  </property>
+  <widget class="QWidget" name="centralWidget">
+   <layout class="QHBoxLayout" name="horizontalLayout">
+    <property name="leftMargin">
+     <number>0</number>
+    </property>
+    <property name="topMargin">
+     <number>0</number>
+    </property>
+    <property name="rightMargin">
+     <number>0</number>
+    </property>
+    <property name="bottomMargin">
+     <number>0</number>
+    </property>
+    <item>
+     <widget class="View" name="view" native="true"/>
+    </item>
+   </layout>
+  </widget>
+ </widget>
+ <layoutdefault spacing="6" margin="11"/>
+ <customwidgets>
+  <customwidget>
+   <class>View</class>
+   <extends>QWidget</extends>
+   <header>view.h</header>
+  </customwidget>
+ </customwidgets>
+ <resources/>
+ <connections/>
+</ui>
diff --git a/wave-sim/src/simulation.cpp b/wave-sim/src/simulation.cpp
new file mode 100644
index 0000000..409fa90
--- /dev/null
+++ b/wave-sim/src/simulation.cpp
@@ -0,0 +1,236 @@
+#include "simulation.h"
+#include "graphics/meshloader.h"
+
+#include <unordered_map>
+#include <unordered_set>
+#include <iostream>
+
+using namespace Eigen;
+
+Simulation::Simulation() {}
+
+int createFaceHash(int a, int b, int c, int n_vertices) {
+    int &low = a;
+    int &middle = b;
+    int &high = c;
+    if (low > middle)
+    {
+        std::swap(low, middle);
+    }
+    if (middle > high)
+    {
+        std::swap(middle, high);
+    }
+    if (low > middle)
+    {
+        std::swap(low, middle);
+    }
+
+    return (n_vertices * n_vertices * low) + (n_vertices * middle) + high;
+}
+
+Eigen::Vector3d Simulation::calculateFaceNormal(Eigen::Vector3d a, Eigen::Vector3d b, Eigen::Vector3d c) {
+    return (b - a).cross(c - a).normalized();
+}
+
+bool Simulation::calculatePointBehindNormal(Eigen::Vector3d a, Eigen::Vector3d b, Eigen::Vector3d c, Eigen::Vector3d p) {
+    Eigen::Vector3d a_to_p = p - a; // Calculates a direction from point on plane to point in question
+    return a_to_p.dot(calculateFaceNormal(a, b, c)) < 0;
+}
+
+void Simulation::init()
+{
+    // STUDENTS: This code loads up the tetrahedral mesh in 'example-meshes/single-tet.mesh'
+    //    (note: your working directory must be set to the root directory of the starter code
+    //    repo for this file to load correctly). You'll probably want to instead have this code
+    //    load up a tet mesh based on e.g. a file path specified with a command line argument.
+    std::vector<Vector3d> vertices;
+    std::vector<Vector4i> tets;
+    this->mainSystem = System();
+    if (MeshLoader::loadTetMesh("./example-meshes/sphere.mesh", vertices, tets)) {
+        // STUDENTS: This code computes the surface mesh of the loaded tet mesh, i.e. the faces
+        //    of tetrahedra which are on the exterior surface of the object. Right now, this is
+        //    hard-coded for the single-tet mesh. You'll need to implement surface mesh extraction
+        //    for arbitrary tet meshes. Think about how you can identify which tetrahedron faces
+        //    are surface faces...
+        std::vector<Vector3i> faces;
+//        std::unordered_map<int, std::pair<Eigen::Vector3i*, int>> includedFaces;
+//        std::unordered_map<int, int> includedFaces;
+        std::unordered_map<int, Eigen::Vector3i> includedFaces;
+        // includes the currently parsed faces and their index
+        std::vector<Vector4i> faceOrders;
+
+        faceOrders.emplace_back(1, 0, 2, 3); // first three are the included points, last is the excluded point
+        faceOrders.emplace_back(2, 0, 3, 1);
+        faceOrders.emplace_back(3, 1, 2, 0);
+        faceOrders.emplace_back(3, 0, 1, 2);
+
+        for(Vector4i t : tets) {
+            for(Vector4i fo : faceOrders) {
+                int hash = createFaceHash(t[fo[0]], t[fo[1]], t[fo[2]], vertices.size()); // Finding which vertex indexes the face has and ordering them from least to smallest
+                if(includedFaces.contains(hash)) {
+                    includedFaces.erase(hash);
+                } else {
+                    Eigen::Vector3i orderedFace;
+                    Vector3d a = vertices.at(t[fo[0]]);
+                    Vector3d b = vertices.at(t[fo[1]]);
+                    Vector3d c = vertices.at(t[fo[2]]);
+                    Vector3d d = vertices.at(t[fo[3]]);
+                    if(calculatePointBehindNormal(a, b, c, d)) {
+                        orderedFace = Eigen::Vector3i(t[fo[0]], t[fo[1]], t[fo[2]]); // Wind it backwards if the excluded point is behind the face
+                    } else {
+                        orderedFace = Eigen::Vector3i(t[fo[2]], t[fo[1]], t[fo[0]]); // Wind it backwards if the excluded point is in front of the face
+                    }
+
+                    includedFaces.emplace(hash, orderedFace);
+                }
+            }
+        }
+        for (const auto & [ key, value ] : includedFaces) {
+            faces.push_back(value);
+        }
+        m_shape.init(vertices, faces, tets);
+        mainSystem.initFromVecs(vertices, tets);
+    }
+    m_shape.setModelMatrix(Affine3f(Eigen::Translation3f(0, 2, 0)));
+
+    initGround();
+    initExtra();
+}
+
+void Simulation::update(double seconds)
+{
+    // STUDENTS: This method should contain all the time-stepping logic for your simulation.
+    //   Specifically, the code you write here should compute new, updated vertex positions for your
+    //   simulation mesh, and it should then call m_shape.setVertices to update the display with those
+    //   newly-updated vertices.
+
+    // STUDENTS: As currently written, the program will just continually compute simulation timesteps as long
+    //    as the program is running (see View::tick in view.cpp) . You might want to e.g. add a hotkey for pausing
+    //    the simulation, and perhaps start the simulation out in a paused state.
+
+    // Note that the "seconds" parameter represents the amount of time that has passed since
+    // the last update
+    if(this->estimationMode == EULER) {
+        mainSystem.updateCalculations();
+        mainSystem.updatePositions(seconds);
+        mainSystem.resolveCollisions();
+    }
+    else if (this->estimationMode == MIDPOINT) {
+        mainSystem.updateCalculations();
+        std::vector<Eigen::Vector3d> originalPositions = mainSystem.getPositions();
+        std::vector<Eigen::Vector3d> originalVelocities = mainSystem.getVelocities();
+        mainSystem.updatePositions(seconds / 2);
+        mainSystem.updateCalculations();
+        mainSystem.setVelocities(originalVelocities);
+        mainSystem.setPositions(originalPositions);
+        mainSystem.updatePositions(seconds);
+        mainSystem.resolveCollisions();
+    //    m_shape.setVertices(mainSystem.getNodePos());
+        m_shape.setVerticesF(mainSystem.getNodePos(), mainSystem.getNodeForces());
+    } else if (this->estimationMode == ADAPTIVE) {
+        double errorTolerance = 1e-4;
+
+        std::vector<Eigen::Vector3d> originalPositions = mainSystem.getPositions();
+        std::vector<Eigen::Vector3d> originalVelocities = mainSystem.getVelocities();
+
+        mainSystem.updateCalculations();
+        mainSystem.updatePositions(seconds);
+        Eigen::VectorXd state1 = mainSystem.getState();
+
+        mainSystem.setVelocities(originalVelocities);
+        mainSystem.setPositions(originalPositions);
+        mainSystem.updatePositions(seconds / 2);
+        mainSystem.updateCalculations();
+        mainSystem.updatePositions(seconds / 2);
+        Eigen::VectorXd state2 = mainSystem.getState();
+
+        double newStepsize = seconds * std::sqrt(errorTolerance / (state1 - state2).norm());
+
+        mainSystem.setVelocities(originalVelocities);
+        mainSystem.setPositions(originalPositions);
+        mainSystem.updateCalculations();
+        mainSystem.updatePositions(newStepsize);
+        mainSystem.resolveCollisions();
+        m_shape.setVerticesF(mainSystem.getNodePos(), mainSystem.getNodeForces());
+    }
+}
+
+void Simulation::draw(Shader *shader)
+{
+    m_shape.draw(shader);
+    m_ground.draw(shader);
+    m_extra.draw(shader);
+}
+
+void Simulation::toggleWire()
+{
+    m_shape.toggleWireframe();
+}
+
+void Simulation::toggleForceRender() {
+    m_shape.toggleForce();
+}
+
+void Simulation::initGround()
+{
+    std::vector<Vector3d> groundVerts;
+    std::vector<Vector3i> groundFaces;
+    groundVerts.emplace_back(-5, 0, -5);
+    groundVerts.emplace_back(-5, 0, 5);
+    groundVerts.emplace_back(5, 0, 5);
+    groundVerts.emplace_back(5, 0, -5);
+    groundFaces.emplace_back(0, 1, 2);
+    groundFaces.emplace_back(0, 2, 3);
+    m_ground.init(groundVerts, groundFaces);
+}
+
+void Simulation::initExtra()
+{
+    std::vector<Vector3d> extraVerts;
+    std::vector<Vector4i> extraTets;
+    Eigen::Vector3d pos = Eigen::Vector3d(0, 0, 0);
+    if (MeshLoader::loadTetMesh("./example-meshes/sphere.mesh", extraVerts, extraTets)) {
+
+        std::vector<Vector3i> faces;
+        std::unordered_map<int, Eigen::Vector3i> includedFaces;
+        std::vector<Vector4i> faceOrders;
+
+        faceOrders.emplace_back(1, 0, 2, 3); // first three are the included points, last is the excluded point
+        faceOrders.emplace_back(2, 0, 3, 1);
+        faceOrders.emplace_back(3, 1, 2, 0);
+        faceOrders.emplace_back(3, 0, 1, 2);
+
+        for(Eigen::Vector3d& ev : extraVerts) {
+            ev += pos;
+        }
+
+        for(Vector4i t : extraTets) {
+            for(Vector4i fo : faceOrders) {
+                int hash = createFaceHash(t[fo[0]], t[fo[1]], t[fo[2]], extraVerts.size()); // Finding which vertex indexes the face has and ordering them from least to smallest
+                if(includedFaces.contains(hash)) {
+                    includedFaces.erase(hash);
+                } else {
+                    Eigen::Vector3i orderedFace;
+                    Vector3d a = extraVerts.at(t[fo[0]]);
+                    Vector3d b = extraVerts.at(t[fo[1]]);
+                    Vector3d c = extraVerts.at(t[fo[2]]);
+                    Vector3d d = extraVerts.at(t[fo[3]]);
+                    if(calculatePointBehindNormal(a, b, c, d)) {
+                        orderedFace = Eigen::Vector3i(t[fo[0]], t[fo[1]], t[fo[2]]); // Wind it backwards if the excluded point is behind the face
+                    } else {
+                        orderedFace = Eigen::Vector3i(t[fo[2]], t[fo[1]], t[fo[0]]); // Wind it backwards if the excluded point is in front of the face
+                    }
+
+                    includedFaces.emplace(hash, orderedFace);
+                }
+            }
+        }
+        for (const auto & [ key, value ] : includedFaces) {
+//            std::cout << key << ": " << value << std::endl;
+            faces.push_back(value);
+        }
+        m_extra.init(extraVerts, faces);
+        m_extra.setColor(0.9, 0.8, 0.1);
+    }
+}
diff --git a/wave-sim/src/simulation.h b/wave-sim/src/simulation.h
new file mode 100644
index 0000000..68f2178
--- /dev/null
+++ b/wave-sim/src/simulation.h
@@ -0,0 +1,42 @@
+#pragma once
+
+#include "graphics/shape.h"
+#include "system.h"
+
+class Shader;
+
+class Simulation
+{
+public:
+    Simulation();
+
+    void init();
+
+    void update(double seconds);
+
+    void draw(Shader *shader);
+
+    void toggleWire();
+
+    void toggleForceRender();
+
+    Eigen::Vector3d calculateFaceNormal(Eigen::Vector3d a, Eigen::Vector3d b, Eigen::Vector3d c);
+
+    bool calculatePointBehindNormal(Eigen::Vector3d a, Eigen::Vector3d b, Eigen::Vector3d c, Eigen::Vector3d p);
+
+    System mainSystem;
+
+    enum EstimationMode { EULER, MIDPOINT, ADAPTIVE };
+
+private:
+    Shape m_shape;
+
+    Shape m_ground;
+
+    Shape m_extra;
+
+    EstimationMode estimationMode = MIDPOINT;
+
+    void initGround();
+    void initExtra();
+};
diff --git a/wave-sim/src/system.cpp b/wave-sim/src/system.cpp
new file mode 100644
index 0000000..f0a4894
--- /dev/null
+++ b/wave-sim/src/system.cpp
@@ -0,0 +1,262 @@
+#include "system.h"
+#include <iostream>
+
+System::System()
+{
+    nodes = std::vector<Node*>();
+    tets = std::vector<Tet*>();
+}
+
+void System::initFromVecs(std::vector<Eigen::Vector3d> v, std::vector<Eigen::Vector4i> t) {
+    nodes = std::vector<Node*>();
+    tets = std::vector<Tet*>();
+
+    for(Eigen::Vector3d vertPos : v) {
+        nodes.push_back( new Node { vertPos } );
+    }
+    for(Eigen::Vector4i tetVerts : t) {
+        std::vector<Node*> definingNodes;
+        definingNodes.push_back(nodes.at(tetVerts[0]));
+        definingNodes.push_back(nodes.at(tetVerts[1]));
+        definingNodes.push_back(nodes.at(tetVerts[2]));
+        definingNodes.push_back(nodes.at(tetVerts[3]));
+        tets.push_back( new Tet(definingNodes) );
+    }
+    this->setParams( Params { } );
+}
+
+std::vector<Eigen::Vector3d> System::getPositions() {
+    std::vector<Eigen::Vector3d> res;
+    for(int i = 0 ; i < nodes.size(); ++i) {
+        res.push_back(nodes.at(i)->pos);
+    }
+
+    return res;
+}
+
+void System::setPositions(std::vector<Eigen::Vector3d> positions) {
+    assert(positions.size() == nodes.size());
+    for(int i = 0 ; i < positions.size(); ++i) {
+        nodes.at(i)->pos = positions.at(i);
+    }
+}
+
+std::vector<Eigen::Vector3d> System::getVelocities() {
+    std::vector<Eigen::Vector3d> res;
+    for(int i = 0 ; i < nodes.size(); ++i) {
+        res.push_back(nodes.at(i)->vel);
+    }
+
+    return res;
+}
+
+void System::setVelocities(std::vector<Eigen::Vector3d> velocities) {
+    assert(velocities.size() == nodes.size());
+    for(int i = 0 ; i < velocities.size(); ++i) {
+        nodes.at(i)->vel = velocities.at(i);
+    }
+}
+
+System::Face::Face(Node* a, Node* b, Node* c, Node* p){
+    area = (b->pos - a->pos).cross(c->pos - a->pos).norm() / 2;
+    normal = (b->pos - a->pos).cross(c->pos - a->pos).normalized();
+    if(normal.dot(p->pos - a->pos) >= 0) {
+        // If the normal faces towards the back point, flip it and wind nodes other way.
+        normal = -normal;
+        nodes.push_back(c);
+        nodes.push_back(b);
+        nodes.push_back(a);
+    }
+    else {
+        nodes.push_back(a);
+        nodes.push_back(b);
+        nodes.push_back(c);
+    }
+    opp = p;
+}
+
+void System::Face::calculateAreaNorms() {
+    Node* a = nodes.at(0);
+    Node* b = nodes.at(1);
+    Node* c = nodes.at(2);
+    this->area = (b->pos - a->pos).cross(c->pos - a->pos).norm() / 2;
+    this->normal = (b->pos - a->pos).cross(c->pos - a->pos).normalized();
+}
+
+System::Tet::Tet(std::vector<Node*> in_nodes)  {
+    nodes = in_nodes;
+
+    Eigen::Vector3d p0 = in_nodes.at(0)->pos;
+    Eigen::Vector3d p1 = in_nodes.at(1)->pos;
+    Eigen::Vector3d p2 = in_nodes.at(2)->pos;
+    Eigen::Vector3d p3 = in_nodes.at(3)->pos;
+
+    Node* a = nodes.at(0);
+    Node* b = nodes.at(1);
+    Node* c = nodes.at(2);
+    Node* d = nodes.at(3); // Origin at d
+
+    // Face 0 defined by 1, 0, 2
+    // Face 1 defined by 2, 0, 3
+    // Face 2 defined by 3, 1, 2
+    // Face 3 defined by 3, 0, 1
+    faces.push_back(new Face(b, a, c, d));
+    faces.push_back(new Face(c, a, d, b));
+    faces.push_back(new Face(d, b, c, a));
+    faces.push_back(new Face(d, a, b, c));
+
+    betaMat << (a->pos - d->pos), (b->pos - d->pos), (c->pos - d->pos);
+
+    posMat = betaMat;
+    betaMat = betaMat.inverse().eval();
+    this->velMat << (a->vel - d->vel), (b->vel - d->vel), (c->vel - d->vel);
+
+    Eigen::Matrix4d massDeterminer;
+
+    massDeterminer << a->pos[0], a->pos[1], a->pos[2], 1,
+                      b->pos[0], b->pos[1], b->pos[2], 1,
+                      c->pos[0], c->pos[1], c->pos[2], 1,
+                      d->pos[0], d->pos[1], d->pos[2], 1;
+
+    mass = this->_rho * (massDeterminer.determinant()) / 6.;
+
+    for(Node* n : this->nodes) {
+        n->mass += mass / 4;
+    }
+
+    for(Face* f : this->faces) {
+        f->calculateAreaNorms();
+//            totalFaceArea += f->area;
+    }
+}
+
+void System::updateCalculations() {
+    for(Node* n : nodes) {
+        n->force = Eigen::Vector3d::Zero();
+        n->force = this->_gravity * n->mass;
+    }
+
+    for(Tet* t : tets) {
+        t->update();
+    }
+}
+
+void System::resolveCollisions() {
+    double groundLevel = -2;
+    double radius = 1;
+
+    for(Node* n : nodes) {
+        if(n->pos[1] < groundLevel && std::abs(n->pos[0]) < 5 && std::abs(n->pos[2]) < 5) {
+//            n->vel[1] *= -1/(4) * (n->pos[1] + 2);
+            n->pos[1] = groundLevel + (groundLevel - n->pos[1]);
+            n->vel[0] /= this->groundTraction;
+            n->vel[1] = std::abs(n->vel[1]) / this->groundAbsorbance;
+        }
+        if(n->pos[0] * n->pos[0] + n->pos[2] * n->pos[2] + std::pow(n->pos[1] + 2, 2) < radius * radius) {
+            Eigen::Vector3d normal = (n->pos - Eigen::Vector3d(0, -2, 0)).normalized();
+//            n->pos = 2 * normal * radius + Eigen::Vector3d(0, -2, 0) - (n->pos - Eigen::Vector3d(0, -2, 0));
+            n->pos = normal * radius + Eigen::Vector3d(0, -2, 0);
+            Eigen::Vector3d normalComponent = n->vel.dot(normal) * normal;
+            Eigen::Vector3d tangentComponent = n->vel - normalComponent;
+            n->vel = -(normalComponent / this->groundAbsorbance) - tangentComponent / this->groundTraction;
+        }
+    }
+}
+
+void System::setParams(Params params) {
+    for(Tet* t : tets) {
+        t->_lambda = params._lambda;
+        t->_mu = params._mu;
+        t->_rho = params._rho;
+        t->_phi = params._phi;
+        t->_psi = params._psi;
+    }
+    this->_gravity = params.gravity;
+    this->groundAbsorbance = params._absorbance;
+    this->groundTraction = params._traction;
+}
+
+void System::updatePositions(double deltaTime) {
+//    double totalFaceArea = 0;
+    for(Node* n : nodes) {
+        n->acc = n->force * (1 / n->mass);
+        n->vel += n->acc * deltaTime;
+        n->pos += n->vel * deltaTime;
+    }
+//    this->resolveCollisions();
+//    std::cout << totalFaceArea << std::endl;
+}
+
+void System::Tet::update() {
+//    for(Face* f : this->faces) {
+//        f->calculateAreaNorms();
+////            totalFaceArea += f->area;
+//    }
+
+    Node* a = nodes.at(0);
+    Node* b = nodes.at(1);
+    Node* c = nodes.at(2);
+    Node* d = nodes.at(3); // Origin at d
+
+    this->posMat << (a->pos - d->pos), (b->pos - d->pos), (c->pos - d->pos);
+    this->velMat << (a->vel - d->vel), (b->vel - d->vel), (c->vel - d->vel);
+
+    this->FMat = posMat * betaMat;
+
+//    Eigen::Matrix3d test;
+//    test << Eigen::Vector3d(1, 2, 3), Eigen::Vector3d(4, 5, 6), Eigen::Vector3d(7, 8, 9);
+
+    if(this->FMat.isApprox(Eigen::Matrix3d::Identity())) {
+        this->FMat = Eigen::Matrix3d::Identity();
+    }
+
+    this->strain = (FMat.transpose() * FMat) - Eigen::Matrix3d::Identity();
+    this->strainRate = ((posMat * betaMat).transpose() * (velMat * betaMat)) + ((velMat * betaMat).transpose() * (posMat * betaMat));
+
+    Eigen::Matrix3d stressFromStrain = (this->_lambda * Eigen::Matrix3d::Identity() * strain.trace()) + (2 * this->_mu * strain);
+    Eigen::Matrix3d stressFromStrainRate = (this->_phi * Eigen::Matrix3d::Identity() * strainRate.trace()) + (2 * this->_psi * strainRate);
+    this->stress = stressFromStrain + stressFromStrainRate;
+
+//    for(Node* n : this->nodes) {
+////        n->force = Eigen::Vector3d(0., -.1, 0.) * n->mass;
+////        n->force = Eigen::Vector3d(0, 0, 0);
+////        n->acc = Ei
+//    }
+
+    for(Face* f : this->faces) {
+        Eigen::Vector3d faceForce = this->FMat * this->stress * f->area * f->normal;
+        for(Node* n : f->nodes) {
+            n->force += faceForce / 3;
+        }
+//        f->opp->force -= faceForce;
+    }
+}
+
+std::vector<Eigen::Vector3d> System::getNodePos() {
+    std::vector<Eigen::Vector3d> res;
+    for(Node* n : nodes) {
+        res.push_back(n->pos);
+    }
+    return res;
+}
+
+std::vector<Eigen::Vector3d> System::getNodeForces() {
+    std::vector<Eigen::Vector3d> res;
+    for(Node* n : nodes) {
+        res.push_back(n->force);
+    }
+    return res;
+}
+
+Eigen::VectorXd System::getState() {
+    std::vector<double> state;
+    for(Node* n : this->nodes) {
+        state.push_back(n->pos[0]);
+        state.push_back(n->pos[1]);
+        state.push_back(n->pos[2]);
+        state.push_back(n->vel[0]);
+        state.push_back(n->vel[1]);
+        state.push_back(n->vel[2]);
+    }
+    return Eigen::Map<Eigen::VectorXd>(state.data(), state.size());
+}
diff --git a/wave-sim/src/system.h b/wave-sim/src/system.h
new file mode 100644
index 0000000..289aca3
--- /dev/null
+++ b/wave-sim/src/system.h
@@ -0,0 +1,99 @@
+#pragma once
+
+#include <Eigen/Dense>
+#include <vector>
+
+class System
+{
+public:
+    System();
+
+    struct Node;
+    struct Tet;
+
+    Eigen::Vector3d _gravity = Eigen::Vector3d(0., -1.0, 0.);
+    double groundAbsorbance = 4e4;
+    double groundTraction = 4e4;
+
+    std::vector<Node*> nodes;
+    std::vector<Tet*> tets;
+
+    std::vector<Eigen::Vector3d> getPositions();
+    void setPositions(std::vector<Eigen::Vector3d> positions);
+
+    std::vector<Eigen::Vector3d> getVelocities();
+    void setVelocities(std::vector<Eigen::Vector3d> velocities);
+
+    void initFromVecs(std::vector<Eigen::Vector3d> v, std::vector<Eigen::Vector4i> t);
+
+    struct Params {
+        const Eigen::Vector3d gravity = Eigen::Vector3d(0, -1, 0);
+        const double _lambda = 1e4; //incompressibility for the whole material
+        const double _mu = 1e4; //rigidity for the whole material
+        const double _rho = 1200.; //density
+        const double _phi = 1e1; //coefficients of viscosity
+        const double _psi = 1e1;
+        const double _traction = 4e4;
+        const double _absorbance = 1e1;
+    };
+
+    void setParams(Params params);
+
+    struct Node {
+        Eigen::Vector3d pos;
+        double mass = 0;
+
+        Eigen::Vector3d vel = Eigen::Vector3d(0., 0., 0.);
+        Eigen::Vector3d force = Eigen::Vector3d(0., 0., 0.);
+        Eigen::Vector3d acc = Eigen::Vector3d(0., 0., 0.);
+    };
+
+    struct Face {
+        double area;
+        Eigen::Vector3d normal;
+
+        std::vector<Node*> nodes;
+        Node* opp;
+
+        // P is unused point we use to see if normal is facing right way.
+        Face(Node* a, Node* b, Node* c, Node* p);
+
+        void calculateAreaNorms();
+    };
+
+    struct Tet {
+        double _lambda = 1e4; //incompressibility for the whole material
+        double _mu = 1e4; //rigidity for the whole material
+        double _rho = 1200.; //density
+        double _phi = 1e1; //coefficients of viscosity
+        double _psi = 1e1;
+
+        std::vector<Node*> nodes;
+//        Eigen::Vector4d faceAreas;
+//        std::vector<Eigen::Vector3d> normals;
+        std::vector<Face*> faces;
+        Eigen::Matrix3d betaMat;
+        Eigen::Matrix3d posMat;
+        Eigen::Matrix3d velMat;
+
+        Eigen::Matrix3d FMat;
+        Eigen::Matrix3d strain;
+        Eigen::Matrix3d strainRate;
+        Eigen::Matrix3d stress;
+
+        double mass;
+
+        Tet(std::vector<Node*> in_nodes);
+
+        void update();
+    };
+
+    void updateCalculations();
+    void resolveCollisions();
+    void updatePositions(double deltaTime);
+
+    std::vector<Eigen::Vector3d> getNodePos();
+    std::vector<Eigen::Vector3d> getNodeForces();
+
+    Eigen::VectorXd getState();
+};