initial

1 files changed, 236 insertions, 0 deletions

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);
+    }
+}