merge with new foam fixes

6 files changed, 437 insertions, 67 deletions

diff --git a/src/arap.cpp b/src/arap.cpp
index 670f06e..2d952ae 100644
--- a/src/arap.cpp
+++ b/src/arap.cpp
@@ -66,7 +66,7 @@ void ARAP::init
 //	minCorner = Vector3f(-1.0f, -1.0f, -1.0f);
 //	maxCorner = Vector3f(1.0f, 1.0f, 1.0f);
 
-    initCausticsShape(10);
+    initCausticsShape(1000);
 }
 
 void ARAP::update(double seconds)
@@ -85,7 +85,7 @@ void ARAP::update(double seconds)
 	m_ocean.fft_prime(m_time);
     m_ocean.update_ocean();
     m_shape.setVertices_and_Normals(m_ocean.get_vertices(), m_ocean.getNormals());
-  //  m_shape.setVertices(m_ocean.get_vertices());
+//      m_shape.setVertices(m_ocean.get_vertices());
 
 //	auto tmp = m_ocean.get_vertices();
 	// print the min and max of the vertices
@@ -94,7 +94,7 @@ void ARAP::update(double seconds)
 //	for (int i = 0; i < tmp.size(); i++) {
 //		min = min.cwiseMin(tmp[i]);
 //		max = max.cwiseMax(tmp[i]);
-//	}
+//	}w
 //	std::cout << "min: " << min << std::endl;
 //std::cout << "max: " << max << std::endl;
 
diff --git a/src/glwidget.cpp b/src/glwidget.cpp
index 8f3c0c1..f8547a6 100755
--- a/src/glwidget.cpp
+++ b/src/glwidget.cpp
@@ -250,7 +250,7 @@ void GLWidget::paintCaustics() {
 
 //    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 
-    glBindFramebuffer(GL_FRAMEBUFFER, m_fbo_texture);
+    glBindFramebuffer(GL_FRAMEBUFFER, m_fbo);
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
     glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
     glEnable( GL_BLEND );
diff --git a/src/graphics/shape.cpp b/src/graphics/shape.cpp
index 1baa233..7073a7f 100644
--- a/src/graphics/shape.cpp
+++ b/src/graphics/shape.cpp
@@ -336,7 +336,7 @@ void Shape::updateMesh(const std::vector<Eigen::Vector3i> &faces,
         Vector3f n = getNormal(face);
 
         for (auto& v: {face[0], face[1], face[2]}) {
-            normals.push_back(Eigen::Vector3f(1,1,1));
+            normals.push_back(n);
             verts.push_back(vertices[v]);
 
             if (m_anchors.find(v) == m_anchors.end()) {
diff --git a/src/graphics/simpleshape.cpp b/src/graphics/simpleshape.cpp
new file mode 100644
index 0000000..567dab3
--- /dev/null
+++ b/src/graphics/simpleshape.cpp
@@ -0,0 +1,272 @@
+#include "simpleshape.h"
+
+#include <iostream>
+#include "graphics/shader.h"
+
+using namespace Eigen;
+using namespace std;
+
+// ================== Constructor
+
+SimpleShape::SimpleShape() :
+    m_surfaceVao(),
+    m_surfaceVbo(),
+    m_surfaceIbo(),
+    m_numSurfaceVertices(),
+    m_verticesSize(),
+    m_red(),
+    m_blue(),
+    m_green(),
+    m_alpha(),
+    m_faces(),
+    m_vertices(),
+    m_anchors(),
+    m_modelMatrix(Matrix4f::Identity()),
+    lastSelected(-1)
+{}
+
+// ================== Initialization and Updating
+
+void SimpleShape::init(const vector<Vector3f> &vertices, const vector<Vector3i> &triangles)
+{
+    m_vertices.clear();
+    copy(vertices.begin(), vertices.end(), back_inserter(m_vertices));
+
+    vector<Vector3f> verts;
+    vector<Vector3f> normals;
+    vector<Vector3f> colors;
+    vector<Vector3i> faces;
+    faces.reserve(triangles.size());
+
+    for (int s = 0; s < triangles.size() * 3; s+=3) faces.push_back(Vector3i(s, s + 1, s + 2));
+    updateMesh(triangles, vertices, verts, normals, colors);
+
+    glGenBuffers(1, &m_surfaceVbo);
+    glGenBuffers(1, &m_surfaceIbo);
+    glGenVertexArrays(1, &m_surfaceVao);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3) + (colors.size() * 3)), nullptr, GL_DYNAMIC_DRAW);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(float) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * verts.size() * 3, sizeof(float) * normals.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3)), sizeof(float) * colors.size() * 3, static_cast<const void *>(colors.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_FLOAT, GL_FALSE, 0, static_cast<GLvoid *>(0));
+    glEnableVertexAttribArray(1);
+    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, reinterpret_cast<GLvoid *>(sizeof(float) * verts.size() * 3));
+    glEnableVertexAttribArray(2);
+    glVertexAttribPointer(2, 3, GL_FLOAT, GL_FALSE, 0, reinterpret_cast<GLvoid *>(sizeof(float) * (verts.size() * 3 + colors.size() * 3)));
+
+    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;
+    m_red   = 0.5f + 0.5f * rand() / ((float) RAND_MAX);
+    m_blue  = 0.5f + 0.5f * rand() / ((float) RAND_MAX);
+    m_green = 0.5f + 0.5f * rand() / ((float) RAND_MAX);
+    m_alpha = 1.0f;
+}
+
+void SimpleShape::setVertices(const vector<Vector3f> &vertices)
+{
+    m_vertices.clear();
+    copy(vertices.begin(), vertices.end(), back_inserter(m_vertices));
+
+    vector<Vector3f> verts;
+    vector<Vector3f> normals;
+    vector<Vector3f> colors;
+
+    updateMesh(m_faces, vertices, verts, normals, colors);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3) + (colors.size() * 3)), nullptr, GL_DYNAMIC_DRAW);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(float) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * verts.size() * 3, sizeof(float) * normals.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3)), sizeof(float) * colors.size() * 3, static_cast<const void *>(colors.data()));
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+// ================== Model Matrix
+
+void SimpleShape::setModelMatrix(const Affine3f &model) { m_modelMatrix = model.matrix(); }
+
+// ================== General Graphics Stuff
+
+void SimpleShape::draw(Shader *shader, GLenum mode)
+{
+    Eigen::Matrix3f m3 = m_modelMatrix.topLeftCorner(3, 3);
+    Eigen::Matrix3f inverseTransposeModel = m3.inverse().transpose();
+
+    switch(mode) {
+    case GL_TRIANGLES:
+    {
+        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);
+        glBindVertexArray(m_surfaceVao);
+        glDrawElements(mode, m_numSurfaceVertices, GL_UNSIGNED_INT, reinterpret_cast<GLvoid *>(0));
+        glBindVertexArray(0);
+        break;
+    }
+    case GL_POINTS:
+    {
+        shader->setUniform("model", m_modelMatrix);
+        shader->setUniform("inverseTransposeModel", inverseTransposeModel);
+        glBindVertexArray(m_surfaceVao);
+        glDrawElements(mode, m_numSurfaceVertices, GL_UNSIGNED_INT, reinterpret_cast<GLvoid *>(0));
+        glBindVertexArray(0);
+        break;
+    }
+    }
+}
+
+SelectMode SimpleShape::select(Shader *shader, int closest_vertex)
+{
+    if (closest_vertex == -1) return SelectMode::None;
+
+    bool vertexIsNowSelected = m_anchors.find(closest_vertex) == m_anchors.end();
+
+    if (vertexIsNowSelected) {
+        m_anchors.insert(closest_vertex);
+    } else {
+        m_anchors.erase(closest_vertex);
+    }
+
+    selectHelper();
+
+    return vertexIsNowSelected ? SelectMode::Anchor : SelectMode::Unanchor;
+}
+
+bool SimpleShape::selectWithSpecifiedMode(Shader *shader, int closest_vertex, SelectMode mode)
+{
+    switch (mode) {
+    case SelectMode::None: {
+        return false;
+    }
+    case SelectMode::Anchor: {
+        if (m_anchors.find(closest_vertex) != m_anchors.end()) return false;
+        m_anchors.insert(closest_vertex);
+        break;
+    }
+    case SelectMode::Unanchor: {
+        if (m_anchors.find(closest_vertex) == m_anchors.end()) return false;
+        m_anchors.erase(closest_vertex);
+        break;
+    }
+    }
+
+    selectHelper();
+
+    return true;
+}
+
+int SimpleShape::getClosestVertex(Vector3f start, Vector3f ray, float threshold)
+{
+    int closest_vertex = -1;
+    int i = 0;
+    float dist = numeric_limits<float>::max();
+    ParametrizedLine line = ParametrizedLine<float, 3>::Through(start, start + ray);
+
+    for (const Vector3f &v : m_vertices) {
+        float d = line.distance(v);
+        if (d<dist) {
+            dist = d;
+            closest_vertex = i;
+        }
+        ++i;
+    }
+
+    if (dist >= threshold) closest_vertex = -1;
+
+    return closest_vertex;
+}
+
+bool SimpleShape::getAnchorPos(int lastSelected,
+                         Eigen::Vector3f& pos,
+                         Eigen::Vector3f  ray,
+                         Eigen::Vector3f  start)
+{
+    bool isAnchor = m_anchors.find(lastSelected) != m_anchors.end();
+    if (isAnchor) {
+        Eigen::Vector3f oldPos = m_vertices[lastSelected];
+        Eigen::ParametrizedLine line = ParametrizedLine<float, 3>::Through(start, start+ray);
+        pos = line.projection(oldPos);
+    }
+    return isAnchor;
+}
+
+// ================== Accessors
+
+const vector<Vector3f>   &SimpleShape::getVertices() { return m_vertices; }
+const vector<Vector3i>   &SimpleShape::getFaces()    { return m_faces;    }
+const unordered_set<int> &SimpleShape::getAnchors()  { return m_anchors;  }
+
+// ================== Helpers
+
+void SimpleShape::selectHelper()
+{
+    vector<Vector3f> verts;
+    vector<Vector3f> normals;
+    vector<Vector3f> colors;
+    updateMesh(m_faces, m_vertices, verts, normals, colors);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_surfaceVbo);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3) + (colors.size() * 3)), nullptr, GL_DYNAMIC_DRAW);
+    glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(float) * verts.size() * 3, static_cast<const void *>(verts.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * verts.size() * 3, sizeof(float) * normals.size() * 3, static_cast<const void *>(normals.data()));
+    glBufferSubData(GL_ARRAY_BUFFER, sizeof(float) * ((verts.size() * 3) + (normals.size() * 3)), sizeof(float) * colors.size() * 3, static_cast<const void *>(colors.data()));
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+Vector3f SimpleShape::getNormal(const Vector3i& face)
+{
+    Vector3f& v1 = m_vertices[face[0]];
+    Vector3f& v2 = m_vertices[face[1]];
+    Vector3f& v3 = m_vertices[face[2]];
+    Vector3f e1 = v2 - v1;
+    Vector3f e2 = v3 - v1;
+    Vector3f n = e1.cross(e2);
+    return n.normalized();
+}
+
+void SimpleShape::updateMesh(const std::vector<Eigen::Vector3i> &faces,
+                       const std::vector<Eigen::Vector3f> &vertices,
+                       std::vector<Eigen::Vector3f>& verts,
+                       std::vector<Eigen::Vector3f>& normals,
+                       std::vector<Eigen::Vector3f>& colors)
+{
+    verts.reserve(faces.size() * 3);
+    normals.reserve(faces.size() * 3);
+
+    for (const Eigen::Vector3i& face : faces) {
+        Vector3f n = getNormal(face);
+
+        for (auto& v: {face[0], face[1], face[2]}) {
+            normals.push_back(n);
+            verts.push_back(vertices[v]);
+
+            if (m_anchors.find(v) == m_anchors.end()) {
+                colors.push_back(Vector3f(1,0,0));
+            } else {
+                colors.push_back(Vector3f(0, 1 - m_green, 1 - m_blue));
+            }
+        }
+    }
+}
+
diff --git a/src/graphics/simpleshape.h b/src/graphics/simpleshape.h
new file mode 100644
index 0000000..03db726
--- /dev/null
+++ b/src/graphics/simpleshape.h
@@ -0,0 +1,89 @@
+#pragma once
+
+#include <GL/glew.h>
+#include <set>
+#include <vector>
+#include <unordered_set>
+
+#define EIGEN_DONT_VECTORIZE
+#define EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
+//#include "Eigen/StdVector"
+//EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(Eigen::Matrix2f)
+//EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(Eigen::Matrix3f)
+//EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(Eigen::Matrix3i)
+#include "Eigen/Dense"
+
+enum SelectMode
+{
+    None     = 0,
+    Anchor   = 1,
+    Unanchor = 2
+};
+
+
+
+
+struct Edge{
+    std::pair<int, int> ij;
+    float weight = 0.f;
+    std::set<int> opposite_two;
+};
+
+struct Cell{
+    int center_index = 0; // this is i
+    std::set<int> neighbor_js; // this is neighboring j's
+    std::set<int> edge_ids; // this is ids if ij edges in edge vector
+};
+
+
+class Shader;
+
+class SimpleShape
+{
+public:
+    SimpleShape();
+
+    void init(const std::vector<Eigen::Vector3f> &vertices, const std::vector<Eigen::Vector3i> &triangles);
+    void setVertices(const std::vector<Eigen::Vector3f> &vertices);
+
+    void setModelMatrix(const Eigen::Affine3f &model);
+
+    void draw(Shader *shader, GLenum mode);
+    SelectMode select(Shader *shader, int vertex);
+    bool selectWithSpecifiedMode(Shader *shader, int vertex, SelectMode mode);
+    int  getClosestVertex(Eigen::Vector3f start, Eigen::Vector3f ray, float threshold);
+    bool getAnchorPos(int lastSelected, Eigen::Vector3f& pos, Eigen::Vector3f ray, Eigen::Vector3f start);
+
+    const std::vector<Eigen::Vector3f>& getVertices();
+    const std::vector<Eigen::Vector3i>& getFaces();
+    const std::unordered_set<int>& getAnchors();
+
+private:
+    GLuint m_surfaceVao;
+    GLuint m_surfaceVbo;
+    GLuint m_surfaceIbo;
+
+    unsigned int m_numSurfaceVertices;
+    unsigned int m_verticesSize;
+    float m_red;
+    float m_blue;
+    float m_green;
+    float m_alpha;
+
+    std::vector<Eigen::Vector3i> m_faces;
+    std::vector<Eigen::Vector3f> m_vertices;
+    std::unordered_set<int>      m_anchors;
+
+    Eigen::Matrix4f m_modelMatrix;
+    int lastSelected = -1;
+
+    // Helpers
+
+    void selectHelper();
+    Eigen::Vector3f getNormal(const Eigen::Vector3i& face);
+    void updateMesh(const std::vector<Eigen::Vector3i> &triangles,
+                    const std::vector<Eigen::Vector3f> &vertices,
+                           std::vector<Eigen::Vector3f>& verts,
+                           std::vector<Eigen::Vector3f>& normals,
+                           std::vector<Eigen::Vector3f>& colors);
+};
diff --git a/src/ocean/ocean_alt.cpp b/src/ocean/ocean_alt.cpp
index 58e76f8..1c9306b 100644
--- a/src/ocean/ocean_alt.cpp
+++ b/src/ocean/ocean_alt.cpp
@@ -316,17 +316,6 @@ Eigen::Vector2d ocean_alt::complex_exp(double exponent){
 void ocean_alt::update_ocean()
 {
     std::vector<Eigen::Vector3f> vertices = std::vector<Eigen::Vector3f>();
-    if (iterations < 10){
-        for (int i = 0; i < N; i++){
-            Eigen::Vector2d amplitude = m_current_h[i];
-            float height = amplitude[0];
-
-            if (height < min) min = height;
-            if (height > max) max = height;
-
-        }
-        iterations ++;
-    }
 	// reset normals & vertices arrays for the single tile
     m_vertices = std::vector<Eigen::Vector3f>(N);
 	m_normals = std::vector<Eigen::Vector3f>(N);
@@ -336,7 +325,18 @@ void ocean_alt::update_ocean()
         Eigen::Vector2d amplitude = m_current_h[i];
         float height = amplitude[0];
 
-        // Eigen::Vector2d slope = m_slopes[i] * .3f;
+		if (iterations++ > 1)
+		{
+			if (height < min) min = height;
+			if (height > max)
+			{
+				max = height;
+//				std::cout << "changed!! max: " << max << std::endl;
+			}
+		}
+
+
+		// Eigen::Vector2d slope = m_slopes[i] * .3f;
         // Eigen::Vector3f s = Eigen::Vector3f(-slope[0], 0.0, -slope[1]);
         // Eigen::Vector3f y = Eigen::Vector3f(0.0, 1.0, 0.0);
 
@@ -368,10 +368,18 @@ void ocean_alt::update_ocean()
 
 
        // m_foam_constants.wavelengths[i] = 2.f* M_PI * m_slopes[i].dot(m_slopes[i]) /  Lx;
-        float h_0 = m_waveIndexConstants[i].h0_prime[0]; // min*.2f;
-        float h_max = max*.001f; // the smaller the constant, the more foam there is
-        float waveheight = (height - h_0 ) / (h_max - h_0);
-        m_foam_constants.wavelengths[i] = waveheight;
+//        float h_0 = m_waveIndexConstants[i].h0_prime[0]; // min*.2f;
+//        float h_max = max*.001f; // the smaller the constant, the more foam there is
+//        float waveheight = (height - h_0 ) / (h_max - h_0);
+//        m_foam_constants.wavelengths[i] = waveheight;
+        float h_0 = 0; // min*.2f;
+        float h_max = max*.35f; // the smaller the constant, the more foam there is
+        m_foam_constants.wavelengths[i] = (height - h_0 ) / (h_max - h_0);
+
+//        if (i < 5){
+//            std::cout << h_0 << ", " << h_max << std::endl;
+//            std::cout << m_foam_constants.wavelengths[i] << std::endl;
+//        }
 
         if (waveheight >= height_threshold){
                  //std::cout << "push" << std::endl;
@@ -381,6 +389,7 @@ void ocean_alt::update_ocean()
 
     }
 
+
     // populate foam constants
     m_foam_constants.positions = vertices;
 }
@@ -427,54 +436,54 @@ std::vector<Eigen::Vector3i> ocean_alt::get_faces()
 {
     // connect the vertices into faces
     std::vector<Eigen::Vector3i> faces = std::vector<Eigen::Vector3i>();
-	for (int i = 0; i < num_tiles_x; i++)
-	{
-		for (int j = 0; j < num_tiles_z; j++)
-		{
-			for (int k = 0; k < N; k++)
-			{
-				int x = k % num_rows;
-				int z = k / num_rows;
-
-				// connect the vertices into faces
-				if (x < num_rows - 1 && z < num_cols - 1)
-				{
-					int tile_index_offset = (j + num_tiles_z * i) * N;
-					int i1 = k + tile_index_offset;
-					int i2 = k + 1 + tile_index_offset;
-					int i3 = k + num_rows + tile_index_offset;
-					int i4 = k + num_rows + 1 + tile_index_offset;
-
-					faces.emplace_back(i2, i1, i3);
-					faces.emplace_back(i2, i3, i4);
-				}
-			}
-		}
-	}
-
-	return faces;
-
-
-//    for (int i = 0; i < N; i++)
-//    {
-//        int x = i / num_rows;
-//        int z = i % num_rows;
+//	for (int i = 0; i < num_tiles_x; i++)
+//	{
+//		for (int j = 0; j < num_tiles_z; j++)
+//		{
+//			for (int k = 0; k < N; k++)
+//			{
+//				int x = k % num_rows;
+//				int z = k / num_rows;
 //
-//        // connect the vertices into faces
-//        if (x < num_rows - 1 && z < num_cols - 1)
-//        {
-//            int i1 = i;
-//            int i2 = i + 1;
-//            int i3 = i + num_rows;
-//            int i4 = i + num_rows + 1;
+//				// connect the vertices into faces
+//				if (x < num_rows - 1 && z < num_cols - 1)
+//				{
+//					int tile_index_offset = (j + num_tiles_z * i) * N;
+//					int i1 = k + tile_index_offset;
+//					int i2 = k + 1 + tile_index_offset;
+//					int i3 = k + num_rows + tile_index_offset;
+//					int i4 = k + num_rows + 1 + tile_index_offset;
 //
-//            faces.emplace_back(i2, i1, i3);
-//            faces.emplace_back(i2, i3, i4);
-//			faces.emplace_back(i1, i2, i3);
-//			faces.emplace_back(i3, i2, i4);
-//        }
-//    }
-//    return faces;
+//					faces.emplace_back(i2, i1, i3);
+//					faces.emplace_back(i2, i3, i4);
+//				}
+//			}
+//		}
+//	}
+//
+//	return faces;
+
+
+    for (int i = 0; i < N; i++)
+    {
+        int x = i / num_rows;
+        int z = i % num_rows;
+
+        // connect the vertices into faces
+        if (x < num_rows - 1 && z < num_cols - 1)
+        {
+            int i1 = i;
+            int i2 = i + 1;
+            int i3 = i + num_rows;
+            int i4 = i + num_rows + 1;
+
+            faces.emplace_back(i2, i1, i3);
+            faces.emplace_back(i2, i3, i4);
+			faces.emplace_back(i1, i2, i3);
+			faces.emplace_back(i3, i2, i4);
+        }
+    }
+    return faces;
 }
 
 Eigen::Vector2d muliply_complex(Eigen::Vector2d a, Eigen::Vector2d b)