introducing CG wrappers

Finite-Element-Tests/tests/CMakeLists.txt

@@ -3,3 +3,4 @@ add_subdirectory(MembranePinch2)
 add_subdirectory(MembranePinch3)
 add_subdirectory(MembranePinch4)
 add_subdirectory(MembranePinch5)
+add_subdirectory(MembranePinch6)

Finite-Element-Tests/tests/MembranePinch5/FiniteElementMesh.h

@@ -2,10 +2,6 @@
 
 #include "AnimatedMesh.h"
 
-#include "CGVectorWrapper.h"
-#include "CGSystemWrapper.h"
-#include "ConjugateGradient.h"
-
 #include <Eigen/Dense>
 
 template<class T>
@@ -102,22 +98,6 @@ struct FiniteElementMesh : public AnimatedMesh<T, 3>
         const int nParticles = m_particleX.size();
         std::vector<Vector2> force(nParticles, Vector2::Zero());
 
-        // std::vector<Vector2> x(nParticles, Vector2::Zero());
-        // std::vector<Vector2> b(nParticles, Vector2::Zero());
-        // std::vector<Vector2> q(nParticles, Vector2::Zero());
-        // std::vector<Vector2> s(nParticles, Vector2::Zero());
-        // std::vector<Vector2> r(nParticles, Vector2::Zero());
-        // CGVectorWrapper<Vector2> xWrapper(x);
-        // CGVectorWrapper<Vector2> bWrapper(b);
-        // CGVectorWrapper<Vector2> qWrapper(q);
-        // CGVectorWrapper<Vector2> sWrapper(s);
-        // CGVectorWrapper<Vector2> rWrapper(r);
-        // CGSystemWrapper<Vector2, FEMType> systemWrapper(*this);
-
-        // ConjugateGradient<T>::Solve(systemWrapper,
-        //     xWrapper, bWrapper, qWrapper, sWrapper, rWrapper,
-        //     1e-3, 50);
-
         addForce(force);
         for(int p = 0; p < nParticles; p++)
             m_particleX[p] += dt * m_particleV[p];

Finite-Element-Tests/tests/MembranePinch6/AnimatedMesh.h

@@ -0,0 +1,114 @@
+#pragma once
+
+#include "pxr/pxr.h"
+
+#include "pxr/usd/sdf/layer.h"
+#include "pxr/usd/sdf/path.h"
+#include "pxr/usd/usd/stage.h"
+#include "pxr/usd/usdGeom/mesh.h"
+#include "pxr/base/vt/array.h"
+
+#include "pxr/base/gf/range3f.h"
+
+#include <Eigen/Dense>
+
+PXR_NAMESPACE_USING_DIRECTIVE
+
+template<class T, int d> // d is the dimension of the mesh elements, e.g. 3 for triangles, 4 for quads
+struct AnimatedMesh
+{
+    using Vector2 = Eigen::Matrix< T , 2 , 1>;
+
+    SdfLayerRefPtr m_layer;
+    UsdStageRefPtr m_stage;
+    UsdGeomMesh m_mesh;
+    UsdAttribute m_pointsAttribute;
+
+    GfRange3f m_extent;
+    int m_lastFrame;
+
+    std::vector<std::array<int, d>> m_meshElements;
+    std::vector<Vector2> m_particleX;
+
+    AnimatedMesh()
+        :m_lastFrame(-1)
+    {}
+
+    void initializeUSD(const std::string filename)
+    {
+        // Create the layer to populate.
+        m_layer = SdfLayer::CreateNew(filename);
+
+        // Create a UsdStage with that root layer.
+        m_stage = UsdStage::Open(m_layer);
+    }
+
+    void initializeTopology()
+    {
+        // Create a mesh for this surface
+        m_mesh = UsdGeomMesh::Define(m_stage, SdfPath("/MeshSurface"));        
+
+        // Create appropriate buffers for vertex counts and indices, and populate them
+        VtIntArray faceVertexCounts, faceVertexIndices;
+        for (const auto& element : m_meshElements) {
+            faceVertexCounts.push_back(element.size());
+            for (const auto& vertex : element)
+                faceVertexIndices.push_back(vertex);
+        }
+
+        // Now set the attributes
+        m_mesh.GetFaceVertexCountsAttr().Set(faceVertexCounts);
+        m_mesh.GetFaceVertexIndicesAttr().Set(faceVertexIndices);
+    }
+
+    void initializeParticles()
+    {
+        // Grab the points (Positions) attribute, and indicate it is time-varying
+        m_pointsAttribute = m_mesh.GetPointsAttr();
+        m_pointsAttribute.SetVariability(SdfVariabilityVarying);
+    }
+
+    void writeFrame(const int frame)
+    {
+        std::cout << "Writing frame " << frame << " ..." << std::endl;
+
+        // Check that there are any particles to write at all
+        if (m_particleX.empty())
+            throw std::logic_error("Empty array of input vertices");
+
+        // Check that frames have been written in sequence
+        if(frame != m_lastFrame+1)
+            throw std::logic_error("Non-consequtive frame sequence requested in writeFrame()");
+        m_lastFrame = frame;
+
+        // Update extent
+        for (const auto& pt : m_particleX)
+            m_extent.UnionWith(GfVec3f(pt[0],pt[1],T()));
+
+        // Copy particleX into VtVec3fArray for Usd
+        VtVec3fArray usdPoints(m_particleX.size());
+        for(int p = 0; p < m_particleX.size(); p++)
+            usdPoints[p] = GfVec3f(m_particleX[p][0],m_particleX[p][1],T());
+
+        // Write the points attribute for the given frame
+        m_pointsAttribute.Set(usdPoints, (double) frame);
+    }
+
+    void writeUSD()
+    {
+        // Set up the timecode
+        m_stage->SetStartTimeCode(0.);
+        m_stage->SetEndTimeCode((double) m_lastFrame);
+
+        // Set the effective extent
+        VtVec3fArray extentArray(2);
+        extentArray[0] = m_extent.GetMin();
+        extentArray[1] = m_extent.GetMax();
+        m_mesh.GetExtentAttr().Set(extentArray);
+
+        // Save USD file
+        m_stage->GetRootLayer()->Save();
+        std::cout << "USD file saved!" << std::endl;
+    }
+};
+

Finite-Element-Tests/tests/MembranePinch6/CMakeLists.txt

@@ -0,0 +1,26 @@
+add_executable (MembranePinch6
+        main.cpp
+        )
+
+find_package(PythonLibs)
+
+set(EIGEN3_INC_DIR /usr/include/eigen3)
+
+target_include_directories(MembranePinch6
+        PUBLIC
+            ${USD_INC_DIR}
+            ${EIGEN3_INC_DIR}
+            ${PYTHON_INCLUDE_PATH}
+        )
+
+target_link_libraries(
+        MembranePinch6
+                ${USD_LIB_DIR}/libgf.so
+                ${USD_LIB_DIR}/libsdf.so
+                ${USD_LIB_DIR}/libtf.so
+                ${USD_LIB_DIR}/libusd.so
+                ${USD_LIB_DIR}/libusdGeom.so
+                ${USD_LIB_DIR}/libvt.so
+                ${USD_LIB_DIR}/libboost_python.so # todo find library
+                ${PYTHON_LIBRARY}
+        )

Finite-Element-Tests/tests/MembranePinch6/FiniteElementMesh.h

@@ -0,0 +1,145 @@
+#pragma once
+
+#include "AnimatedMesh.h"
+
+#include "CGVectorWrapper.h"
+#include "CGSystemWrapper.h"
+#include "ConjugateGradient.h"
+
+#include <Eigen/Dense>
+
+template<class T>
+struct FiniteElementMesh : public AnimatedMesh<T, 3>
+{
+    using Base = AnimatedMesh<T, 3>;
+    using Base::m_meshElements;
+    using Base::m_particleX;
+    using Vector2 = typename Base::Vector2;
+    using Matrix22 = Eigen::Matrix< T , 2 , 2>;
+
+    int m_nFrames;
+    int m_subSteps;
+    T m_frameDt;
+
+    const T m_density;
+    const T m_mu;
+    const T m_lambda;
+    const T m_rayleighCoefficient;
+
+    std::vector<T> m_particleMass;
+    std::vector<Vector2> m_particleV;
+    std::vector<Matrix22> m_DmInverse;
+    std::vector<T> m_restVolume;
+
+    FiniteElementMesh(const T density, const T mu, const T lambda, const T rayleighCoefficient)
+        :m_density(density), m_mu(mu), m_lambda(lambda), m_rayleighCoefficient(rayleighCoefficient)
+    {}
+
+    void initializeUndeformedConfiguration()
+    {
+        // Initialize rest shape and particle mass (based on constant density)
+        m_particleMass.resize(m_particleX.size(), T()); // Initialize all particle masses to zero
+        for(const auto& element: m_meshElements)
+        {
+            Matrix22 Dm;
+            for(int j = 0; j < 2; j++)
+                Dm.col(j) = m_particleX[element[j+1]]-m_particleX[element[0]];
+            T restVolume = .5 * Dm.determinant();
+            if(restVolume < 0)
+                throw std::logic_error("Inverted element");
+            m_DmInverse.emplace_back(Dm.inverse());
+            m_restVolume.push_back(restVolume);
+            T elementMass = m_density * restVolume;
+            for(const int v: element)
+                m_particleMass[v] += (1./3.) * elementMass;
+        }
+    }
+
+    void addElasticForce(std::vector<Vector2>& f) const
+    {
+        for(int e = 0; e < m_meshElements.size(); e++)
+        {
+            const auto& element = m_meshElements[e];
+
+            // Linear Elasticity
+            Matrix22 Ds;
+            for(int j = 0; j < 2; j++)
+                Ds.col(j) = m_particleX[element[j+1]]-m_particleX[element[0]];
+            Matrix22 F = Ds * m_DmInverse[e];
+
+            Matrix22 strain = .5 * (F + F.transpose()) - Matrix22::Identity();
+            Matrix22 P = 2. * m_mu * strain + m_lambda * strain.trace() * Matrix22::Identity();
+
+            Matrix22 H = -m_restVolume[e] * P * m_DmInverse[e].transpose();
+
+            for(int j = 0; j < 2; j++){
+                f[element[j+1]] += H.col(j);
+                f[element[0]] -= H.col(j);
+            }
+        }
+    }
+
+    void addProductWithStiffnessMatrix(std::vector<Vector2>& w, std::vector<Vector2>& f, const T scale) const
+    {
+        for(int e = 0; e < m_meshElements.size(); e++)
+        {
+            const auto& element = m_meshElements[e];
+
+            // Linear Damping
+            Matrix22 Ds_dot;
+            for(int j = 0; j < 2; j++)
+                Ds_dot.col(j) = w[element[j+1]]-w[element[0]];
+            Matrix22 F_dot = Ds_dot * m_DmInverse[e];
+
+            Matrix22 strain_rate = .5 * (F_dot + F_dot.transpose());
+            Matrix22 P_damping = scale * (2. * m_mu * strain_rate + m_lambda * strain_rate.trace() * Matrix22::Identity());
+
+            Matrix22 H_damping = -m_restVolume[e] * P_damping * m_DmInverse[e].transpose();
+
+            for(int j = 0; j < 2; j++){
+                f[element[j+1]] += H_damping.col(j);
+                f[element[0]] -= H_damping.col(j);
+            }
+        }
+    }
+
+    void simulateSubstep(const T dt)
+    {
+        using FEMType = FiniteElementMesh<T>;
+
+        const int nParticles = m_particleX.size();
+        std::vector<Vector2> force(nParticles, Vector2::Zero());
+
+        std::vector<Vector2> x(nParticles, Vector2::Zero());
+        std::vector<Vector2> b(nParticles, Vector2::Zero());
+        std::vector<Vector2> q(nParticles, Vector2::Zero());
+        std::vector<Vector2> s(nParticles, Vector2::Zero());
+        std::vector<Vector2> r(nParticles, Vector2::Zero());
+        CGVectorWrapper<Vector2> xWrapper(x);
+        CGVectorWrapper<Vector2> bWrapper(b);
+        CGVectorWrapper<Vector2> qWrapper(q);
+        CGVectorWrapper<Vector2> sWrapper(s);
+        CGVectorWrapper<Vector2> rWrapper(r);
+        CGSystemWrapper<Vector2, FEMType> systemWrapper(*this);
+
+        // ConjugateGradient<T>::Solve(systemWrapper,
+        //     xWrapper, bWrapper, qWrapper, sWrapper, rWrapper,
+        //     1e-3, 50);
+
+        addElasticForce(force);
+        addProductWithStiffnessMatrix(m_particleV, force, m_rayleighCoefficient);
+        for(int p = 0; p < nParticles; p++)
+            m_particleX[p] += dt * m_particleV[p];
+        for(int p = 0; p < nParticles; p++)
+            m_particleV[p] += (dt / m_particleMass[p]) * force[p];
+    }
+
+    void simulateFrame(const int frame)
+    {
+        T stepDt = m_frameDt / (T) m_subSteps;
+
+        for(int step = 1; step <= m_subSteps; step++)
+            simulateSubstep(stepDt);
+    }
+};
+