Add cloth simulation

Integrate FastMassSpring, implemented by Samer Itani(@sam007961),
Original paper: "Fast Simulation of Mass-Spring Systems" by Liu, T., Bargteil, A. W., Obrien, J. F., & Kavan, L.

ACKNOWLEDGEMENTS.html

@@ -1272,3 +1272,28 @@ https://www.reddit.com/r/gamedev/comments/5iuf3h/i_am_writting_a_3d_monster_mode
     other computer software, distribute, and sublicense such enhancements or
     derivative works thereof, in binary and source code form.
 </pre>
+
+<h1>FastMassSpring</h1>
+<pre>
+    MIT License
+
+    Copyright (c) 2019 Samer Itani
+
+    Permission is hereby granted, free of charge, to any person obtaining a copy
+    of this software and associated documentation files (the "Software"), to deal
+    in the Software without restriction, including without limitation the rights
+    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+    copies of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be included in all
+    copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+    SOFTWARE.
+</pre>

dust3d.pro

@@ -482,10 +482,20 @@ HEADERS += src/contourtopartconverter.h
 SOURCES += src/remesher.cpp
 HEADERS += src/remesher.h
 
+SOURCES += src/clothsimulator.cpp
+HEADERS += src/clothsimulator.h
+
+SOURCES += src/componentlayer.cpp
+HEADERS += src/componentlayer.h
+
 SOURCES += src/main.cpp
 
 HEADERS += src/version.h
 
+INCLUDEPATH += thirdparty/FastMassSpring/ClothApp
+SOURCES += thirdparty/FastMassSpring/ClothApp/MassSpringSolver.cpp
+HEADERS += thirdparty/FastMassSpring/ClothApp/MassSpringSolver.h
+
 INCLUDEPATH += thirdparty/instant-meshes
 INCLUDEPATH += thirdparty/instant-meshes/instant-meshes-dust3d/src
 INCLUDEPATH += thirdparty/instant-meshes/instant-meshes-dust3d/ext/tbb/include

languages/dust3d_zh_CN.ts

@@ -748,6 +748,10 @@ Tips:
         <source>Poly</source>
         <translation>面数</translation>
     </message>
+    <message>
+        <source>Layer</source>
+        <translation>层</translation>
+    </message>
 </context>
 <context>
     <name>PartWidget</name>
@@ -1070,6 +1074,14 @@ Tips:
         <source>Extreme High Poly</source>
         <translation>极高面</translation>
     </message>
+    <message>
+        <source>Body</source>
+        <translation>身体</translation>
+    </message>
+    <message>
+        <source>Cloth</source>
+        <translation>衣服</translation>
+    </message>
 </context>
 <context>
     <name>RigWidget</name>

src/clothsimulator.cpp

@@ -0,0 +1,240 @@
+#include <MassSpringSolver.h>
+#include <set>
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Polyhedron_incremental_builder_3.h>
+#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
+#include <CGAL/AABB_face_graph_triangle_primitive.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+#include "clothsimulator.h"
+#include "booleanmesh.h"
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef K::Point_3 Point;
+typedef K::Triangle_3 Triangle;
+typedef CGAL::Polyhedron_3<K> Polyhedron;
+typedef Polyhedron::HalfedgeDS HalfedgeDS;
+typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
+typedef CGAL::AABB_traits<K, Primitive> Traits;
+typedef CGAL::AABB_tree<Traits> Tree;
+typedef CGAL::Side_of_triangle_mesh<Polyhedron, K> Point_inside;
+
+template <class HDS>
+class Build_mesh : public CGAL::Modifier_base<HDS> {
+public:
+    Build_mesh(const std::vector<QVector3D> *vertices,
+            const std::vector<std::vector<size_t>> *faces) :
+        m_vertices(vertices),
+        m_faces(faces)
+    {
+    };
+    void operator()(HDS& hds)
+    {
+        // Postcondition: hds is a valid polyhedral surface.
+        CGAL::Polyhedron_incremental_builder_3<HDS> B(hds, false);
+        B.begin_surface(m_vertices->size(), m_faces->size());
+        typedef typename HDS::Vertex   Vertex;
+        typedef typename Vertex::Point Point;
+        for (const auto &it: *m_vertices)
+            B.add_vertex(Point(it.x(), it.y(), it.z()));
+        for (const auto &it: *m_faces) {
+            B.begin_facet();
+            B.add_vertex_to_facet(it[0]);
+            B.add_vertex_to_facet(it[1]);
+            B.add_vertex_to_facet(it[2]);
+            B.end_facet();
+        }
+        B.end_surface();
+    };
+private:
+    const std::vector<QVector3D> *m_vertices = nullptr;
+    const std::vector<std::vector<size_t>> *m_faces = nullptr;
+};
+
+// System parameters
+namespace SystemParam {
+    static const int n = 61; // must be odd, n * n = n_vertices | 61
+    static const float w = 2.0f; // width | 2.0f
+    static const float h = 0.008f; // time step, smaller for better results | 0.008f = 0.016f/2
+    static const float r = w / (n - 1) * 1.05f; // spring rest legnth
+    static const float k = 1.0f; // spring stiffness | 1.0f;
+    static const float m = 0.25f / (n * n); // point mass | 0.25f
+    static const float a = 0.993f; // damping, close to 1.0 | 0.993f
+    static const float g = 9.8f * m; // gravitational force | 9.8f
+}
+
+// Point - mesh collision node
+class CgMeshCollisionNode : public CgPointNode {
+private:
+    Tree *m_aabbTree = nullptr;
+    Point_inside *m_insideTester = nullptr;
+    Polyhedron m_polyhedron;
+public:
+    CgMeshCollisionNode(mass_spring_system *system, float *vbuff,
+            const std::vector<QVector3D> &collisionVertices,
+            const std::vector<std::vector<size_t>> &collisionTriangles) :
+        CgPointNode(system, vbuff)
+    {
+        if (!collisionTriangles.empty()) {
+            Build_mesh<HalfedgeDS> mesh(&collisionVertices, &collisionTriangles);
+            m_polyhedron.delegate(mesh);
+            m_aabbTree = new Tree(faces(m_polyhedron).first, faces(m_polyhedron).second, m_polyhedron);
+            m_aabbTree->accelerate_distance_queries();
+            m_insideTester = new Point_inside(*m_aabbTree);
+        }
+    }
+    
+    ~CgMeshCollisionNode()
+    {
+        delete m_insideTester;
+        delete m_aabbTree;
+    };
+    
+    bool query(unsigned int i) const
+    {
+        return false;
+    };
+    
+    void satisfy()
+    {
+        for (unsigned int i = 0; i < system->n_points; i++) {
+            auto offset = 3 * i;
+            Point point(vbuff[offset + 0],
+                vbuff[offset + 1],
+                vbuff[offset + 2]);
+            if (nullptr != m_insideTester &&
+                    (*m_insideTester)(point) != CGAL::ON_UNBOUNDED_SIDE) {
+                Point closestPoint = m_aabbTree->closest_point(point);
+                vbuff[offset + 0] = closestPoint.x();
+                vbuff[offset + 1] = closestPoint.y();
+                vbuff[offset + 2] = closestPoint.z();
+            }
+        }
+    }
+};
+
+ClothSimulator::ClothSimulator(const std::vector<QVector3D> &vertices,
+        const std::vector<std::vector<size_t>> &faces,
+        const std::vector<QVector3D> &collisionVertices,
+        const std::vector<std::vector<size_t>> &collisionTriangles) :
+    m_vertices(vertices),
+    m_faces(faces),
+    m_collisionVertices(collisionVertices),
+    m_collisionTriangles(collisionTriangles)
+{
+}
+
+ClothSimulator::~ClothSimulator()
+{
+    delete m_massSpringSystem;
+    delete m_massSpringSolver;
+    delete m_rootNode;
+    delete m_deformationNode;
+    delete m_meshCollisionNode;
+}
+
+void ClothSimulator::convertMeshToCloth()
+{
+    m_clothPointSources.reserve(m_vertices.size());
+    m_clothPointBuffer.reserve(m_vertices.size() * 3);
+
+    std::map<size_t, size_t> oldVertexToNewMap;
+    auto addPoint = [&](size_t index) {
+        auto findNew = oldVertexToNewMap.find(index);
+        if (findNew != oldVertexToNewMap.end())
+            return findNew->second;
+        const auto &position = m_vertices[index];
+        m_clothPointBuffer.push_back(position.x());
+        m_clothPointBuffer.push_back(position.y());
+        m_clothPointBuffer.push_back(position.z());
+        size_t newIndex = m_clothPointSources.size();
+        m_clothPointSources.push_back(index);
+        oldVertexToNewMap.insert({index, newIndex});
+        return newIndex;
+    };
+
+    std::set<std::pair<size_t, size_t>> oldEdges;
+    for (const auto &it: m_faces) {
+        for (size_t i = 0; i < it.size(); ++i) {
+            size_t j = (i + 1) % it.size();
+            if (oldEdges.find(std::make_pair(it[i], it[j])) != oldEdges.end())
+                continue;
+            m_clothSprings.push_back({addPoint(it[i]), addPoint(it[j])});
+            oldEdges.insert({it[i], it[j]});
+            oldEdges.insert({it[j], it[i]});
+        }
+    }
+}
+
+void ClothSimulator::getCurrentVertices(std::vector<QVector3D> *currentVertices)
+{
+    *currentVertices = m_vertices;
+    for (size_t newIndex = 0; newIndex < m_clothPointSources.size(); ++newIndex) {
+        size_t oldIndex = m_clothPointSources[newIndex];
+        auto offset = newIndex * 3;
+        (*currentVertices)[oldIndex] = QVector3D(m_clothPointBuffer[offset + 0],
+            m_clothPointBuffer[offset + 1],
+            m_clothPointBuffer[offset + 2]);
+    }
+}
+
+void ClothSimulator::step()
+{
+    if (nullptr == m_massSpringSolver)
+        return;
+        
+    m_massSpringSolver->solve(5);
+    m_massSpringSolver->solve(5);
+    
+    CgSatisfyVisitor visitor;
+    visitor.satisfy(*m_rootNode);
+}
+
+void ClothSimulator::create()
+{
+    convertMeshToCloth();
+    
+    if (m_clothPointSources.empty())
+        return;
+    
+    mass_spring_system::VectorXf masses(SystemParam::m * mass_spring_system::VectorXf::Ones((unsigned int)m_clothSprings.size()));
+    
+    mass_spring_system::EdgeList springList(m_clothSprings.size());
+    mass_spring_system::VectorXf restLengths(m_clothSprings.size());
+    mass_spring_system::VectorXf stiffnesses(m_clothSprings.size());
+    
+    for (size_t i = 0; i < m_clothSprings.size(); ++i) {
+        const auto &source = m_clothSprings[i];
+        springList[i] = mass_spring_system::Edge(source.first, source.second);
+        restLengths[i] = (m_vertices[m_clothPointSources[source.first]] - m_vertices[m_clothPointSources[source.second]]).length();
+        stiffnesses[i] = SystemParam::k;
+    }
+    
+    mass_spring_system::VectorXf fext = Eigen::Vector3f(0, -SystemParam::g, 0).replicate(m_clothPointSources.size(), 1);
+    
+    m_massSpringSystem = new mass_spring_system(m_clothPointSources.size(), m_clothSprings.size(), SystemParam::h, springList, restLengths,
+        stiffnesses, masses, fext, SystemParam::a);
+
+    m_massSpringSolver = new MassSpringSolver(m_massSpringSystem, m_clothPointBuffer.data());
+    
+    // deformation constraint parameters
+    const float tauc = 0.12f; // critical spring deformation | 0.12f
+    const unsigned int deformIter = 15; // number of iterations | 15
+    
+    std::vector<unsigned int> structSprintIndexList;
+    structSprintIndexList.reserve(springList.size());
+    m_deformationNode =
+        new CgSpringDeformationNode(m_massSpringSystem, m_clothPointBuffer.data(), tauc, deformIter);
+    m_deformationNode->addSprings(structSprintIndexList);
+    
+    m_rootNode = new CgRootNode(m_massSpringSystem, m_clothPointBuffer.data());
+    m_rootNode->addChild(m_deformationNode);
+    
+    m_meshCollisionNode = new CgMeshCollisionNode(m_massSpringSystem, m_clothPointBuffer.data(),
+        m_collisionVertices,
+        m_collisionTriangles);
+    m_rootNode->addChild(m_meshCollisionNode);
+}

src/clothsimulator.h

@@ -0,0 +1,41 @@
+#ifndef DUST3D_CLOTH_SIMULATOR_H
+#define DUST3D_CLOTH_SIMULATOR_H
+#include <QObject>
+#include <QVector3D>
+#include <vector>
+
+struct mass_spring_system;
+class MassSpringSolver;
+class CgRootNode;
+class CgSpringDeformationNode;
+class CgMeshCollisionNode;
+
+class ClothSimulator : public QObject
+{
+    Q_OBJECT
+public:
+    ClothSimulator(const std::vector<QVector3D> &vertices,
+        const std::vector<std::vector<size_t>> &faces,
+        const std::vector<QVector3D> &collisionVertices,
+        const std::vector<std::vector<size_t>> &collisionTriangles);
+    ~ClothSimulator();
+    void create();
+    void step();
+    void getCurrentVertices(std::vector<QVector3D> *currentVertices);
+private:
+    std::vector<QVector3D> m_vertices;
+    std::vector<std::vector<size_t>> m_faces;
+    std::vector<QVector3D> m_collisionVertices;
+    std::vector<std::vector<size_t>> m_collisionTriangles;
+    std::vector<float> m_clothPointBuffer;
+    std::vector<size_t> m_clothPointSources;
+    std::vector<std::pair<size_t, size_t>> m_clothSprings;
+    mass_spring_system *m_massSpringSystem = nullptr;
+    MassSpringSolver *m_massSpringSolver = nullptr;
+    CgRootNode *m_rootNode = nullptr;
+    CgSpringDeformationNode *m_deformationNode = nullptr;
+    CgMeshCollisionNode *m_meshCollisionNode = nullptr;
+    void convertMeshToCloth();
+};
+
+#endif

src/componentlayer.cpp

@@ -0,0 +1,6 @@
+#include <QObject>
+#include "componentlayer.h"
+
+IMPL_ComponentLayerToString
+IMPL_ComponentLayerFromString
+IMPL_ComponentLayerToDispName

src/componentlayer.h

@@ -0,0 +1,49 @@
+#ifndef DUST3D_COMPONENT_LAYER_H
+#define DUST3D_COMPONENT_LAYER_H
+#include <QString>
+
+enum class ComponentLayer
+{
+    Body = 0,
+    Cloth,
+    Count
+};
+ComponentLayer ComponentLayerFromString(const char *layerString);
+#define IMPL_ComponentLayerFromString                                       \
+ComponentLayer ComponentLayerFromString(const char *layerString)            \
+{                                                                           \
+    QString layer = layerString;                                            \
+    if (layer == "Body")                                                    \
+        return ComponentLayer::Body;                                        \
+    if (layer == "Cloth")                                                   \
+        return ComponentLayer::Cloth;                                       \
+    return ComponentLayer::Body;                                            \
+}
+const char *ComponentLayerToString(ComponentLayer layer);
+#define IMPL_ComponentLayerToString                                         \
+const char *ComponentLayerToString(ComponentLayer layer)                    \
+{                                                                           \
+    switch (layer) {                                                        \
+        case ComponentLayer::Body:                                          \
+            return "Body";                                                  \
+        case ComponentLayer::Cloth:                                         \
+            return "Cloth";                                                 \
+        default:                                                            \
+            return "Body";                                                  \
+    }                                                                       \
+}
+QString ComponentLayerToDispName(ComponentLayer layer);
+#define IMPL_ComponentLayerToDispName                                       \
+QString ComponentLayerToDispName(ComponentLayer layer)                      \
+{                                                                           \
+    switch (layer) {                                                        \
+        case ComponentLayer::Body:                                          \
+            return QObject::tr("Body");                                     \
+        case ComponentLayer::Cloth:                                         \
+            return QObject::tr("Cloth");                                    \
+        default:                                                            \
+            return QObject::tr("Body");                                     \
+    }                                                                       \
+}
+
+#endif

src/document.cpp

@@ -1201,6 +1201,8 @@ void Document::toSnapshot(Snapshot *snapshot, const std::set<QUuid> &limitNodeId
                 component["smoothSeam"] = QString::number(componentIt.second.smoothSeam);
             if (componentIt.second.polyCount != PolyCount::Original)
                 component["polyCount"] = PolyCountToString(componentIt.second.polyCount);
+            if (componentIt.second.layer != ComponentLayer::Body)
+                component["layer"] = ComponentLayerToString(componentIt.second.layer);
             QStringList childIdList;
             for (const auto &childId: componentIt.second.childrenIds) {
                 childIdList.append(childId.toString());
@@ -1708,6 +1710,7 @@ void Document::addFromSnapshot(const Snapshot &snapshot, bool fromPaste)
         if (smoothSeamIt != componentKv.second.end())
             component.setSmoothSeam(smoothSeamIt->second.toFloat());
         component.polyCount = PolyCountFromString(valueOfKeyInMapOrEmpty(componentKv.second, "polyCount").toUtf8().constData());
+        component.layer = ComponentLayerFromString(valueOfKeyInMapOrEmpty(componentKv.second, "layer").toUtf8().constData());
         //qDebug() << "Add component:" << component.id << " old:" << componentKv.first << "name:" << component.name;
         if ("partId" == linkDataType) {
             QUuid partId = oldNewIdMap[QUuid(linkData)];
@@ -2497,6 +2500,20 @@ void Document::setComponentPolyCount(QUuid componentId, PolyCount count)
     emit skeletonChanged();
 }
 
+void Document::setComponentLayer(QUuid componentId, ComponentLayer layer)
+{
+    Component *component = (Component *)findComponent(componentId);
+    if (nullptr == component)
+        return;
+    if (component->layer == layer)
+        return;
+    
+    component->layer = layer;
+    component->dirty = true;
+    emit componentLayerChanged(componentId);
+    emit skeletonChanged();
+}
+
 void Document::createNewComponentAndMoveThisIn(QUuid componentId)
 {
     auto component = componentMap.find(componentId);

src/document.h

@@ -30,6 +30,7 @@
 #include "preferences.h"
 #include "paintmode.h"
 #include "proceduralanimation.h"
+#include "componentlayer.h"
 
 class MaterialPreviewsGenerator;
 class MotionsGenerator;
@@ -67,6 +68,7 @@ public:
     float smoothAll = 0.0;
     float smoothSeam = 0.0;
     PolyCount polyCount = PolyCount::Original;
+    ComponentLayer layer = ComponentLayer::Body;
     std::vector<QUuid> childrenIds;
     QString linkData() const
     {
@@ -396,6 +398,7 @@ signals:
     void componentSmoothAllChanged(QUuid componentId);
     void componentSmoothSeamChanged(QUuid componentId);
     void componentPolyCountChanged(QUuid componentId);
+    void componentLayerChanged(QUuid componentId);
     void nodeRemoved(QUuid nodeId);
     void edgeRemoved(QUuid edgeId);
     void nodeRadiusChanged(QUuid nodeId);
@@ -642,6 +645,7 @@ public slots:
     void setComponentSmoothAll(QUuid componentId, float toSmoothAll);
     void setComponentSmoothSeam(QUuid componentId, float toSmoothSeam);
     void setComponentPolyCount(QUuid componentId, PolyCount count);
+    void setComponentLayer(QUuid componentId, ComponentLayer layer);
     void hideOtherComponents(QUuid componentId);
     void lockOtherComponents(QUuid componentId);
     void hideAllComponents();

src/documentwindow.cpp

@@ -1014,6 +1014,7 @@ DocumentWindow::DocumentWindow() :
     connect(partTreeWidget, &PartTreeWidget::setComponentSmoothAll, m_document, &Document::setComponentSmoothAll);
     connect(partTreeWidget, &PartTreeWidget::setComponentSmoothSeam, m_document, &Document::setComponentSmoothSeam);
     connect(partTreeWidget, &PartTreeWidget::setComponentPolyCount, m_document, &Document::setComponentPolyCount);
+    connect(partTreeWidget, &PartTreeWidget::setComponentLayer, m_document, &Document::setComponentLayer);
     connect(partTreeWidget, &PartTreeWidget::moveComponent, m_document, &Document::moveComponent);
     connect(partTreeWidget, &PartTreeWidget::removeComponent, m_document, &Document::removeComponent);
     connect(partTreeWidget, &PartTreeWidget::hideOtherComponents, m_document, &Document::hideOtherComponents);

src/meshgenerator.cpp

@@ -18,6 +18,7 @@
 #include "triangulatefaces.h"
 #include "remesher.h"
 #include "polycount.h"
+#include "clothsimulator.h"
 
 MeshGenerator::MeshGenerator(Snapshot *snapshot) :
     m_snapshot(snapshot)
@@ -840,6 +841,11 @@ CombineMode MeshGenerator::componentCombineMode(const std::map<QString, QString>
             combineMode = CombineMode::Inversion;
         if (componentRemeshed(component))
             combineMode = CombineMode::Uncombined;
+        if (combineMode == CombineMode::Normal) {
+            if (ComponentLayer::Body != ComponentLayerFromString(valueOfKeyInMapOrEmpty(*component, "layer").toUtf8().constData())) {
+                combineMode = CombineMode::Uncombined;
+            }
+        }
     }
     return combineMode;
 }
@@ -879,6 +885,13 @@ QString MeshGenerator::componentColorName(const std::map<QString, QString> *comp
     return QString();
 }
 
+ComponentLayer MeshGenerator::componentLayer(const std::map<QString, QString> *component)
+{
+    if (nullptr == component)
+        return ComponentLayer::Body;
+    return ComponentLayerFromString(valueOfKeyInMapOrEmpty(*component, "layer").toUtf8().constData());
+}
+
 MeshCombiner::Mesh *MeshGenerator::combineComponentMesh(const QString &componentIdString, CombineMode *combineMode)
 {
     MeshCombiner::Mesh *mesh = nullptr;
@@ -1035,6 +1048,16 @@ MeshCombiner::Mesh *MeshGenerator::combineComponentMesh(const QString &component
         mesh = nullptr;
     }
     
+    if (componentId.isNull()) {
+        // Prepare cloth collision shap
+        if (nullptr != mesh && !mesh->isNull()) {
+            mesh->fetch(m_clothCollisionVertices, m_clothCollisionTriangles);
+        } else {
+            // TODO: when no body is valid, may add ground plane as collision shape
+            // ... ...
+        }
+    }
+    
     if (nullptr != mesh) {
         float polyCountValue = 1.0f;
         bool remeshed = componentId.isNull() ? componentRemeshed(&m_snapshot->canvas, &polyCountValue) : componentRemeshed(component, &polyCountValue);
@@ -1386,6 +1409,7 @@ void MeshGenerator::generate()
     
     // Recursively check uncombined components
     collectUncombinedComponent(QUuid().toString());
+    collectClothComponent(QUuid().toString());
     
     // Collect errored parts
     for (const auto &it: m_cacheContext->parts) {
@@ -1514,6 +1538,8 @@ void MeshGenerator::collectUncombinedComponent(const QString &componentIdString)
 {
     const auto &component = findComponent(componentIdString);
     if (CombineMode::Uncombined == componentCombineMode(component)) {
+        if (ComponentLayer::Body != componentLayer(component))
+            return;
         const auto &componentCache = m_cacheContext->components[componentIdString];
         if (nullptr == componentCache.mesh || componentCache.mesh->isNull()) {
             qDebug() << "Uncombined mesh is null";
@@ -1553,6 +1579,71 @@ void MeshGenerator::collectUncombinedComponent(const QString &componentIdString)
     }
 }
 
+void MeshGenerator::collectClothComponent(const QString &componentIdString)
+{
+    const auto &component = findComponent(componentIdString);
+    if (ComponentLayer::Cloth == componentLayer(component)) {
+        const auto &componentCache = m_cacheContext->components[componentIdString];
+        if (nullptr == componentCache.mesh || componentCache.mesh->isNull()) {
+            return;
+        }
+        if (m_clothCollisionTriangles.empty())
+            return;
+        
+        std::vector<QVector3D> uncombinedVertices;
+        std::vector<std::vector<size_t>> uncombinedFaces;
+        componentCache.mesh->fetch(uncombinedVertices, uncombinedFaces);
+        
+        std::map<PositionKey, std::pair<QUuid, QUuid>> positionMap;
+        for (const auto &it: componentCache.outcomeNodeVertices) {
+            positionMap.insert({PositionKey(it.first), it.second});
+        }
+        std::vector<std::pair<QUuid, QUuid>> uncombinedVertexSources(uncombinedVertices.size());
+        for (size_t i = 0; i < uncombinedVertices.size(); ++i) {
+            auto findSource = positionMap.find(PositionKey(uncombinedVertices[i]));
+            if (findSource == positionMap.end())
+                continue;
+            uncombinedVertexSources[i] = findSource->second;
+        }
+        
+        ClothSimulator clothSimulator(uncombinedVertices,
+            uncombinedFaces,
+            m_clothCollisionVertices,
+            m_clothCollisionTriangles);
+        clothSimulator.create();
+        for (size_t i = 0; i < 30; ++i)
+            clothSimulator.step();
+        clothSimulator.getCurrentVertices(&uncombinedVertices);
+    
+        auto vertexStartIndex = m_outcome->vertices.size();
+        auto updateVertexIndices = [=](std::vector<std::vector<size_t>> &faces) {
+            for (auto &it: faces) {
+                for (auto &subIt: it)
+                    subIt += vertexStartIndex;
+            }
+        };
+        updateVertexIndices(uncombinedFaces);
+        
+        m_outcome->vertices.insert(m_outcome->vertices.end(), uncombinedVertices.begin(), uncombinedVertices.end());
+        m_outcome->triangles.insert(m_outcome->triangles.end(), uncombinedFaces.begin(), uncombinedFaces.end());
+        m_outcome->triangleAndQuads.insert(m_outcome->triangleAndQuads.end(), uncombinedFaces.begin(), uncombinedFaces.end());
+        
+        m_outcome->nodes.insert(m_outcome->nodes.end(), componentCache.outcomeNodes.begin(), componentCache.outcomeNodes.end());
+        for (size_t i = 0; i < uncombinedVertices.size(); ++i) {
+            const auto &source = uncombinedVertexSources[i];
+            if (source.first.isNull())
+                continue;
+            m_outcome->nodeVertices.push_back(std::make_pair(uncombinedVertices[i], source));
+        }
+        return;
+    }
+    for (const auto &childIdString: valueOfKeyInMapOrEmpty(*component, "children").split(",")) {
+        if (childIdString.isEmpty())
+            continue;
+        collectClothComponent(childIdString);
+    }
+}
+
 void MeshGenerator::generateSmoothTriangleVertexNormals(const std::vector<QVector3D> &vertices, const std::vector<std::vector<size_t>> &triangles,
     const std::vector<QVector3D> &triangleNormals,
     std::vector<std::vector<QVector3D>> *triangleVertexNormals)

src/meshgenerator.h

@@ -11,6 +11,7 @@
 #include "snapshot.h"
 #include "combinemode.h"
 #include "meshloader.h"
+#include "componentlayer.h"
 
 class GeneratedPart
 {
@@ -100,6 +101,8 @@ private:
     std::map<QUuid, StrokeMeshBuilder::CutFaceTransform> *m_cutFaceTransforms = nullptr;
     std::map<QUuid, std::map<QString, QVector2D>> *m_nodesCutFaces = nullptr;
     quint64 m_id = 0;
+    std::vector<QVector3D> m_clothCollisionVertices;
+    std::vector<std::vector<size_t>> m_clothCollisionTriangles;
     
     void collectParts();
     bool checkIsComponentDirty(const QString &componentIdString);
@@ -125,7 +128,9 @@ private:
         GeneratedComponent &componentCache);
     MeshCombiner::Mesh *combineMultipleMeshes(const std::vector<std::tuple<MeshCombiner::Mesh *, CombineMode, QString>> &multipleMeshes, bool recombine=true);
     QString componentColorName(const std::map<QString, QString> *component);
+    ComponentLayer componentLayer(const std::map<QString, QString> *component);
     void collectUncombinedComponent(const QString &componentIdString);
+    void collectClothComponent(const QString &componentIdString);
     void cutFaceStringToCutTemplate(const QString &cutFaceString, std::vector<QVector2D> &cutTemplate);
     void remesh(const std::vector<OutcomeNode> &inputNodes,
         const std::vector<QVector3D> &inputVertices,

src/parttreewidget.cpp

@@ -315,6 +315,17 @@ void PartTreeWidget::showContextMenu(const QPoint &pos)
     }
     QWidget *widget = new QWidget;
     if (nullptr != component) {
+        QComboBox *componentLayerSelectBox = new QComboBox;
+        for (size_t i = 0; i < (size_t)ComponentLayer::Count; ++i) {
+            ComponentLayer layer = (ComponentLayer)i;
+            componentLayerSelectBox->addItem(ComponentLayerToDispName(layer));
+        }
+        componentLayerSelectBox->setCurrentIndex((int)component->layer);
+        connect(componentLayerSelectBox, static_cast<void (QComboBox::*)(int)>(&QComboBox::currentIndexChanged), this, [=](int index) {
+            emit setComponentLayer(component->id, (ComponentLayer)index);
+            emit groupOperationAdded();
+        });
+        
         QComboBox *combineModeSelectBox = new QComboBox;
         for (size_t i = 0; i < (size_t)CombineMode::Count; ++i) {
             CombineMode mode = (CombineMode)i;
@@ -368,6 +379,7 @@ void PartTreeWidget::showContextMenu(const QPoint &pos)
         if (nullptr != partTargetSelectBox)
             componentSettingsLayout->addRow(tr("Target"), partTargetSelectBox);
         componentSettingsLayout->addRow(tr("Mode"), combineModeSelectBox);
+        componentSettingsLayout->addRow(tr("Layer"), componentLayerSelectBox);
     
         QVBoxLayout *newLayout = new QVBoxLayout;
         newLayout->addLayout(layout);

src/parttreewidget.h

@@ -23,6 +23,7 @@ signals:
     void setComponentSmoothAll(QUuid componentId, float toSmoothAll);
     void setComponentSmoothSeam(QUuid componentId, float toSmoothSeam);
     void setComponentPolyCount(QUuid componentId, PolyCount count);
+    void setComponentLayer(QUuid componentId, ComponentLayer layer);
     void setPartTarget(QUuid partId, PartTarget target);
     void setPartBase(QUuid partId, PartBase base);
     void moveComponent(QUuid componentId, QUuid toParentId);

thirdparty/FastMassSpring/.gitattributes

@@ -0,0 +1,63 @@
+###############################################################################
+# Set default behavior to automatically normalize line endings.
+###############################################################################
+* text=auto
+
+###############################################################################
+# Set default behavior for command prompt diff.
+#
+# This is need for earlier builds of msysgit that does not have it on by
+# default for csharp files.
+# Note: This is only used by command line
+###############################################################################
+#*.cs     diff=csharp
+
+###############################################################################
+# Set the merge driver for project and solution files
+#
+# Merging from the command prompt will add diff markers to the files if there
+# are conflicts (Merging from VS is not affected by the settings below, in VS
+# the diff markers are never inserted). Diff markers may cause the following 
+# file extensions to fail to load in VS. An alternative would be to treat
+# these files as binary and thus will always conflict and require user
+# intervention with every merge. To do so, just uncomment the entries below
+###############################################################################
+#*.sln       merge=binary
+#*.csproj    merge=binary
+#*.vbproj    merge=binary
+#*.vcxproj   merge=binary
+#*.vcproj    merge=binary
+#*.dbproj    merge=binary
+#*.fsproj    merge=binary
+#*.lsproj    merge=binary
+#*.wixproj   merge=binary
+#*.modelproj merge=binary
+#*.sqlproj   merge=binary
+#*.wwaproj   merge=binary
+
+###############################################################################
+# behavior for image files
+#
+# image files are treated as binary by default.
+###############################################################################
+#*.jpg   binary
+#*.png   binary
+#*.gif   binary
+
+###############################################################################
+# diff behavior for common document formats
+# 
+# Convert binary document formats to text before diffing them. This feature
+# is only available from the command line. Turn it on by uncommenting the 
+# entries below.
+###############################################################################
+#*.doc   diff=astextplain
+#*.DOC   diff=astextplain
+#*.docx  diff=astextplain
+#*.DOCX  diff=astextplain
+#*.dot   diff=astextplain
+#*.DOT   diff=astextplain
+#*.pdf   diff=astextplain
+#*.PDF   diff=astextplain
+#*.rtf   diff=astextplain
+#*.RTF   diff=astextplain

thirdparty/FastMassSpring/.gitignore

@@ -0,0 +1,242 @@
+## Ignore Visual Studio temporary files, build results, and
+## files generated by popular Visual Studio add-ons.
+
+# User-specific files
+*.suo
+*.user
+*.userosscache
+*.sln.docstates
+
+# User-specific files (MonoDevelop/Xamarin Studio)
+*.userprefs
+
+# Build results
+[Dd]ebug/
+[Dd]ebugPublic/
+[Rr]elease/
+[Rr]eleases/
+[Xx]64/
+[Xx]86/
+[Bb]uild/
+bld/
+[Bb]in/
+[Oo]bj/
+
+# Visual Studio 2015 cache/options directory
+.vs/
+# Uncomment if you have tasks that create the project's static files in wwwroot
+#wwwroot/
+
+# MSTest test Results
+[Tt]est[Rr]esult*/
+[Bb]uild[Ll]og.*
+
+# NUNIT
+*.VisualState.xml
+TestResult.xml
+
+# Build Results of an ATL Project
+[Dd]ebugPS/
+[Rr]eleasePS/
+dlldata.c
+
+# DNX
+project.lock.json
+artifacts/
+
+*_i.c
+*_p.c
+*_i.h
+*.ilk
+*.meta
+*.obj
+*.pch
+*.pdb
+*.pgc
+*.pgd
+*.rsp
+*.sbr
+*.tlb
+*.tli
+*.tlh
+*.tmp
+*.tmp_proj
+*.log
+*.vspscc
+*.vssscc
+.builds
+*.pidb
+*.svclog
+*.scc
+
+# Chutzpah Test files
+_Chutzpah*
+
+# Visual C++ cache files
+ipch/
+*.aps
+*.ncb
+*.opendb
+*.opensdf
+*.sdf
+*.cachefile
+*.VC.db
+
+# Visual Studio profiler
+*.psess
+*.vsp
+*.vspx
+*.sap
+
+# TFS 2012 Local Workspace
+$tf/
+
+# Guidance Automation Toolkit
+*.gpState
+
+# ReSharper is a .NET coding add-in
+_ReSharper*/
+*.[Rr]e[Ss]harper
+*.DotSettings.user
+
+# JustCode is a .NET coding add-in
+.JustCode
+
+# TeamCity is a build add-in
+_TeamCity*
+
+# DotCover is a Code Coverage Tool
+*.dotCover
+
+# NCrunch
+_NCrunch_*
+.*crunch*.local.xml
+nCrunchTemp_*
+
+# MightyMoose
+*.mm.*
+AutoTest.Net/
+
+# Web workbench (sass)
+.sass-cache/
+
+# Installshield output folder
+[Ee]xpress/
+
+# DocProject is a documentation generator add-in
+DocProject/buildhelp/
+DocProject/Help/*.HxT
+DocProject/Help/*.HxC
+DocProject/Help/*.hhc
+DocProject/Help/*.hhk
+DocProject/Help/*.hhp
+DocProject/Help/Html2
+DocProject/Help/html
+
+# Click-Once directory
+publish/
+
+# Publish Web Output
+*.[Pp]ublish.xml
+*.azurePubxml
+
+# TODO: Un-comment the next line if you do not want to checkin
+# your web deploy settings because they may include unencrypted
+# passwords
+#*.pubxml
+*.publishproj
+
+# NuGet Packages
+*.nupkg
+# The packages folder can be ignored because of Package Restore
+**/packages/*
+# except build/, which is used as an MSBuild target.
+!**/packages/build/
+# Uncomment if necessary however generally it will be regenerated when needed
+#!**/packages/repositories.config
+# NuGet v3's project.json files produces more ignoreable files
+*.nuget.props
+*.nuget.targets
+
+# Microsoft Azure Build Output
+csx/
+*.build.csdef
+
+# Microsoft Azure Emulator
+ecf/
+rcf/
+
+# Windows Store app package directory
+AppPackages/
+BundleArtifacts/
+
+# Visual Studio cache files
+# files ending in .cache can be ignored
+*.[Cc]ache
+# but keep track of directories ending in .cache
+!*.[Cc]ache/
+
+# Others
+ClientBin/
+[Ss]tyle[Cc]op.*
+~$*
+*~
+*.dbmdl
+*.dbproj.schemaview
+*.pfx
+*.publishsettings
+node_modules/
+orleans.codegen.cs
+
+# RIA/Silverlight projects
+Generated_Code/
+
+# Backup & report files from converting an old project file
+# to a newer Visual Studio version. Backup files are not needed,
+# because we have git ;-)
+_UpgradeReport_Files/
+Backup*/
+UpgradeLog*.XML
+UpgradeLog*.htm
+
+# SQL Server files
+*.mdf
+*.ldf
+
+# Business Intelligence projects
+*.rdl.data
+*.bim.layout
+*.bim_*.settings
+
+# Microsoft Fakes
+FakesAssemblies/
+
+# GhostDoc plugin setting file
+*.GhostDoc.xml
+
+# Node.js Tools for Visual Studio
+.ntvs_analysis.dat
+
+# Visual Studio 6 build log
+*.plg
+
+# Visual Studio 6 workspace options file
+*.opt
+
+# Visual Studio LightSwitch build output
+**/*.HTMLClient/GeneratedArtifacts
+**/*.DesktopClient/GeneratedArtifacts
+**/*.DesktopClient/ModelManifest.xml
+**/*.Server/GeneratedArtifacts
+**/*.Server/ModelManifest.xml
+_Pvt_Extensions
+
+# LightSwitch generated files
+GeneratedArtifacts/
+ModelManifest.xml
+
+# Paket dependency manager
+.paket/paket.exe
+
+# FAKE - F# Make
+.fake/

thirdparty/FastMassSpring/ClothApp/ClothApp.vcxproj

@@ -0,0 +1,204 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" ToolsVersion="14.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations">
+    <ProjectConfiguration Include="Debug|Win32">
+      <Configuration>Debug</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|Win32">
+      <Configuration>Release</Configuration>
+      <Platform>Win32</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Debug|x64">
+      <Configuration>Debug</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+    <ProjectConfiguration Include="Release|x64">
+      <Configuration>Release</Configuration>
+      <Platform>x64</Platform>
+    </ProjectConfiguration>
+  </ItemGroup>
+  <PropertyGroup Label="Globals">
+    <ProjectGuid>{bbca6691-6c35-4bfe-9cb9-13f16d29990f}</ProjectGuid>
+    <Keyword>Win32Proj</Keyword>
+    <RootNamespace>ClothApp</RootNamespace>
+    <WindowsTargetPlatformVersion>8.1</WindowsTargetPlatformVersion>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v140</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v140</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>true</UseDebugLibraries>
+    <PlatformToolset>v140</PlatformToolset>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
+    <ConfigurationType>Application</ConfigurationType>
+    <UseDebugLibraries>false</UseDebugLibraries>
+    <PlatformToolset>v140</PlatformToolset>
+    <WholeProgramOptimization>true</WholeProgramOptimization>
+    <CharacterSet>Unicode</CharacterSet>
+  </PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <ImportGroup Label="ExtensionSettings">
+  </ImportGroup>
+  <ImportGroup Label="Shared">
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
+  </ImportGroup>
+  <PropertyGroup Label="UserMacros" />
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <LinkIncremental>true</LinkIncremental>
+    <IncludePath>$(SolutionDir)packages\nupengl.core.0.1.0.1\build\native\include;$(IncludePath)</IncludePath>
+    <Linkage-freeglut>static</Linkage-freeglut>
+    <Linkage-glew>static</Linkage-glew>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <LinkIncremental>true</LinkIncremental>
+    <IncludePath>C:\Program Files\Eigen\include;C:\Program Files\OpenMesh 7.1\include;$(IncludePath)</IncludePath>
+    <LibraryPath>C:\Program Files\OpenMesh 7.1\lib;$(LibraryPath)</LibraryPath>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <LinkIncremental>false</LinkIncremental>
+    <IncludePath>$(SolutionDir)packages\nupengl.core.0.1.0.1\build\native\include;$(IncludePath)</IncludePath>
+  </PropertyGroup>
+  <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <LinkIncremental>false</LinkIncremental>
+    <IncludePath>C:\Program Files\Eigen\include;C:\Program Files\OpenMesh 7.1\include;$(IncludePath)</IncludePath>
+    <LibraryPath>C:\Program Files\OpenMesh 7.1\lib;$(LibraryPath)</LibraryPath>
+  </PropertyGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
+    <ClCompile>
+      <PrecompiledHeader>
+      </PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <PreprocessorDefinitions>WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <SDLCheck>true</SDLCheck>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalLibraryDirectories>%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
+      <AdditionalDependencies>%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
+    <ClCompile>
+      <PrecompiledHeader>
+      </PrecompiledHeader>
+      <WarningLevel>Level3</WarningLevel>
+      <Optimization>Disabled</Optimization>
+      <PreprocessorDefinitions>_USE_MATH_DEFINES;_SCL_SECURE_NO_WARNINGS</PreprocessorDefinitions>
+      <SDLCheck>true</SDLCheck>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>OpenMeshCored.lib;OpenMeshToolsd.lib;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <PrecompiledHeader>
+      </PrecompiledHeader>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <PreprocessorDefinitions>WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+      <SDLCheck>true</SDLCheck>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
+    <ClCompile>
+      <WarningLevel>Level3</WarningLevel>
+      <PrecompiledHeader>
+      </PrecompiledHeader>
+      <Optimization>MaxSpeed</Optimization>
+      <FunctionLevelLinking>true</FunctionLevelLinking>
+      <IntrinsicFunctions>true</IntrinsicFunctions>
+      <PreprocessorDefinitions>_USE_MATH_DEFINES;_SCL_SECURE_NO_WARNINGS;NDEBUG</PreprocessorDefinitions>
+      <SDLCheck>true</SDLCheck>
+    </ClCompile>
+    <Link>
+      <SubSystem>Console</SubSystem>
+      <EnableCOMDATFolding>true</EnableCOMDATFolding>
+      <OptimizeReferences>true</OptimizeReferences>
+      <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalDependencies>OpenMeshCore.lib;OpenMeshTools.lib;%(AdditionalDependencies)</AdditionalDependencies>
+    </Link>
+  </ItemDefinitionGroup>
+  <ItemGroup>
+    <None Include="packages.config" />
+    <None Include="shaders\phong.fshader">
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild>
+    </None>
+    <None Include="shaders\pick.fshader">
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild>
+    </None>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="App.cpp" />
+    <ClCompile Include="MassSpringSolver.cpp" />
+    <ClCompile Include="Mesh.cpp" />
+    <ClCompile Include="Renderer.cpp" />
+    <ClCompile Include="Shader.cpp" />
+    <ClCompile Include="UserInteraction.cpp" />
+  </ItemGroup>
+  <ItemGroup>
+    <FxCompile Include="shaders\basic.vshader">
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild>
+    </FxCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <ClInclude Include="MassSpringSolver.h" />
+    <ClInclude Include="Mesh.h" />
+    <ClInclude Include="Renderer.h" />
+    <ClInclude Include="Shader.h" />
+    <ClInclude Include="UserInteraction.h" />
+  </ItemGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <ImportGroup Label="ExtensionTargets">
+    <Import Project="..\packages\nupengl.core.redist.0.1.0.1\build\native\nupengl.core.redist.targets" Condition="Exists('..\packages\nupengl.core.redist.0.1.0.1\build\native\nupengl.core.redist.targets')" />
+    <Import Project="..\packages\nupengl.core.0.1.0.1\build\native\nupengl.core.targets" Condition="Exists('..\packages\nupengl.core.0.1.0.1\build\native\nupengl.core.targets')" />
+  </ImportGroup>
+  <Target Name="EnsureNuGetPackageBuildImports" BeforeTargets="PrepareForBuild">
+    <PropertyGroup>
+      <ErrorText>This project references NuGet package(s) that are missing on this computer. Use NuGet Package Restore to download them.  For more information, see http://go.microsoft.com/fwlink/?LinkID=322105. The missing file is {0}.</ErrorText>
+    </PropertyGroup>
+    <Error Condition="!Exists('..\packages\nupengl.core.redist.0.1.0.1\build\native\nupengl.core.redist.targets')" Text="$([System.String]::Format('$(ErrorText)', '..\packages\nupengl.core.redist.0.1.0.1\build\native\nupengl.core.redist.targets'))" />
+    <Error Condition="!Exists('..\packages\nupengl.core.0.1.0.1\build\native\nupengl.core.targets')" Text="$([System.String]::Format('$(ErrorText)', '..\packages\nupengl.core.0.1.0.1\build\native\nupengl.core.targets'))" />
+  </Target>
+</Project>

thirdparty/FastMassSpring/ClothApp/ClothApp.vcxproj.filters

@@ -0,0 +1,71 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup>
+    <Filter Include="Source Files">
+      <UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
+      <Extensions>cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
+    </Filter>
+    <Filter Include="Header Files">
+      <UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
+      <Extensions>h;hh;hpp;hxx;hm;inl;inc;xsd</Extensions>
+    </Filter>
+    <Filter Include="Resource Files">
+      <UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
+      <Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
+    </Filter>
+    <Filter Include="Shader Files">
+      <UniqueIdentifier>{55e50294-8ec8-4842-bc08-005ce41662d4}</UniqueIdentifier>
+    </Filter>
+  </ItemGroup>
+  <ItemGroup>
+    <None Include="packages.config" />
+    <None Include="shaders\phong.fshader">
+      <Filter>Shader Files</Filter>
+    </None>
+    <None Include="shaders\pick.fshader">
+      <Filter>Shader Files</Filter>
+    </None>
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="App.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="Renderer.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="Mesh.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="MassSpringSolver.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="UserInteraction.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="Shader.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <FxCompile Include="shaders\basic.vshader">
+      <Filter>Shader Files</Filter>
+    </FxCompile>
+  </ItemGroup>
+  <ItemGroup>
+    <ClInclude Include="Mesh.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="Renderer.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="MassSpringSolver.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="UserInteraction.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="Shader.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+  </ItemGroup>
+</Project>

thirdparty/FastMassSpring/ClothApp/MassSpringSolver.cpp

@@ -0,0 +1,384 @@
+#include "MassSpringSolver.h"
+#include <iostream>
+
+// S Y S T E M //////////////////////////////////////////////////////////////////////////////////////
+mass_spring_system::mass_spring_system(
+	unsigned int n_points,
+	unsigned int n_springs,
+	float time_step,
+	EdgeList spring_list,
+	VectorXf rest_lengths,
+	VectorXf stiffnesses,
+	VectorXf masses,
+	VectorXf fext,
+	float damping_factor
+)
+	: n_points(n_points), n_springs(n_springs),
+	time_step(time_step), spring_list(spring_list),
+	rest_lengths(rest_lengths), stiffnesses(stiffnesses), masses(masses),
+	fext(fext), damping_factor(damping_factor) {}
+
+// S O L V E R //////////////////////////////////////////////////////////////////////////////////////
+MassSpringSolver::MassSpringSolver(mass_spring_system* system, float* vbuff) 
+	: system(system), current_state(vbuff, system->n_points * 3), 
+	prev_state(current_state), spring_directions(system->n_springs * 3) {
+	
+	float h2 = system->time_step * system->time_step; // shorthand
+
+	// compute M, L, J
+	TripletList LTriplets, JTriplets;
+	// L
+	L.resize(3 * system->n_points, 3 * system->n_points);
+	unsigned int k = 0; // spring counter
+	for (Edge& i : system->spring_list) {
+		for (int j = 0; j < 3; j++) {
+			LTriplets.push_back(
+				Triplet(3 * i.first + j, 3 * i.first  + j,  1 * system->stiffnesses[k]));
+			LTriplets.push_back(
+				Triplet(3 * i.first + j, 3 * i.second + j, -1 * system->stiffnesses[k]));
+			LTriplets.push_back(
+				Triplet(3 * i.second + j, 3 * i.first + j, -1 * system->stiffnesses[k]));
+			LTriplets.push_back(
+				Triplet(3 * i.second + j, 3 * i.second + j, 1 * system->stiffnesses[k]));
+		}
+		k++;
+	}
+	L.setFromTriplets(LTriplets.begin(), LTriplets.end());
+
+	// J
+	J.resize(3 * system->n_points, 3 * system->n_springs);
+	k = 0; // spring counter
+	for (Edge& i : system->spring_list) {
+		for (unsigned int j = 0; j < 3; j++) {
+			JTriplets.push_back(
+				Triplet(3 * i.first  + j, 3 * k + j,  1 * system->stiffnesses[k]));
+			JTriplets.push_back(
+				Triplet(3 * i.second + j, 3 * k + j, -1 * system->stiffnesses[k]));
+		}
+		k++;
+	}
+	J.setFromTriplets(JTriplets.begin(), JTriplets.end());
+
+	// M
+	TripletList MTriplets;
+	M.resize(3 * system->n_points, 3 * system->n_points);
+	for (unsigned int i = 0; i < system->n_points; i++) {
+		for (int j = 0; j < 3; j++) {
+			MTriplets.push_back(Triplet(3 * i + j, 3 * i + j, system->masses[i]));
+		}
+	}
+	M.setFromTriplets(MTriplets.begin(), MTriplets.end());
+	
+	// pre-factor system matrix
+	SparseMatrix A = M + h2 * L;
+	system_matrix.compute(A);
+}
+
+void MassSpringSolver::globalStep() {
+	float h2 = system->time_step * system->time_step; // shorthand
+
+	// compute right hand side
+	VectorXf b = inertial_term
+		+ h2 * J * spring_directions
+		+ h2 * system->fext;
+
+	// solve system and update state
+	current_state = system_matrix.solve(b);
+}
+
+void MassSpringSolver::localStep() {
+	unsigned int j = 0;
+	for (Edge& i : system->spring_list) {
+		Vector3f p12(
+			current_state[3 * i.first + 0] - current_state[3 * i.second + 0],
+			current_state[3 * i.first + 1] - current_state[3 * i.second + 1],
+			current_state[3 * i.first + 2] - current_state[3 * i.second + 2]
+		);
+
+		p12.normalize();
+		spring_directions[3 * j + 0] = 	system->rest_lengths[j] * p12[0];
+		spring_directions[3 * j + 1] =	system->rest_lengths[j] * p12[1];
+		spring_directions[3 * j + 2] =	system->rest_lengths[j] * p12[2];
+		j++;
+	}
+}
+
+void MassSpringSolver::solve(unsigned int n) {
+	float a = system->damping_factor; // shorthand
+
+	// update inertial term
+	inertial_term = M * ((a + 1) * (current_state) - a * prev_state);
+
+	// save current state in previous state
+	prev_state = current_state;
+
+	// perform steps
+	for (unsigned int i = 0; i < n; i++) {
+		localStep();
+		globalStep();
+	}
+}
+
+void MassSpringSolver::timedSolve(unsigned int ms) {
+	// TODO
+}
+
+
+// B U I L D E R ////////////////////////////////////////////////////////////////////////////////////
+void MassSpringBuilder::uniformGrid(
+	unsigned int n,
+	float time_step,
+	float rest_length,
+	float stiffness,
+	float mass,
+	float damping_factor,
+	float gravity
+) {
+	// n must be odd
+	assert(n % 2 == 1);
+
+	// shorthand
+	const double root2 = 1.41421356237;
+
+	// compute n_points and n_springs
+	unsigned int n_points = n * n;
+	unsigned int n_springs = (n - 1) * (5 * n - 2);
+
+	// build mass list
+	VectorXf masses(mass * VectorXf::Ones(n_springs));
+
+	// build spring list and spring parameters
+	EdgeList spring_list(n_springs);
+	structI.reserve(2 * (n - 1) * n);
+	shearI.reserve(2 * (n - 1) * (n - 1));
+	bendI.reserve(n * (n - 1));
+
+	VectorXf rest_lengths(n_springs);
+	VectorXf stiffnesses(n_springs);
+	unsigned int k = 0; // spring counter
+	for(unsigned int i = 0; i < n; i++) {
+		for(unsigned int j = 0; j < n; j++) {
+			// bottom right corner
+			if(i == n - 1 && j == n - 1) {
+				continue;
+			}
+
+			if (i == n - 1) {
+				// structural spring
+				spring_list[k] = Edge(n * i + j, n * i + j + 1);
+				rest_lengths[k] = rest_length;
+				stiffnesses[k] = stiffness;
+				structI.push_back(k++);
+
+				// bending spring
+				if (j % 2 == 0) {
+					spring_list[k] = Edge(n * i + j, n * i + j + 2);
+					rest_lengths[k] = 2 * rest_length;
+					stiffnesses[k] = stiffness;
+					bendI.push_back(k++);
+				}
+				continue;
+			}
+
+			// right edge
+			if (j == n - 1) {
+				// structural spring
+				spring_list[k] = Edge(n * i + j, n * (i + 1) + j);
+				rest_lengths[k] = rest_length;
+				stiffnesses[k] = stiffness;
+				structI.push_back(k++);
+
+				// bending spring
+				if (i % 2 == 0){
+					spring_list[k] = Edge(n * i + j, n * (i + 2) + j);
+					rest_lengths[k] = 2 * rest_length;
+					stiffnesses[k] = stiffness;
+					bendI.push_back(k++);
+				}
+				continue;
+			}
+
+			// structural springs
+			spring_list[k] = Edge(n * i + j, n * i + j + 1);
+			rest_lengths[k] = rest_length;
+			stiffnesses[k] = stiffness;
+			structI.push_back(k++);
+
+			spring_list[k] = Edge(n * i + j, n * (i + 1) + j);
+			rest_lengths[k] = rest_length;
+			stiffnesses[k] = stiffness;
+			structI.push_back(k++);
+
+			// shearing springs
+			spring_list[k] = Edge(n * i + j, n * (i + 1) + j + 1);
+			rest_lengths[k] = root2 * rest_length;
+			stiffnesses[k] = stiffness;
+			shearI.push_back(k++);
+
+			spring_list[k] = Edge(n * (i + 1) + j, n * i + j + 1);
+			rest_lengths[k] = root2 * rest_length;
+			stiffnesses[k] = stiffness;
+			shearI.push_back(k++);
+
+			// bending springs
+			if (j % 2 == 0) {
+				spring_list[k] = Edge(n * i + j, n * i + j + 2);
+				rest_lengths[k] = 2 * rest_length;
+				stiffnesses[k] = stiffness;
+				bendI.push_back(k++);
+			}
+			if (i % 2 == 0) {
+				spring_list[k] = Edge(n * i + j, n * (i + 2) + j);
+				rest_lengths[k] = 2 * rest_length;
+				stiffnesses[k] = stiffness;
+				bendI.push_back(k++);
+			}
+		}
+	}
+
+	// compute external forces
+	VectorXf fext = Vector3f(0, 0, -gravity).replicate(n_points, 1);
+
+	result = new mass_spring_system(n_points, n_springs, time_step, spring_list, rest_lengths,
+		stiffnesses, masses, fext, damping_factor);
+}
+MassSpringBuilder::IndexList MassSpringBuilder::getStructIndex() { return structI; }
+MassSpringBuilder::IndexList MassSpringBuilder::getShearIndex() { return shearI; }
+MassSpringBuilder::IndexList MassSpringBuilder::getBendIndex() { return bendI; }
+mass_spring_system* MassSpringBuilder::getResult() { return result; }
+
+// C O N S T R A I N T //////////////////////////////////////////////////////////////////////////////
+CgNode::CgNode(mass_spring_system* system, float* vbuff) : system(system), vbuff(vbuff) {}
+
+// point node
+CgPointNode::CgPointNode(mass_spring_system* system, float* vbuff) : CgNode(system, vbuff) {}
+bool CgPointNode::accept(CgNodeVisitor& visitor) { return visitor.visit(*this); }
+
+// spring node
+CgSpringNode::CgSpringNode(mass_spring_system* system, float* vbuff) : CgNode(system, vbuff) {}
+bool CgSpringNode::accept(CgNodeVisitor& visitor) {
+	for (CgNode* child : children) {
+		if (!child->accept(visitor)) return false;
+	}
+	return visitor.visit(*this);
+}
+void CgSpringNode::addChild(CgNode* node) { children.push_back(node); }
+void CgSpringNode::removeChild(CgNode* node) { 
+	children.erase(find(children.begin(), children.end(), node)); 
+}
+
+// root node
+CgRootNode::CgRootNode(mass_spring_system* system, float* vbuff) : CgSpringNode(system, vbuff) {}
+void CgRootNode::satisfy() { return; }
+bool CgRootNode::accept(CgNodeVisitor& visitor) {
+	for (CgNode* child : children) {
+		if (!child->accept(visitor)) return false;
+	}
+	return true;
+}
+
+// point fix node
+CgPointFixNode::CgPointFixNode(mass_spring_system* system, float* vbuff) 
+	: CgPointNode(system, vbuff) {}
+bool CgPointFixNode::query(unsigned int i) const { return fix_map.find(3 * i) != fix_map.end(); }
+void CgPointFixNode::satisfy() {
+	for (auto fix : fix_map)
+		for (int i = 0; i < 3; i++)
+			vbuff[fix.first + i] = fix.second[i];
+}
+void CgPointFixNode::fixPoint(unsigned int i) {
+	assert(i >= 0 && i < system->n_points);
+	fix_map[3 * i] = Vector3f(vbuff[3 * i], vbuff[3 * i + 1], vbuff[3 * i + 2]);
+}
+void CgPointFixNode::releasePoint(unsigned int i) { fix_map.erase(3 * i); }
+
+// spring deformation node
+CgSpringDeformationNode::CgSpringDeformationNode(mass_spring_system* system, float* vbuff,
+	float tauc, unsigned int n_iter) : CgSpringNode(system, vbuff), tauc(tauc), n_iter(n_iter) {}
+void CgSpringDeformationNode::satisfy() {
+	for (int k = 0; k < n_iter; k++) {
+		for (unsigned int i : items) {
+			Edge spring = system->spring_list[i];
+			CgQueryFixedPointVisitor visitor;
+
+			Vector3f p12(
+				vbuff[3 * spring.first + 0] - vbuff[3 * spring.second + 0],
+				vbuff[3 * spring.first + 1] - vbuff[3 * spring.second + 1],
+				vbuff[3 * spring.first + 2] - vbuff[3 * spring.second + 2]
+			);
+
+			float len = p12.norm();
+			float rlen = system->rest_lengths[i];
+			float diff = (len - (1 + tauc) * rlen) / len;
+			float rate = (len - rlen) / rlen;
+
+			// check deformation
+			if (rate <= tauc) continue;
+
+			// check if points are fixed
+			float f1, f2;
+			f1 = f2 = 0.5f;
+
+			// if first point is fixed
+			if (visitor.queryPoint(*this, spring.first)) { f1 = 0.0f; f2 = 1.0f; }
+
+			// if second point is fixed
+			if (visitor.queryPoint(*this, spring.second)) {
+				f1 = (f1 != 0.0f ? 1.0f : 0.0f);
+				f2 = 0.0f;
+			}
+
+			for (int j = 0; j < 3; j++) {
+				vbuff[3 * spring.first + j] -= p12[j] * f1 * diff;
+				vbuff[3 * spring.second + j] += p12[j] * f2 * diff;
+			}
+		}
+	}
+}
+void CgSpringDeformationNode::addSprings(std::vector<unsigned int> springs) {
+	items.insert(springs.begin(), springs.end());
+}
+
+// sphere collision node
+CgSphereCollisionNode::CgSphereCollisionNode(mass_spring_system* system, float* vbuff,
+	float radius, Vector3f center) : CgPointNode(system, vbuff), radius(radius), center(center) {}
+bool CgSphereCollisionNode::query(unsigned int i) const { return false; }
+void CgSphereCollisionNode::satisfy() {
+	for (int i = 0; i < system->n_points; i++) {
+		Vector3f p(
+			vbuff[3 * i + 0] - center[0],
+			vbuff[3 * i + 1] - center[1],
+			vbuff[3 * i + 2] - center[2]
+		);
+
+		if (p.norm() < radius) {
+			p.normalize();
+			p = radius * p;
+		}
+		else continue;
+
+		for (int j = 0; j < 3; j++) {
+			vbuff[3 * i + j] = p[j] + center[j];
+		}
+	}
+}
+
+// node visitor
+bool CgNodeVisitor::visit(CgPointNode& node) { return true; }
+bool CgNodeVisitor::visit(CgSpringNode& node) { return true; }
+
+// query fixed point visitor
+bool CgQueryFixedPointVisitor::visit(CgPointNode& node) {
+	queryResult = node.query(i);
+	return !queryResult;
+}
+bool CgQueryFixedPointVisitor::queryPoint(CgNode& root, unsigned int i) {
+	this->i = i;
+	root.accept(*this);
+	return queryResult;
+}
+
+// satisfy visitor
+bool CgSatisfyVisitor::visit(CgPointNode& node) { node.satisfy(); return true; }
+bool CgSatisfyVisitor::visit(CgSpringNode& node) { node.satisfy(); return true; }
+void CgSatisfyVisitor::satisfy(CgNode& root) { root.accept(*this); }

thirdparty/FastMassSpring/ClothApp/MassSpringSolver.h

@@ -0,0 +1,230 @@
+#pragma once
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <vector>
+#include <unordered_map>
+#include <unordered_set>
+
+// Mass-Spring System struct
+struct mass_spring_system { 
+	typedef Eigen::SparseMatrix<float> SparseMatrix;
+	typedef Eigen::VectorXf VectorXf;
+	typedef std::pair<unsigned int, unsigned int> Edge;
+	typedef std::vector<Edge> EdgeList;
+	
+	// parameters
+	unsigned int n_points; // number of points
+	unsigned int n_springs; // number of springs
+	float time_step; // time step
+	EdgeList spring_list; // spring edge list
+	VectorXf rest_lengths; // spring rest lengths
+	VectorXf stiffnesses; // spring stiffnesses
+	VectorXf masses; // point masses
+	VectorXf fext; // external forces
+	float damping_factor; // damping factor
+	
+	mass_spring_system(
+		unsigned int n_points,       // number of points
+		unsigned int n_springs,      // number of springs
+		float time_step,             // time step
+		EdgeList spring_list,        // spring edge list
+		VectorXf rest_lengths,       // spring rest lengths
+		VectorXf stiffnesses,        // spring stiffnesses
+		VectorXf masses,             // point masses
+		VectorXf fext,               // external forces
+		float damping_factor         // damping factor
+	);
+};
+
+// Mass-Spring System Solver class
+class MassSpringSolver {
+private:
+	typedef Eigen::Vector3f Vector3f;
+	typedef Eigen::VectorXf VectorXf;
+	typedef Eigen::SparseMatrix<float> SparseMatrix;
+	typedef Eigen::SimplicialLLT<Eigen::SparseMatrix<float> > Cholesky;
+	typedef Eigen::Map<Eigen::VectorXf> Map;
+	typedef std::pair<unsigned int, unsigned int> Edge;
+	typedef Eigen::Triplet<float> Triplet;
+	typedef std::vector<Triplet> TripletList;
+
+	// system
+	mass_spring_system* system;
+	Cholesky system_matrix;
+
+	// M, L, J matrices
+	SparseMatrix M;
+	SparseMatrix L;
+	SparseMatrix J;
+
+	// state
+	Map current_state; // q(n), current state
+	VectorXf prev_state; // q(n - 1), previous state
+	VectorXf spring_directions; // d, spring directions
+	VectorXf inertial_term; // M * y, y = (a + 1) * q(n) - a * q(n - 1)
+
+	// steps
+	void globalStep();
+	void localStep();
+
+public:
+	MassSpringSolver(mass_spring_system* system, float* vbuff);
+
+	// solve iterations
+	void solve(unsigned int n);
+	void timedSolve(unsigned int ms);
+};
+
+// Mass-Spring System Builder Class
+class MassSpringBuilder {
+private:
+	typedef Eigen::Vector3f Vector3f;
+	typedef Eigen::VectorXf VectorXf;	
+	typedef std::pair<unsigned int, unsigned int> Edge;
+	typedef std::vector<Edge> EdgeList;
+	typedef Eigen::Triplet<float> Triplet;
+	typedef std::vector<Triplet> TripletList;
+	typedef std::vector<unsigned int> IndexList;
+
+	IndexList structI, shearI, bendI;
+	mass_spring_system* result;
+
+public:
+	void uniformGrid(
+		unsigned int n,          // grid width
+		float time_step,         // time step
+		float rest_length,       // spring rest length (non-diagonal)
+		float stiffness,         // spring stiffness
+		float mass,              // node mass
+		float damping_factor,    // damping factor
+		float gravity            // gravitationl force (-z axis)
+	);
+
+
+	// indices
+	IndexList getStructIndex(); // structural springs
+	IndexList getShearIndex(); // shearing springs
+	IndexList getBendIndex(); // bending springs
+
+	mass_spring_system* getResult();
+};
+
+// Constraint Graph
+class CgNodeVisitor; // Constraint graph node visitor
+
+// Constraint graph node
+class CgNode {
+protected:
+	mass_spring_system* system;
+	float* vbuff;
+
+public:
+	CgNode(mass_spring_system* system, float* vbuff);
+
+	virtual void satisfy() = 0; // satisfy constraint
+	virtual bool accept(CgNodeVisitor& visitor) = 0; // accept visitor
+
+};
+
+// point constraint node
+class CgPointNode : public CgNode {
+public:
+	CgPointNode(mass_spring_system* system, float* vbuff);
+	virtual bool query(unsigned int i) const = 0; // check if point with index i is constrained
+	virtual bool accept(CgNodeVisitor& visitor);
+
+};
+
+// spring constraint node
+class CgSpringNode : public CgNode {
+protected:
+	typedef std::vector<CgNode*> NodeList;
+	NodeList children;
+
+public:
+	CgSpringNode(mass_spring_system* system, float* vbuff);
+
+	virtual bool accept(CgNodeVisitor& visitor);
+	void addChild(CgNode* node);
+	void removeChild(CgNode* node);
+};
+
+// root node 
+class CgRootNode : public CgSpringNode {
+public:
+	CgRootNode(mass_spring_system* system, float* vbuff);
+
+	virtual void satisfy();
+	virtual bool accept(CgNodeVisitor& visitor);
+};
+
+class CgPointFixNode : public CgPointNode {
+protected:
+	typedef Eigen::Vector3f Vector3f;
+	std::unordered_map<unsigned int, Vector3f> fix_map;
+public:
+	CgPointFixNode(mass_spring_system* system, float* vbuff);
+	virtual void satisfy();
+
+	virtual bool query(unsigned int i) const;
+	virtual void fixPoint(unsigned int i); // add point at index i to list
+	virtual void releasePoint(unsigned int i); // remove point at index i from list
+};
+
+// spring node
+class CgSpringDeformationNode : public CgSpringNode {
+private:
+	typedef std::pair<unsigned int, unsigned int> Edge;
+	typedef Eigen::Vector3f Vector3f;
+	std::unordered_set<unsigned int> items;
+	float tauc; // critical deformation rate
+	unsigned int n_iter; // number of iterations
+
+public:
+	CgSpringDeformationNode(mass_spring_system* system, float* vbuff, float tauc, unsigned int n_iter);
+	virtual void satisfy();
+
+	void addSprings(std::vector<unsigned int> springs);
+};
+
+// sphere collision node
+class CgSphereCollisionNode : public CgPointNode {
+private:
+	typedef Eigen::Vector3f Vector3f;
+
+	float radius;
+	Vector3f center;
+
+public:
+	CgSphereCollisionNode(mass_spring_system* system, float* vbuff, float radius, Vector3f center);
+	virtual bool query(unsigned int i) const;
+	virtual void satisfy();
+};
+
+// node visitor
+class CgNodeVisitor {
+public:
+	
+	virtual bool visit(CgPointNode& node);
+	virtual bool visit(CgSpringNode& node);
+};
+
+// fixed point query visitor
+class CgQueryFixedPointVisitor : public CgNodeVisitor {
+private:
+	unsigned int i;
+	bool queryResult;
+public:
+	virtual bool visit(CgPointNode& node);
+
+	bool queryPoint(CgNode& root, unsigned int i);
+};
+
+// satisfy visitor
+class CgSatisfyVisitor : public CgNodeVisitor {
+public:
+	virtual bool visit(CgPointNode& node);
+	virtual bool visit(CgSpringNode& node);
+
+	void satisfy(CgNode& root);
+};

thirdparty/FastMassSpring/ClothApp/Mesh.cpp

@@ -0,0 +1,76 @@
+#include "Mesh.h"
+
+// M E S H /////////////////////////////////////////////////////////////////////////////////////
+float* Mesh::vbuff() { return VERTEX_DATA(this); }
+float* Mesh::nbuff() { return NORMAL_DATA(this); }
+float* Mesh::tbuff() { return TEXTURE_DATA(this); }
+unsigned int* Mesh::ibuff() { return &_ibuff[0]; }
+void Mesh::useIBuff(std::vector<unsigned int>& _ibuff) { this->_ibuff = _ibuff; }
+
+unsigned int Mesh::vbuffLen() { return (unsigned int)n_vertices() * 3; }
+unsigned int Mesh::nbuffLen() { return (unsigned int)n_vertices() * 3; }
+unsigned int Mesh::tbuffLen() { return (unsigned int)n_vertices() * 2; }
+unsigned int Mesh::ibuffLen() { return (unsigned int)_ibuff.size(); }
+
+
+// M E S H  B U I L D E R /////////////////////////////////////////////////////////////////////
+void MeshBuilder::uniformGrid(float w, int n) {
+	result = new Mesh;
+	unsigned int ibuffLen = 6 * (n - 1) * (n - 1);
+	std::vector<unsigned int> ibuff(ibuffLen);
+
+	// request mesh properties
+	result->request_vertex_normals();
+	result->request_vertex_normals();
+	result->request_vertex_texcoords2D();
+
+	// generate mesh
+	unsigned int idx = 0; // vertex index
+	const float d = w / (n - 1); // step distance
+	const float ud = 1.0f / (n - 1); // unit step distance
+	const OpenMesh::Vec3f o = OpenMesh::Vec3f(-w/2.0f, w/2.0f, 0.0f); // origin
+	const OpenMesh::Vec3f ux = OpenMesh::Vec3f(1.0f, 0.0f, 0.0f); // unit x direction
+	const OpenMesh::Vec3f uy = OpenMesh::Vec3f(0.0f, -1.0f, 0.0f); // unit y direction
+	std::vector<OpenMesh::VertexHandle> handle_table(n * n); // table storing vertex handles for easy grid connectivity establishment
+
+	for (int i = 0; i < n; i++) {
+		for (int j = 0; j < n; j++) {
+			handle_table[j + i * n] = result->add_vertex(o + d*j*ux + d*i*uy); // add vertex
+			result->set_texcoord2D(handle_table[j + i * n], OpenMesh::Vec2f(ud*j, ud*i)); // add texture coordinates
+
+			//add connectivity
+			if (i > 0 && j < n - 1) {
+				result->add_face(
+					handle_table[j + i * n],
+					handle_table[j + 1 + (i - 1) * n],
+					handle_table[j + (i - 1) * n]
+				);
+
+				ibuff[idx++] = j + i * n;
+				ibuff[idx++] = j + 1 + (i - 1) * n;
+				ibuff[idx++] = j + (i - 1) * n;
+			}
+
+			if (j > 0 && i > 0) {
+				result->add_face(
+					handle_table[j + i * n],
+					handle_table[j + (i - 1) * n],
+					handle_table[j - 1 + i * n]
+				);
+
+				ibuff[idx++] = j + i * n;
+				ibuff[idx++] = j + (i - 1) * n;
+				ibuff[idx++] = j - 1 + i * n;
+			}
+		}
+	}
+
+	// calculate normals
+	result->request_face_normals();
+	result->update_normals();
+	result->release_face_normals();
+
+	// set index buffer
+	result->useIBuff(ibuff);
+}
+Mesh* MeshBuilder::getResult() { return result; }

thirdparty/FastMassSpring/ClothApp/Mesh.h

@@ -0,0 +1,42 @@
+#pragma once
+#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
+#include <vector>
+
+// Mesh type
+typedef OpenMesh::TriMesh_ArrayKernelT<> _Mesh;
+
+// Macros for extracting buffers from OpenMesh
+#define VERTEX_DATA(mesh) (float*) &(mesh->point(*mesh->vertices_begin()))
+#define NORMAL_DATA(mesh) (float*) &(mesh->normal(*mesh->vertices_begin()))
+#define TEXTURE_DATA(mesh) (float*) &(mesh->texcoord2D(*mesh->vertices_begin()))
+
+// Mesh class
+class Mesh : public _Mesh {
+private:
+	std::vector<unsigned int> _ibuff;
+
+public:
+	// pointers to buffers
+	float* vbuff();
+	float* nbuff();
+	float* tbuff();
+	unsigned int* ibuff();
+
+	// buffer sizes
+	unsigned int vbuffLen();
+	unsigned int nbuffLen();
+	unsigned int tbuffLen();
+	unsigned int ibuffLen();
+
+	// set index buffer
+	void useIBuff(std::vector<unsigned int>& _ibuff);
+};
+
+class MeshBuilder {
+private:
+	Mesh* result;
+
+public:
+	void uniformGrid(float w, int n);
+	Mesh* getResult();
+};

thirdparty/FastMassSpring/ClothApp/Renderer.cpp

@@ -0,0 +1,49 @@
+#include "Renderer.h"
+#include <glm/gtc/type_ptr.hpp>
+
+Renderer::Renderer() {}
+
+void Renderer::setProgram(GLProgram* program) {
+	assert(program != nullptr);
+	assert((*program) > 0);
+	this->program = program;
+}
+
+void Renderer::setProgramInput(ProgramInput* input) {
+	this->input = input;
+}
+
+void Renderer::setModelview(const glm::mat4& mv) {
+	assert(program != nullptr);
+	assert((*program) > 0);
+	glUseProgram(*program);
+	program->setModelView(mv);
+	glUseProgram(0);
+}
+
+void Renderer::setProjection(const glm::mat4& p) {
+	assert(program != nullptr);
+	assert((*program) > 0);
+	glUseProgram(*program);
+	program->setProjection(p);
+	glUseProgram(0);
+}
+
+void Renderer::setElementCount(unsigned int n_elements) { this->n_elements = n_elements; }
+
+void Renderer::draw() {
+	assert(program != nullptr);
+	assert((*program) > 0);
+	glUseProgram(*program);
+	glBindVertexArray(*input);
+	glEnableVertexAttribArray(0);
+	glEnableVertexAttribArray(1);
+	glEnableVertexAttribArray(2);
+	glDrawElements(GL_TRIANGLES, n_elements, GL_UNSIGNED_INT, 0);
+	glDisableVertexAttribArray(0);
+	glDisableVertexAttribArray(1);
+	glDisableVertexAttribArray(2);
+	glBindVertexArray(0);
+	glUseProgram(0);
+}
+

thirdparty/FastMassSpring/ClothApp/Renderer.h

@@ -0,0 +1,22 @@
+#pragma once
+#include <glm/gtc/matrix_transform.hpp>
+#include "Shader.h"
+
+
+class Renderer {
+protected:
+	GLProgram* program;
+	ProgramInput* input;
+	unsigned int n_elements;
+
+public:
+	Renderer();
+
+	void setProgram(GLProgram* program);
+	void setProgramInput(ProgramInput* input);
+	void setModelview(const glm::mat4& mv);
+	void setProjection(const glm::mat4& p);
+	void setElementCount(unsigned int n_elements);
+
+	void draw();
+};

thirdparty/FastMassSpring/ClothApp/Shader.cpp

@@ -0,0 +1,190 @@
+#include "Shader.h"
+#include <vector>
+#include <glm/gtc/type_ptr.hpp>
+#include <iostream>
+
+// GLSHADER ///////////////////////////////////////////////////////////////////////////////////
+GLShader::GLShader(GLenum shaderType) : handle(glCreateShader(shaderType)) {};
+
+GLShader::operator GLuint() const {
+	return handle;
+}
+
+GLShader::~GLShader() {
+	glDeleteShader(handle);
+}
+
+void GLShader::compile(const char* source) {
+	GLint compiled = 0;  // Compiled flag
+	const char *ptrs[] = { source };
+	const GLint lens[] = { std::strlen(source) };
+	glShaderSource(handle, 1, ptrs, lens);
+	glCompileShader(handle);
+	glGetShaderiv(handle, GL_COMPILE_STATUS, &compiled);
+	if (!compiled) {
+		GLint logSize = 0;
+		glGetShaderiv(handle, GL_INFO_LOG_LENGTH, &logSize);
+		std::vector<GLchar> errorLog(logSize);
+		glGetShaderInfoLog(handle, logSize, &logSize, &errorLog[0]);
+		std::cerr << &errorLog[0] << std::endl;
+		throw std::runtime_error("Failed to compile shader.");
+	}
+}
+
+void GLShader::compile(std::ifstream& source) {
+	std::vector<char> text;
+	source.seekg(0, std::ios_base::end);
+	std::streampos fileSize = source.tellg();
+	text.resize(fileSize);
+
+	source.seekg(0, std::ios_base::beg);
+	source.read(&text[0], fileSize);
+	compile(&text[0]);
+}
+
+// GLPROGRAM //////////////////////////////////////////////////////////////////////////////////
+GLProgram::GLProgram() : handle(glCreateProgram()) {}
+
+void GLProgram::link(const GLShader& vshader, const GLShader& fshader) {
+	GLint linked = 0; // Linked flag
+	glAttachShader(handle, vshader);
+	glAttachShader(handle, fshader);
+	glLinkProgram(handle);
+	glDetachShader(handle, vshader);
+	glDetachShader(handle, fshader);
+	glGetProgramiv(handle, GL_LINK_STATUS, &linked);
+	if (!linked)
+		throw std::runtime_error("Failed to link shaders.");
+
+	// get camera uniforms
+	uModelViewMatrix = glGetUniformLocation(handle, "uModelViewMatrix");
+	uProjectionMatrix = glGetUniformLocation(handle, "uProjectionMatrix");
+
+	// post link
+	postLink();
+
+}
+
+void GLProgram::postLink() {}
+
+void GLProgram::setUniformMat4(GLuint unif, glm::mat4 m) {
+	glUseProgram(handle);
+	glUniformMatrix4fv(unif,
+		1, GL_FALSE, glm::value_ptr(m[0]));
+	glUseProgram(0);
+}
+
+
+void GLProgram::setModelView(glm::mat4 m) {
+	setUniformMat4(uModelViewMatrix, m);
+}
+
+void GLProgram::setProjection(glm::mat4 m) {
+	setUniformMat4(uProjectionMatrix, m);
+}
+
+GLProgram::operator GLuint() const { return handle; }
+
+GLProgram::~GLProgram() { glDeleteProgram(handle); }
+
+// PROGRAM INPUT //////////////////////////////////////////////////////////////////////////////
+
+ProgramInput::ProgramInput() {
+	// generate buffers
+	glGenBuffers(4, &vbo[0]);
+
+	// generate vertex array object
+	glGenVertexArrays(1, &handle);
+	glBindVertexArray(handle);
+
+	// positions
+	glBindBuffer(GL_ARRAY_BUFFER, vbo[0]);
+	glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, 0);
+
+	// normals
+	glBindBuffer(GL_ARRAY_BUFFER, vbo[1]);
+	glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, 0);
+
+	// texture coordinates
+	glBindBuffer(GL_ARRAY_BUFFER, vbo[2]);
+	glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, 0, 0);
+
+	// indices
+	glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vbo[3]);
+
+	glBindVertexArray(0);
+
+}
+
+void ProgramInput::bufferData(unsigned int index, void* buff, size_t size) {
+	glBindBuffer(GL_ARRAY_BUFFER, vbo[index]);
+	glBufferData(GL_ARRAY_BUFFER, size,
+		buff, GL_STATIC_DRAW);
+}
+void ProgramInput::setPositionData(float* buff, unsigned int len) {
+	bufferData(0, buff, sizeof(float) * len);
+}
+
+void ProgramInput::setNormalData(float* buff, unsigned int len) {
+	bufferData(1, buff, sizeof(float) * len);
+}
+
+void ProgramInput::setTextureData(float* buff, unsigned int len) {
+	bufferData(2, buff, sizeof(float) * len);
+}
+
+void ProgramInput::setIndexData(unsigned int* buff, unsigned int len) {
+	bufferData(3, buff, sizeof(unsigned int) * len);
+}
+
+ProgramInput::operator GLuint() const {
+	return handle;
+}
+
+
+ProgramInput::~ProgramInput() {
+	glDeleteBuffers(4, vbo);
+	glDeleteVertexArrays(1, &handle);
+}
+
+// SHADER PROGRAMS ////////////////////////////////////////////////////////////////////////////
+PhongShader::PhongShader() : GLProgram() {}
+void PhongShader::postLink() {
+	// get uniforms
+	uAlbedo = glGetUniformLocation(handle, "uAlbedo");
+	uAmbient = glGetUniformLocation(handle, "uAmbient");
+	uLight = glGetUniformLocation(handle, "uLight");
+}
+void PhongShader::setAlbedo(const glm::vec3& albedo) {
+	assert(uAlbedo >= 0);
+	glUseProgram(*this);
+	glUniform3f(uAlbedo, albedo[0], albedo[1], albedo[2]);
+	glUseProgram(0);
+}
+void PhongShader::setAmbient(const glm::vec3& ambient) {
+	assert(uAmbient >= 0);
+	glUseProgram(*this);
+	glUniform3f(uAmbient, ambient[0], ambient[1], ambient[2]);
+	glUseProgram(0);
+}
+void PhongShader::setLight(const glm::vec3& light) {
+	assert(uLight >= 0);
+	glUseProgram(*this);
+	glUniform3f(uLight, light[0], light[1], light[2]);
+	glUseProgram(0);
+}
+
+
+PickShader::PickShader() : GLProgram() {}
+
+void PickShader::postLink() {
+	// get uniforms
+	uTessFact = glGetUniformLocation(handle, "uTessFact");
+}
+
+void PickShader::setTessFact(unsigned int n) {
+	assert(uTessFact >= 0);
+	glUseProgram(*this);
+	glUniform1i(uTessFact, n);
+	glUseProgram(0);
+}

thirdparty/FastMassSpring/ClothApp/Shader.h

@@ -0,0 +1,87 @@
+#pragma once
+#include <fstream>
+#include <GL\glew.h>
+#include <glm/glm.hpp>
+
+class NonCopyable {
+private:
+	NonCopyable(const NonCopyable& other) = delete;
+	NonCopyable& operator=(const NonCopyable& other) = delete;
+
+public:
+	NonCopyable() {}
+};
+
+class GLShader : public NonCopyable {
+private:
+	GLuint handle; // Shader handle
+
+public:
+	GLShader(GLenum shaderType);
+	/*GLShader(GLenum shaderType, const char* source);
+	GLShader(GLenum shaderType, std::ifstream& source);*/
+	void compile(const char* source);
+	void compile(std::ifstream& source);
+	operator GLuint() const; // cast to GLuint
+	~GLShader();
+};
+
+class GLProgram : public NonCopyable {
+protected:
+	GLuint handle;
+	GLuint uModelViewMatrix, uProjectionMatrix;
+	void setUniformMat4(GLuint unif, glm::mat4 m);
+
+public:
+	GLProgram();
+	virtual void link(const GLShader& vshader, const GLShader& fshader);
+	virtual void postLink();
+	operator GLuint() const; // cast to GLuint
+	void setModelView(glm::mat4 m);
+	void setProjection(glm::mat4 m);
+
+	~GLProgram();
+};
+
+class ProgramInput : public NonCopyable {
+private:
+	GLuint handle; // vertex array object handle
+	GLuint vbo[4]; // vertex buffer object handles | position, normal, texture, index
+	void bufferData(unsigned int index, void* buff, size_t size);
+
+public:
+	ProgramInput();
+
+	void setPositionData(float* buff, unsigned int len);
+	void setNormalData(float* buff, unsigned int len);
+	void setTextureData(float* buff, unsigned int len);
+	void setIndexData(unsigned int* buff, unsigned int len);
+
+	operator GLuint() const; // cast to GLuint
+
+	~ProgramInput();
+};
+
+class PhongShader : public GLProgram {
+private:
+	// Albedo | Ambient Light | Light Direction
+	GLuint uAlbedo, uAmbient, uLight;
+
+public:
+	PhongShader();
+	virtual void postLink();
+	void setAlbedo(const glm::vec3& albedo);
+	void setAmbient(const glm::vec3& ambient);
+	void setLight(const glm::vec3& light);
+};
+
+
+class PickShader : public GLProgram {
+private:
+	GLuint uTessFact;
+
+public:
+	PickShader();
+	virtual void postLink();
+	void setTessFact(unsigned int n);
+};

thirdparty/FastMassSpring/ClothApp/UserInteraction.cpp

@@ -0,0 +1,48 @@
+#include "UserInteraction.h"
+#include <GL/glew.h>
+#include <iostream>
+#include <cmath>
+
+UserInteraction::UserInteraction(Renderer* renderer, CgPointFixNode* fixer, float* vbuff) 
+	: renderer(renderer), vbuff(vbuff), fixer(fixer), i(-1) {}
+
+void UserInteraction::setModelview(const glm::mat4& mv) { renderer->setModelview(mv); }
+void UserInteraction::setProjection(const glm::mat4& p) { renderer->setProjection(p); }
+
+void UserInteraction::grabPoint(int mouse_x, int mouse_y){
+	// render scene
+	glClearColor(0, 0, 0, 0);
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+	glDisable(GL_FRAMEBUFFER_SRGB);
+	renderer->draw();
+	glFlush();
+
+	// read color
+	color c(3);
+	glReadPixels(mouse_x, mouse_y, 1, 1, GL_RGB, GL_UNSIGNED_BYTE, &c[0]);
+	i = colorToIndex(c);
+	if (i != -1) fixer->fixPoint(i);
+
+	// return to normal state
+	glClearColor(0.25f, 0.25f, 0.25f, 0);
+	glEnable(GL_FRAMEBUFFER_SRGB);
+}
+
+void UserInteraction::releasePoint() { if (i == -1) return; fixer->releasePoint(i); i = -1; }
+void UserInteraction::movePoint(vec3 v) {
+	if (i == -1) return;
+	fixer->releasePoint(i);
+	for(int j = 0; j < 3; j++)
+		vbuff[3 * i + j] += v[j];
+	fixer->fixPoint(i);
+}
+
+GridMeshUI::GridMeshUI(Renderer* renderer, CgPointFixNode* fixer, float* vbuff, unsigned int n)
+	: UserInteraction(renderer, fixer, vbuff), n(n) {}
+
+int GridMeshUI::colorToIndex(color c) const {
+	if (c[2] != 51) return -1;
+	int vx = std::round((n - 1) * c[0] / 255.0);
+	int vy = std::round((n - 1) * c[1] / 255.0);
+	return n * vy + vx;
+}

thirdparty/FastMassSpring/ClothApp/UserInteraction.h

@@ -0,0 +1,35 @@
+#pragma once
+#include <glm/common.hpp>
+
+#include "MassSpringSolver.h"
+#include "Renderer.h"
+
+class UserInteraction {
+protected:
+	typedef glm::vec3 vec3;
+	typedef std::vector<unsigned char> color;
+
+	int i; // index of fixed point
+	float* vbuff; // vertex buffer
+	CgPointFixNode* fixer; // point fixer
+	Renderer* renderer; // pick shader renderer
+	virtual int colorToIndex(color c) const = 0;
+
+public:
+	UserInteraction(Renderer* renderer, CgPointFixNode* fixer, float* vbuff);
+
+	void setModelview(const glm::mat4& mv);
+	void setProjection(const glm::mat4& p);
+
+	void grabPoint(int mouse_x, int mouse_y); // grab point with color c
+	void movePoint(vec3 v); // move grabbed point along mouse
+	void releasePoint(); // release grabbed point;
+};
+
+class GridMeshUI : public UserInteraction {
+protected:
+	const unsigned int n; // grid width
+	virtual int colorToIndex(color c) const;
+public:
+	GridMeshUI(Renderer* renderer, CgPointFixNode* fixer, float* vbuff, unsigned int n);
+};

thirdparty/FastMassSpring/ClothApp/app.cpp

@@ -0,0 +1,471 @@
+#include <GL/glew.h>
+#include <GL/glut.h>
+#include <glm/gtc/matrix_transform.hpp>
+#include <glm/gtc/type_ptr.hpp>
+#include <stdexcept>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "Shader.h"
+#include "Mesh.h"
+#include "Renderer.h"
+#include "MassSpringSolver.h"
+#include "UserInteraction.h"
+
+// G L O B A L S ///////////////////////////////////////////////////////////////////
+
+// Window
+static int g_windowWidth = 640, g_windowHeight = 640;
+static bool g_mouseClickDown = false;
+static bool g_mouseLClickButton, g_mouseRClickButton, g_mouseMClickButton;
+static int g_mouseClickX;
+static int g_mouseClickY;
+
+// User Interaction
+static UserInteraction* UI;
+static Renderer* g_pickRenderer;
+
+// Constants
+static const float PI = glm::pi<float>();
+
+// Shader Handles
+static PhongShader* g_phongShader; // linked phong shader
+static PickShader* g_pickShader; // linked pick shader
+
+// Shader parameters
+static const glm::vec3 g_albedo(0.0f, 0.3f, 0.7f);
+static const glm::vec3 g_ambient(0.01f, 0.01f, 0.01f);
+static const glm::vec3 g_light(1.0f, 1.0f, -1.0f);
+
+// Mesh
+static Mesh* g_clothMesh; // halfedge data structure
+
+// Render Target
+static ProgramInput* g_render_target; // vertex, normal, texutre, index
+
+// Animation
+static const int g_fps = 60; // frames per second  | 60
+static const int g_iter = 5; // iterations per time step | 10
+static const int g_frame_time = 15; // approximate time for frame calculations | 15
+static const int g_animation_timer = (int) ((1.0f / g_fps) * 1000 - g_frame_time);
+
+// Mass Spring System
+static mass_spring_system* g_system;
+static MassSpringSolver* g_solver;
+
+// System parameters
+namespace SystemParam {
+	static const int n = 61; // must be odd, n * n = n_vertices | 61
+	static const float w = 2.0f; // width | 2.0f
+	static const float h = 0.008f; // time step, smaller for better results | 0.008f = 0.016f/2
+	static const float r = w / (n - 1) * 1.05f; // spring rest legnth
+	static const float k = 1.0f; // spring stiffness | 1.0f;
+	static const float m = 0.25f / (n * n); // point mass | 0.25f
+	static const float a = 0.993f; // damping, close to 1.0 | 0.993f
+	static const float g = 9.8f * m; // gravitational force | 9.8f
+}
+
+// Constraint Graph
+static CgRootNode* g_cgRootNode;
+
+// Scene parameters
+static const float g_camera_distance = 4.2f;
+
+// Scene matrices
+static glm::mat4 g_ModelViewMatrix;
+static glm::mat4 g_ProjectionMatrix;
+
+// F U N C T I O N S //////////////////////////////////////////////////////////////
+// state initialization
+static void initGlutState(int, char**);
+static void initGLState();
+
+static void initShaders(); // Read, compile and link shaders
+static void initCloth(); // Generate cloth mesh
+static void initScene(); // Generate scene matrices
+
+// demos
+static void demo_hang(); // curtain hanging from top corners
+static void demo_drop(); // curtain dropping on sphere
+static void(*g_demo)() = demo_drop;
+
+// glut callbacks
+static void display();
+static void reshape(int, int);
+static void mouse(int, int, int, int);
+static void motion(int, int);
+
+// draw cloth function
+static void drawCloth();
+static void animateCloth(int value);
+
+// scene update
+static void updateProjection();
+static void updateRenderTarget();
+
+// cleaning
+static void cleanUp();
+
+// error checks
+void checkGlErrors();
+
+
+
+// M A I N //////////////////////////////////////////////////////////////////////////
+int main(int argc, char** argv) {
+	try {
+		initGlutState(argc, argv);
+		glewInit();
+		initGLState();
+
+		initShaders();
+		initCloth();
+		initScene();
+
+		glutTimerFunc(g_animation_timer, animateCloth, 0);
+		glutMainLoop();
+
+		cleanUp();
+		return 0;
+	}
+	catch (const std::runtime_error& e) {
+		std::cout << "Exception caught: " << e.what() << std::endl;
+		return -1;
+	}
+}
+
+
+// S T A T E  I N I T I A L I Z A T O N /////////////////////////////////////////////
+static void initGlutState(int argc, char** argv) {
+	glutInit(&argc, argv);
+	glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
+	glutInitWindowSize(g_windowWidth, g_windowHeight);
+	glutCreateWindow("Cloth App");
+
+	glutDisplayFunc(display);
+	glutReshapeFunc(reshape);
+	glutMouseFunc(mouse);
+	glutMotionFunc(motion);
+}
+
+static void initGLState() {
+	glClearColor(0.25f, 0.25f, 0.25f, 0);
+	glClearDepth(1.);
+	glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
+	glPixelStorei(GL_PACK_ALIGNMENT, 1);
+	glEnable(GL_DEPTH_TEST);
+	glDepthFunc(GL_LESS);
+	glReadBuffer(GL_BACK);
+	glEnable(GL_FRAMEBUFFER_SRGB);
+
+	checkGlErrors();
+}
+
+static void initShaders() {
+	GLShader basic_vert(GL_VERTEX_SHADER);
+	GLShader phong_frag(GL_FRAGMENT_SHADER);
+	GLShader pick_frag(GL_FRAGMENT_SHADER);
+
+	basic_vert.compile(std::ifstream("./shaders/basic.vshader"));
+	phong_frag.compile(std::ifstream("./shaders/phong.fshader"));
+	pick_frag.compile(std::ifstream("./shaders/pick.fshader"));
+
+	g_phongShader = new PhongShader;
+	g_pickShader = new PickShader;
+	g_phongShader->link(basic_vert, phong_frag);
+	g_pickShader->link(basic_vert, pick_frag);
+
+	checkGlErrors();
+}
+
+static void initCloth() {
+	// short hand
+	const int n = SystemParam::n;
+	const float w = SystemParam::w;
+
+	// generate mesh
+	MeshBuilder meshBuilder;
+	meshBuilder.uniformGrid(w, n);
+	g_clothMesh = meshBuilder.getResult();
+
+	// fill program input
+	g_render_target = new ProgramInput;
+	g_render_target->setPositionData(g_clothMesh->vbuff(), g_clothMesh->vbuffLen());
+	g_render_target->setNormalData(g_clothMesh->nbuff(), g_clothMesh->nbuffLen());
+	g_render_target->setTextureData(g_clothMesh->tbuff(), g_clothMesh->tbuffLen());
+	g_render_target->setIndexData(g_clothMesh->ibuff(), g_clothMesh->ibuffLen());
+
+	// check errors
+	checkGlErrors();
+
+	// build demo system
+	g_demo();
+}
+
+static void initScene() {
+	g_ModelViewMatrix = glm::lookAt(
+		glm::vec3(0.618, -0.786, 0.3f) * g_camera_distance,
+		glm::vec3(0.0f, 0.0f, -1.0f),
+		glm::vec3(0.0f, 0.0f, 1.0f)
+	) * glm::translate(glm::mat4(1), glm::vec3(0.0f, 0.0f, SystemParam::w / 4));
+	updateProjection();
+}
+
+static void demo_hang() {
+	// short hand
+	const int n = SystemParam::n;
+
+	// initialize mass spring system
+	MassSpringBuilder massSpringBuilder;
+	massSpringBuilder.uniformGrid(
+		SystemParam::n,
+		SystemParam::h,
+		SystemParam::r,
+		SystemParam::k,
+		SystemParam::m,
+		SystemParam::a,
+		SystemParam::g
+	);
+	g_system = massSpringBuilder.getResult();
+
+	// initialize mass spring solver
+	g_solver = new MassSpringSolver(g_system, g_clothMesh->vbuff());
+
+	// deformation constraint parameters
+	const float tauc = 0.4f; // critical spring deformation | 0.4f
+	const unsigned int deformIter = 15; // number of iterations | 15
+
+	// initialize constraints
+	// spring deformation constraint
+	CgSpringDeformationNode* deformationNode =
+		new CgSpringDeformationNode(g_system, g_clothMesh->vbuff(), tauc, deformIter);
+	deformationNode->addSprings(massSpringBuilder.getShearIndex());
+	deformationNode->addSprings(massSpringBuilder.getStructIndex());
+
+	// fix top corners
+	CgPointFixNode* cornerFixer = new CgPointFixNode(g_system, g_clothMesh->vbuff());
+	cornerFixer->fixPoint(0);
+	cornerFixer->fixPoint(n - 1);
+
+	// initialize user interaction
+	g_pickRenderer = new Renderer();
+	g_pickRenderer->setProgram(g_pickShader);
+	g_pickRenderer->setProgramInput(g_render_target);
+	g_pickRenderer->setElementCount(g_clothMesh->ibuffLen());
+	g_pickShader->setTessFact(SystemParam::n);
+	CgPointFixNode* mouseFixer = new CgPointFixNode(g_system, g_clothMesh->vbuff());
+	UI = new GridMeshUI(g_pickRenderer, mouseFixer, g_clothMesh->vbuff(), n);
+
+	// build constraint graph
+	g_cgRootNode = new CgRootNode(g_system, g_clothMesh->vbuff());
+
+	// first layer
+	g_cgRootNode->addChild(deformationNode);
+
+	// second layer
+	deformationNode->addChild(cornerFixer);
+	deformationNode->addChild(mouseFixer);
+}
+
+static void demo_drop() {
+	// short hand
+	const int n = SystemParam::n;
+
+	// initialize mass spring system
+	MassSpringBuilder massSpringBuilder;
+	massSpringBuilder.uniformGrid(
+		SystemParam::n,
+		SystemParam::h,
+		SystemParam::r,
+		SystemParam::k,
+		SystemParam::m,
+		SystemParam::a,
+		SystemParam::g
+	);
+	g_system = massSpringBuilder.getResult();
+
+	// initialize mass spring solver
+	g_solver = new MassSpringSolver(g_system, g_clothMesh->vbuff());
+
+	// sphere collision constraint parameters
+	const float radius = 0.64f; // sphere radius | 0.64f
+	const Eigen::Vector3f center(0, 0, -1);// sphere center | (0, 0, -1)
+
+	// deformation constraint parameters
+	const float tauc = 0.12f; // critical spring deformation | 0.12f
+	const unsigned int deformIter = 15; // number of iterations | 15
+
+	// initialize constraints
+	// sphere collision constraint
+	CgSphereCollisionNode* sphereCollisionNode =
+		new CgSphereCollisionNode(g_system, g_clothMesh->vbuff(), radius, center);
+
+	// spring deformation constraint
+	CgSpringDeformationNode* deformationNode =
+		new CgSpringDeformationNode(g_system, g_clothMesh->vbuff(), tauc, deformIter);
+	deformationNode->addSprings(massSpringBuilder.getShearIndex());
+	deformationNode->addSprings(massSpringBuilder.getStructIndex());
+
+	// initialize user interaction
+	g_pickRenderer = new Renderer();
+	g_pickRenderer->setProgram(g_pickShader);
+	g_pickRenderer->setProgramInput(g_render_target);
+	g_pickRenderer->setElementCount(g_clothMesh->ibuffLen());
+	g_pickShader->setTessFact(SystemParam::n);
+	CgPointFixNode* mouseFixer = new CgPointFixNode(g_system, g_clothMesh->vbuff());
+	UI = new GridMeshUI(g_pickRenderer, mouseFixer, g_clothMesh->vbuff(), n);
+
+	// build constraint graph
+	g_cgRootNode = new CgRootNode(g_system, g_clothMesh->vbuff());
+
+	// first layer
+	g_cgRootNode->addChild(deformationNode);
+	g_cgRootNode->addChild(sphereCollisionNode);
+
+	// second layer
+	deformationNode->addChild(mouseFixer);
+}
+
+// G L U T  C A L L B A C K S //////////////////////////////////////////////////////
+static void display() {
+	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+	drawCloth();
+	glutSwapBuffers();
+
+	checkGlErrors();
+}
+
+static void reshape(int w, int h) {
+	g_windowWidth = w;
+	g_windowHeight = h;
+	glViewport(0, 0, w, h);
+	updateProjection();
+	glutPostRedisplay();
+}
+
+static void mouse(const int button, const int state, const int x, const int y) {
+	g_mouseClickX = x;
+	g_mouseClickY = g_windowHeight - y - 1;
+
+	g_mouseLClickButton |= (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN);
+	g_mouseRClickButton |= (button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN);
+	g_mouseMClickButton |= (button == GLUT_MIDDLE_BUTTON && state == GLUT_DOWN);
+
+	g_mouseLClickButton &= !(button == GLUT_LEFT_BUTTON && state == GLUT_UP);
+	g_mouseRClickButton &= !(button == GLUT_RIGHT_BUTTON && state == GLUT_UP);
+	g_mouseMClickButton &= !(button == GLUT_MIDDLE_BUTTON && state == GLUT_UP);
+
+	g_mouseClickDown = g_mouseLClickButton || g_mouseRClickButton || g_mouseMClickButton;
+
+	// TODO: move to UserInteraction class: add renderer member variable
+	// pick point
+	if (g_mouseLClickButton) {
+		UI->setModelview(g_ModelViewMatrix);
+		UI->setProjection(g_ProjectionMatrix);
+		UI->grabPoint(g_mouseClickX, g_mouseClickY);
+	}
+	else UI->releasePoint();
+}
+
+static void motion(const int x, const int y) {
+	const float dx = float(x - g_mouseClickX);
+	const float dy = float (-(g_windowHeight - y - 1 - g_mouseClickY));
+
+	if (g_mouseLClickButton) {
+		//glm::vec3 ux(g_ModelViewMatrix * glm::vec4(1, 0, 0, 0));
+		//glm::vec3 uy(g_ModelViewMatrix * glm::vec4(0, 1, 0, 0));
+		glm::vec3 ux(0, 1, 0);
+		glm::vec3 uy(0, 0, -1);
+		UI->movePoint(0.01f * (dx * ux + dy * uy));
+	}
+
+	g_mouseClickX = x;
+	g_mouseClickY = g_windowHeight - y - 1;
+}
+
+// C L O T H ///////////////////////////////////////////////////////////////////////
+static void drawCloth() {
+	Renderer renderer;
+	renderer.setProgram(g_phongShader);
+	renderer.setModelview(g_ModelViewMatrix);
+	renderer.setProjection(g_ProjectionMatrix);
+	g_phongShader->setAlbedo(g_albedo);
+	g_phongShader->setAmbient(g_ambient);
+	g_phongShader->setLight(g_light);
+	renderer.setProgramInput(g_render_target);
+	renderer.setElementCount(g_clothMesh->ibuffLen());
+	renderer.draw();
+}
+
+static void animateCloth(int value) {
+
+	// solve two time-steps
+	g_solver->solve(g_iter);
+	g_solver->solve(g_iter);
+
+	// fix points
+	CgSatisfyVisitor visitor;
+	visitor.satisfy(*g_cgRootNode);
+
+	// update normals
+	g_clothMesh->request_face_normals();
+	g_clothMesh->update_normals();
+	g_clothMesh->release_face_normals();
+
+	// update target
+	updateRenderTarget();
+
+	// redisplay
+	glutPostRedisplay();
+
+	// reset timer
+	glutTimerFunc(g_animation_timer, animateCloth, 0);
+}
+
+// S C E N E  U P D A T E ///////////////////////////////////////////////////////////
+static void updateProjection() {
+	g_ProjectionMatrix = glm::perspective(PI / 4.0f,
+		g_windowWidth * 1.0f / g_windowHeight, 0.01f, 1000.0f);
+}
+
+static void updateRenderTarget() {
+	// update vertex positions
+	g_render_target->setPositionData(g_clothMesh->vbuff(), g_clothMesh->vbuffLen());
+
+	// update vertex normals
+	g_render_target->setNormalData(g_clothMesh->nbuff(), g_clothMesh->vbuffLen());
+
+}
+
+// C L E A N  U P //////////////////////////////////////////////////////////////////
+static void cleanUp() {
+	// delete mesh
+	delete g_clothMesh;
+
+	// delete UI
+	delete g_pickRenderer;
+	delete UI;
+
+	// delete render target
+	delete g_render_target;
+
+	// delete mass-spring system
+	delete g_system;
+	delete g_solver;
+
+	// delete constraint graph
+	// TODO
+}
+
+// E R R O R S /////////////////////////////////////////////////////////////////////
+void checkGlErrors() {
+	const GLenum errCode = glGetError();
+
+	if (errCode != GL_NO_ERROR) {
+		std::string error("GL Error: ");
+		error += reinterpret_cast<const char*>(gluErrorString(errCode));
+		std::cerr << error << std::endl;
+		throw std::runtime_error(error);
+	}
+}

thirdparty/FastMassSpring/ClothApp/packages.config

@@ -0,0 +1,5 @@
+<?xml version="1.0" encoding="utf-8"?>
+<packages>
+  <package id="nupengl.core" version="0.1.0.1" targetFramework="native" />
+  <package id="nupengl.core.redist" version="0.1.0.1" targetFramework="native" />
+</packages>

thirdparty/FastMassSpring/ClothApp/shaders/basic.vshader

@@ -0,0 +1,19 @@
+#version 430
+
+uniform mat4 uModelViewMatrix;
+// uniform mat4 uNormalMatrix;
+uniform mat4 uProjectionMatrix;
+
+layout(location=0) in vec3 aPosition;
+layout(location=1) in vec3 aNormal;
+layout(location=2) in vec2 aTexCoord;
+
+out vec3 vNormal;
+out vec2 vTexCoord;
+
+void main(){
+    vNormal = aNormal;
+    vTexCoord = aTexCoord;
+    vec4 position = vec4(aPosition, 1.0);
+    gl_Position = uProjectionMatrix * uModelViewMatrix * position;
+}

thirdparty/FastMassSpring/ClothApp/shaders/phong.fshader

@@ -0,0 +1,22 @@
+#version 430
+
+in vec3 vNormal;
+out vec4 fragColor;
+
+uniform vec3 uAlbedo;
+uniform vec3 uAmbient;
+uniform vec3 uLight;
+void main(){
+    vec3 normal = normalize(vNormal);
+	vec3 albedo = uAlbedo;
+    if(!gl_FrontFacing) {
+		normal = -normal;
+		albedo = 1.0 - albedo;
+	}
+
+    vec3 toLight = normalize(-uLight);
+    float diffuse = max(0, dot(toLight, normal));
+
+    vec3 color = diffuse * albedo + uAmbient * albedo;
+    fragColor = vec4(color, 1.0);
+}

thirdparty/FastMassSpring/ClothApp/shaders/pick.fshader

@@ -0,0 +1,12 @@
+#version 430
+
+in vec2 vTexCoord;
+out vec4 fragColor;
+
+uniform int uTessFact;
+
+void main(){
+    float vx = round(vTexCoord.x * (uTessFact - 1));
+    float vy = round(vTexCoord.y * (uTessFact - 1));
+    fragColor = vec4(vx / (uTessFact - 1), vy / (uTessFact - 1), 0.2, 1.0);
+}

thirdparty/FastMassSpring/FastMassSpring.sln

@@ -0,0 +1,33 @@
+
+Microsoft Visual Studio Solution File, Format Version 12.00
+# Visual Studio 14
+VisualStudioVersion = 14.0.25420.1
+MinimumVisualStudioVersion = 10.0.40219.1
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "ClothApp", "ClothApp\ClothApp.vcxproj", "{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}"
+EndProject
+Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{D32F8712-F80F-46F9-85E1-FEB2BC718330}"
+	ProjectSection(SolutionItems) = preProject
+		readme.md = readme.md
+	EndProjectSection
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|x64 = Debug|x64
+		Debug|x86 = Debug|x86
+		Release|x64 = Release|x64
+		Release|x86 = Release|x86
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Debug|x64.ActiveCfg = Debug|x64
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Debug|x64.Build.0 = Debug|x64
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Debug|x86.ActiveCfg = Debug|Win32
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Debug|x86.Build.0 = Debug|Win32
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Release|x64.ActiveCfg = Release|x64
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Release|x64.Build.0 = Release|x64
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Release|x86.ActiveCfg = Release|Win32
+		{BBCA6691-6C35-4BFE-9CB9-13F16D29990F}.Release|x86.Build.0 = Release|Win32
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+EndGlobal

thirdparty/FastMassSpring/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2019 Samer Itani
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

thirdparty/FastMassSpring/readme.md

@@ -0,0 +1,16 @@
+A C++ implementation of *Fast Simulation of Mass-Spring Systems* [1], rendered with OpenGL.
+The dynamic inverse procedure described in [2] was implemented to constrain spring deformations and prevent the "super-elastic" effect when using large time-steps.
+
+### Dependencies
+* **OpenGL, freeGLUT, GLEW, GLM** for rendering.
+* **OpenMesh** for computing normals.
+* **Eigen** for sparse matrix algebra.
+
+### Demonstration
+![curtain hang](https://media.giphy.com/media/5EC1drLIyBuHxRC4I8/giphy.gif)
+![curtain drop](https://media.giphy.com/media/1zRbfGDmHTcn9RoUjy/giphy.gif)
+
+### References
+[1] Liu, T., Bargteil, A. W., Obrien, J. F., & Kavan, L. (2013). Fast simulation of mass-spring systems. *ACM Transactions on Graphics,32*(6), 1-7. doi:10.1145/2508363.2508406
+
+[2] Provot, X. (1995). Deformation constraints in a mass-spring modelto describe rigid cloth behavior. *InGraphics Interface* 1995,147–154.