diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..8a5f669 --- /dev/null +++ b/.clang-format @@ -0,0 +1,246 @@ +--- +Language: Cpp +# BasedOnStyle: Mozilla +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignArrayOfStructures: None +AlignConsecutiveAssignments: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: true +AlignConsecutiveBitFields: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveDeclarations: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveMacros: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCompound: false + AlignFunctionPointers: false + PadOperators: false +AlignConsecutiveShortCaseStatements: + Enabled: false + AcrossEmptyLines: false + AcrossComments: false + AlignCaseColons: false +AlignEscapedNewlines: Right +AlignOperands: Align +AlignTrailingComments: + Kind: Always + OverEmptyLines: 0 +AllowAllArgumentsOnNextLine: true +AllowAllParametersOfDeclarationOnNextLine: false +AllowBreakBeforeNoexceptSpecifier: Never +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortCompoundRequirementOnASingleLine: true +AllowShortEnumsOnASingleLine: true +AllowShortFunctionsOnASingleLine: Inline +AllowShortIfStatementsOnASingleLine: Never +AllowShortLambdasOnASingleLine: All +AllowShortLoopsOnASingleLine: false +AlwaysBreakAfterDefinitionReturnType: TopLevel +AlwaysBreakAfterReturnType: TopLevel +AlwaysBreakBeforeMultilineStrings: false +AlwaysBreakTemplateDeclarations: Yes +AttributeMacros: + - __capability +BinPackArguments: false +BinPackParameters: false +BitFieldColonSpacing: Both +BraceWrapping: + AfterCaseLabel: false + AfterClass: true + AfterControlStatement: Never + AfterEnum: true + AfterExternBlock: true + AfterFunction: true + AfterNamespace: false + AfterObjCDeclaration: false + AfterStruct: true + AfterUnion: true + BeforeCatch: false + BeforeElse: false + BeforeLambdaBody: false + BeforeWhile: false + IndentBraces: false + SplitEmptyFunction: true + SplitEmptyRecord: false + SplitEmptyNamespace: true +BreakAdjacentStringLiterals: true +BreakAfterAttributes: Leave +BreakAfterJavaFieldAnnotations: false +BreakArrays: true +BreakBeforeBinaryOperators: None +BreakBeforeConceptDeclarations: Always +BreakBeforeBraces: Mozilla +BreakBeforeInlineASMColon: OnlyMultiline +BreakBeforeTernaryOperators: true +BreakConstructorInitializers: BeforeComma +BreakInheritanceList: BeforeComma +BreakStringLiterals: true +ColumnLimit: 120 +CommentPragmas: '^ IWYU pragma:' +CompactNamespaces: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: false +DerivePointerAlignment: false +DisableFormat: false +EmptyLineAfterAccessModifier: Never +EmptyLineBeforeAccessModifier: LogicalBlock +ExperimentalAutoDetectBinPacking: false +FixNamespaceComments: false +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IfMacros: + - KJ_IF_MAYBE +IncludeBlocks: Preserve +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 2 + SortPriority: 0 + CaseSensitive: false + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 3 + SortPriority: 0 + CaseSensitive: false + - Regex: '.*' + Priority: 1 + SortPriority: 0 + CaseSensitive: false +IncludeIsMainRegex: '(Test)?$' +IncludeIsMainSourceRegex: '' +IndentAccessModifiers: false +IndentCaseBlocks: false +IndentCaseLabels: true +IndentExternBlock: AfterExternBlock +IndentGotoLabels: true +IndentPPDirectives: None +IndentRequiresClause: true +IndentWidth: 4 +IndentWrappedFunctionNames: false +InsertBraces: false +InsertNewlineAtEOF: false +InsertTrailingCommas: None +IntegerLiteralSeparator: + Binary: 0 + BinaryMinDigits: 0 + Decimal: 0 + DecimalMinDigits: 0 + Hex: 0 + HexMinDigits: 0 +JavaScriptQuotes: Leave +JavaScriptWrapImports: true +KeepEmptyLinesAtTheStartOfBlocks: true +KeepEmptyLinesAtEOF: false +LambdaBodyIndentation: Signature +LineEnding: DeriveLF +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBinPackProtocolList: Auto +ObjCBlockIndentWidth: 4 +ObjCBreakBeforeNestedBlockParam: true +ObjCSpaceAfterProperty: true +ObjCSpaceBeforeProtocolList: false +PackConstructorInitializers: BinPack +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 19 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakOpenParenthesis: 0 +PenaltyBreakScopeResolution: 500 +PenaltyBreakString: 1000 +PenaltyBreakTemplateDeclaration: 10 +PenaltyExcessCharacter: 1000000 +PenaltyIndentedWhitespace: 0 +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Left +PPIndentWidth: -1 +QualifierAlignment: Leave +ReferenceAlignment: Pointer +ReflowComments: true +RemoveBracesLLVM: false +RemoveParentheses: Leave +RemoveSemicolon: false +RequiresClausePosition: OwnLine +RequiresExpressionIndentation: OuterScope +SeparateDefinitionBlocks: Leave +ShortNamespaceLines: 1 +SkipMacroDefinitionBody: false +SortIncludes: CaseSensitive +SortJavaStaticImport: Before +SortUsingDeclarations: LexicographicNumeric +SpaceAfterCStyleCast: false +SpaceAfterLogicalNot: false +SpaceAfterTemplateKeyword: false +SpaceAroundPointerQualifiers: Default +SpaceBeforeAssignmentOperators: true +SpaceBeforeCaseColon: false +SpaceBeforeCpp11BracedList: false +SpaceBeforeCtorInitializerColon: true +SpaceBeforeInheritanceColon: true +SpaceBeforeJsonColon: false +SpaceBeforeParens: ControlStatements +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterForeachMacros: true + AfterFunctionDefinitionName: false + AfterFunctionDeclarationName: false + AfterIfMacros: true + AfterOverloadedOperator: false + AfterPlacementOperator: true + AfterRequiresInClause: false + AfterRequiresInExpression: false + BeforeNonEmptyParentheses: false +SpaceBeforeRangeBasedForLoopColon: true +SpaceBeforeSquareBrackets: false +SpaceInEmptyBlock: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: Never +SpacesInContainerLiterals: true +SpacesInLineCommentPrefix: + Minimum: 1 + Maximum: -1 +SpacesInParens: Never +SpacesInParensOptions: + InCStyleCasts: false + InConditionalStatements: false + InEmptyParentheses: false + Other: false +SpacesInSquareBrackets: false +Standard: Latest +StatementAttributeLikeMacros: + - Q_EMIT +StatementMacros: + - Q_UNUSED + - QT_REQUIRE_VERSION +TabWidth: 4 +UseTab: Never +VerilogBreakBetweenInstancePorts: true +WhitespaceSensitiveMacros: + - BOOST_PP_STRINGIZE + - CF_SWIFT_NAME + - NS_SWIFT_NAME + - PP_STRINGIZE + - STRINGIZE +... + diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 9514767..71ab9b8 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -12,7 +12,7 @@ jobs: steps: - name: clone - uses: actions/checkout@v4 # Action to check out your repository + uses: actions/checkout@v4 - name: dependencies run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index d26e645..34fecf3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ project(cloth_sim_tutorial LANGUAGES CXX) option(CLOTH_DEBUG "run with debug options" OFF) option(CLOTH_TESTING_ONLY "only compile unit tests" OFF) -set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake) add_definitions(-DCLOTH_ROOT_DIR="${CMAKE_CURRENT_SOURCE_DIR}") diff --git a/README.md b/README.md index 20c8b21..0cefcb6 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ This results in a block-sparse linear system for which the velocity deltas are s The B&W method is conceptually straight forward. However, there are numerous gotchas not apparent in the original references. In particular, working out the Jacobians of some models can be nontrivial. Even when correct, the derivatives can still produce hard-to-debug instabilities. -The [dynamic deformables course notes](https://doi.org/10.1145/3388769.3407490) describes this well, and I highly recommend reviewing the later chapters on cloth simulation for a thorough review of the topic. +The [dynamic deformables course notes](https://doi.org/10.1145/3388769.3407490) describes this well, and I highly recommend reviewing the later chapters on cloth simulation. Nonetheless, implicit time integration remains a valuable and sometimes necessary tool for modern graphics applications. Recently, optimization-based methods have grown in popularity. This approach formulates the Eqs. of motion as an [iteratively minimized objective function](https://en.wikipedia.org/wiki/Optimization_problem). @@ -26,6 +26,8 @@ Below we describe a simple implicit solver for mass-spring systems using gradien This is meant solely as a tutorial to optimization-based time stepping, and not for practical applications. For that, you'll want to consider higher-order optimization algorithms or [local-global type solvers](https://github.com/alecjacobson/computer-graphics-mass-spring-systems). +

+ ## Implementation The objective we aim to minimize is the sum of internal elastic energy potentials plus a quadratic penalty of linear momentum: @@ -34,7 +36,7 @@ $`\bar{x} = x + hv + h^2 M^{-1} f_{ext}`$, $`g(x) = \frac{1}{2h^2}\|M^{1/2}(x-\bar{x})\|^2 + \sum E(x)`$, -for vertex locations $x$, time step (sec) $h$, diagonal mass matrix $M$, velocities $v$, and external forces (gravity, wind, etc...) +for vertex locations $x$, time step (seconds) $h$, diagonal mass matrix $M$, velocities $v$, and external forces (gravity, wind, etc...) $`f_{ext}`$. The term $`\bar{x}`$ is the *explicit predictor* that is computed at the beginning of the time step. Computing a frame of animation amounts to iteratively minimizing the above objective until some convergence criteria is met (e.g., the norm of the gradient is below some threshold). @@ -91,4 +93,5 @@ There are many ways to improve the solver illustrated above: - Collision against triangle meshes with a [signed distance field](https://github.com/InteractiveComputerGraphics/TriangleMeshDistance) - Better energy models with triangular elements ([Dynamic Deformables appendix D](https://doi.org/10.1145/3388769.3407490)) - Preconditioning or accelerated versions of gradient descent +- Non-conservative forces like friction - Optimizers with better convergence like [L-BFGS](https://en.wikipedia.org/wiki/Limited-memory_BFGS) or [projected Newton](https://doi.org/10.1145/1073368.1073394) diff --git a/data/sphere.png b/data/sphere.png new file mode 100644 index 0000000..08b8b8e Binary files /dev/null and b/data/sphere.png differ diff --git a/src/ClothAssert.hpp b/src/ClothAssert.hpp index 7975731..11b7441 100644 --- a/src/ClothAssert.hpp +++ b/src/ClothAssert.hpp @@ -4,29 +4,28 @@ #ifndef CLOTHASSERT_HPP #define CLOTHASSERT_HPP 1 -#include -#include +#include #include +#include -namespace cloth -{ +namespace cloth { -static inline void AssertHandler(bool cond, const std::string& file, const int& line) +static inline void +AssertHandler(bool cond, const std::string& file, const int& line) { - if (!cond) - { - std::string err_msg = "Assertion failed in "+file+" line "+std::to_string(line); - throw std::runtime_error(err_msg.c_str()); - } + if (!cond) { + std::string err_msg = "Assertion failed in " + file + " line " + std::to_string(line); + throw std::runtime_error(err_msg.c_str()); + } } -static inline void AssertHandlerMsg(bool cond, const std::string& file, const int& line, const std::string &msg) +static inline void +AssertHandlerMsg(bool cond, const std::string& file, const int& line, const std::string& msg) { - if (!cond) - { - std::string err_msg = "Assertion failed in "+file+" line "+std::to_string(line)+": "+msg; - throw std::runtime_error(err_msg.c_str()); - } + if (!cond) { + std::string err_msg = "Assertion failed in " + file + " line " + std::to_string(line) + ": " + msg; + throw std::runtime_error(err_msg.c_str()); + } } } @@ -37,12 +36,9 @@ static inline void AssertHandlerMsg(bool cond, const std::string& file, const in #define ClothAssert_withmsg(cond, msg) AssertHandlerMsg(cond, std::string(__FILE__), __LINE__, msg) #define ClothAssert_nomsg(cond) AssertHandler(cond, std::string(__FILE__), __LINE__) -#define ClothAssert_stripargs(xx,cond,msg,FUNC, ...) FUNC +#define ClothAssert_stripargs(xx, cond, msg, FUNC, ...) FUNC -#define ClothAssert(...) \ - ClothAssert_stripargs(,##__VA_ARGS__,\ - ClothAssert_withmsg(__VA_ARGS__),\ - ClothAssert_nomsg(__VA_ARGS__),\ - ) +#define ClothAssert(...) \ + ClothAssert_stripargs(, ##__VA_ARGS__, ClothAssert_withmsg(__VA_ARGS__), ClothAssert_nomsg(__VA_ARGS__), ) #endif diff --git a/src/ClothMesh.hpp b/src/ClothMesh.hpp index 214ab9d..1a60b88 100644 --- a/src/ClothMesh.hpp +++ b/src/ClothMesh.hpp @@ -8,131 +8,116 @@ #include #include -namespace cloth -{ +namespace cloth { class ClothMesh { -public: - - Eigen::MatrixXd V; // rest vertices [nx3] - Eigen::MatrixXi F; // faces [fx3] - Eigen::MatrixXi E; // edges [ex3] - Eigen::VectorXd E0; // edge rest length - Eigen::MatrixXi B; // bend edges [bx3] - Eigen::VectorXd B0; // bend rest length - Eigen::VectorXd masses; // per-vertex masses - Eigen::VectorXi V_offset; // per-mesh index into V - Eigen::VectorXi F_offset; // per-mesh index into F - - // Adds a mesh to the collection of deformables and returns mesh number - inline int add_mesh(const Eigen::MatrixXd &inV, const Eigen::MatrixXi &inF) - { - int mesh_idx = -1; - ClothAssert(inV.rows() > 0 && inV.cols()==3); - ClothAssert(inF.rows() > 0 && inF.cols()==3); - - // First mesh? - if (V_offset.size()==0) - { - using namespace Eigen; - V_offset = VectorXi::Zero(2); - F_offset = VectorXi::Zero(2); - V_offset[1] = inV.rows(); - F_offset[1] = inF.rows(); - V = inV; - F = inF; - mesh_idx = 0; - } - else - { - mesh_idx = V_offset.size()-1; - V_offset.conservativeResize(mesh_idx+2); - F_offset.conservativeResize(mesh_idx+2); - V_offset[mesh_idx+1] = inV.rows(); - F_offset[mesh_idx+1] = inF.rows(); - V.conservativeResize(V.rows() + inV.rows(), 3); - V.bottomRows(inV.rows()) = inV; - F.conservativeResize(F.rows() + inF.rows(), 3); - F.bottomRows(inF.rows()) = inF; - } - - masses = Eigen::VectorXd::Constant(V.rows(), 0.1/V.rows()); - - compute_edges(); - return mesh_idx; - } - - // Computes edges and cross-edges from faces - inline void compute_edges() - { - using namespace Eigen; - std::map, std::pair> edges; // [edge, bend] - typedef std::map, std::pair>::iterator EdgeIter; - int n_bend = 0; // keep track of valid bend edges (non-border) - - int nf = F.rows(); - for (int i=0; i edge = std::make_pair(F(i,j), F(i,(j+1)%3)); - int adj = F(i,(j+2)%3); - if (edge.second < edge.first) - { - std::swap(edge.first, edge.second); - } - - EdgeIter it = edges.find(edge); - if (it == edges.end()) - { - edges[edge] = std::make_pair(adj, -1); - } - else - { - n_bend++; - std::pair &adj_pair = it->second; - adj_pair.second = adj; - if (adj_pair.second < adj_pair.first) - { - std::swap(adj_pair.first, adj_pair.second); - } - } - } - } - - // Loop over map and create eigen buffers - E = MatrixXi::Zero(edges.size(), 2); - E0 = VectorXd::Constant(edges.size(), -1); - B = MatrixXi::Zero(n_bend, 2); - B0 = VectorXd::Constant(n_bend, -1); - - int edge_idx = 0; - int bend_idx = 0; - for(EdgeIter it = edges.begin(); it != edges.end(); ++it, ++edge_idx) - { - const std::pair &e = it->first; - E(edge_idx,0) = e.first; - E(edge_idx,1) = e.second; - E0[edge_idx] = (V.row(e.first) - V.row(e.second)).norm(); - - const std::pair &b = it->second; - if (b.first >= 0 && b.second >= 0) - { - ClothAssert(bend_idx < B.rows()); - B(bend_idx,0) = b.first; - B(bend_idx,1) = b.second; - B0[bend_idx] = (V.row(b.first) - V.row(b.second)).norm(); - bend_idx++; - } - } - - // Check to make sure we set every edge - ClothAssert(E0.minCoeff() >= 0); - ClothAssert(B0.minCoeff() >= 0); - - } // end compute edges - + public: + Eigen::MatrixXd V; // rest vertices [nx3] + Eigen::MatrixXi F; // faces [fx3] + Eigen::MatrixXi E; // edges [ex3] + Eigen::VectorXd E0; // edge rest length + Eigen::MatrixXi B; // bend edges [bx3] + Eigen::VectorXd B0; // bend rest length + Eigen::VectorXd masses; // per-vertex masses + Eigen::VectorXi V_offset; // per-mesh index into V + Eigen::VectorXi F_offset; // per-mesh index into F + + // Adds a mesh to the collection of deformables and returns mesh number + inline int add_mesh(const Eigen::MatrixXd& inV, const Eigen::MatrixXi& inF) + { + int mesh_idx = -1; + ClothAssert(inV.rows() > 0 && inV.cols() == 3); + ClothAssert(inF.rows() > 0 && inF.cols() == 3); + + // First mesh? + if (V_offset.size() == 0) { + using namespace Eigen; + V_offset = VectorXi::Zero(2); + F_offset = VectorXi::Zero(2); + V_offset[1] = inV.rows(); + F_offset[1] = inF.rows(); + V = inV; + F = inF; + mesh_idx = 0; + } else { + mesh_idx = V_offset.size() - 1; + V_offset.conservativeResize(mesh_idx + 2); + F_offset.conservativeResize(mesh_idx + 2); + V_offset[mesh_idx + 1] = inV.rows(); + F_offset[mesh_idx + 1] = inF.rows(); + V.conservativeResize(V.rows() + inV.rows(), 3); + V.bottomRows(inV.rows()) = inV; + F.conservativeResize(F.rows() + inF.rows(), 3); + F.bottomRows(inF.rows()) = inF; + } + + masses = Eigen::VectorXd::Constant(V.rows(), 0.1 / V.rows()); + + compute_edges(); + return mesh_idx; + } + + // Computes edges and cross-edges from faces + inline void compute_edges() + { + using namespace Eigen; + std::map, std::pair> edges; // [edge, bend] + typedef std::map, std::pair>::iterator EdgeIter; + int n_bend = 0; // keep track of valid bend edges (non-border) + + int nf = F.rows(); + for (int i = 0; i < nf; ++i) { + for (int j = 0; j < 3; ++j) { + std::pair edge = std::make_pair(F(i, j), F(i, (j + 1) % 3)); + int adj = F(i, (j + 2) % 3); + if (edge.second < edge.first) { + std::swap(edge.first, edge.second); + } + + EdgeIter it = edges.find(edge); + if (it == edges.end()) { + edges[edge] = std::make_pair(adj, -1); + } else { + n_bend++; + std::pair& adj_pair = it->second; + adj_pair.second = adj; + if (adj_pair.second < adj_pair.first) { + std::swap(adj_pair.first, adj_pair.second); + } + } + } + } + + // Loop over map and create eigen buffers + E = MatrixXi::Zero(edges.size(), 2); + E0 = VectorXd::Constant(edges.size(), -1); + B = MatrixXi::Zero(n_bend, 2); + B0 = VectorXd::Constant(n_bend, -1); + + int edge_idx = 0; + int bend_idx = 0; + for (EdgeIter it = edges.begin(); it != edges.end(); ++it, ++edge_idx) { + const std::pair& e = it->first; + E(edge_idx, 0) = e.first; + E(edge_idx, 1) = e.second; + E0[edge_idx] = (V.row(e.first) - V.row(e.second)).norm(); + + const std::pair& b = it->second; + if (b.first >= 0 && b.second >= 0) { + ClothAssert(bend_idx < B.rows()); + B(bend_idx, 0) = b.first; + B(bend_idx, 1) = b.second; + B0[bend_idx] = (V.row(b.first) - V.row(b.second)).norm(); + bend_idx++; + } + } + + // Check to make sure we set every edge + ClothAssert(E0.minCoeff() >= 0); + ClothAssert(B0.minCoeff() >= 0); + + } // end compute edges }; } diff --git a/src/Objective.hpp b/src/Objective.hpp index 4e0d1ff..6ced5f2 100644 --- a/src/Objective.hpp +++ b/src/Objective.hpp @@ -7,100 +7,84 @@ #include "ClothAssert.hpp" #include "ClothMesh.hpp" -namespace cloth -{ +namespace cloth { class Objective { -public: - - // Computes and returns the Hookean energy of a spring with rest length r and stiffness k. - // If x.rows == grad.rows, the gradients (f, -f) are added to grad. - double spring_gradient( - const Eigen::MatrixXd &x, - int e0, - int e1, - double r, - double k, - Eigen::MatrixXd &grad) const - { - using namespace Eigen; - Vector3d x0 = x.row(e0); - Vector3d x1 = x.row(e1); - Vector3d edge = (x0-x1); - double l = edge.norm(); - if (grad.rows() == x.rows() && l > 1e-12) - { - edge /= l; - edge *= (l-r)*k; - ClothAssert(edge.allFinite()); - grad.row(e0) += edge; - grad.row(e1) -= edge; - } - double energy = 0.5 * k * (l-r) * (l-r); // (k/2)||l-l0||^2 - ClothAssert(std::isfinite(energy)); - return energy; - } - - // Spring energies that push vertices out of collision. - // We'll cheat a bit and hard code the sphere radius and center since we know it. - double collision_gradient( - const Eigen::MatrixXd &x, - int idx, - double k, - Eigen::MatrixXd &grad) const - { - using namespace Eigen; - Vector3d cent(0,0,0); - double rad = 0.51; - Vector3d xi = x.row(idx); - Vector3d dir = (xi-cent); - double l = dir.norm(); - if (l > rad) { return 0; } - if (grad.rows() == x.rows() && l > 1e-12) - { - dir /= l; - dir *= (l-rad)*k; - ClothAssert(dir.allFinite()); - grad.row(idx) += dir; - } - double energy = 0.5 * k * (l-rad) * (l-rad); // (k/2)||l-l0||^2 - ClothAssert(std::isfinite(energy)); - return energy; - } + public: + // Computes and returns the Hookean energy of a spring with rest length r and stiffness k. + // If x.rows == grad.rows, the gradients (f, -f) are added to grad. + double spring_gradient(const Eigen::MatrixXd& x, int e0, int e1, double r, double k, Eigen::MatrixXd& grad) const + { + using namespace Eigen; + Vector3d x0 = x.row(e0); + Vector3d x1 = x.row(e1); + Vector3d edge = (x0 - x1); + double l = edge.norm(); + if (grad.rows() == x.rows() && l > 1e-12) { + edge /= l; + edge *= (l - r) * k; + ClothAssert(edge.allFinite()); + grad.row(e0) += edge; + grad.row(e1) -= edge; + } + double energy = 0.5 * k * (l - r) * (l - r); // (k/2)||l-l0||^2 + ClothAssert(std::isfinite(energy)); + return energy; + } - // Momentum potential: ||x-xbar||^2_M / (2dt^2) - // Like springs above, if xrows == grad.rows, gradient is added. - // Unlike springs, the momentum is a global term. - inline double momentum_gradient( - const Eigen::MatrixXd &x, - const Eigen::MatrixXd &x_start, - const Eigen::MatrixXd &v, - const Eigen::VectorXd &masses, - double dt, - Eigen::MatrixXd &grad) const - { - using namespace Eigen; - // Some of this can be precomputed to speed up evaluations - MatrixXd x_bar = x_start + dt * v; - MatrixXd dx = x - x_bar; - MatrixXd M_dx = masses.asDiagonal() * dx; - double dot = 0; - for (int i=0; i<3; ++i) - { - dot += dx.col(i).dot(M_dx.col(i)); - } + // Spring energies that push vertices out of collision. + // We'll cheat a bit and hard code the sphere radius and center since we know it. + double collision_gradient(const Eigen::MatrixXd& x, int idx, double k, Eigen::MatrixXd& grad) const + { + using namespace Eigen; + Vector3d cent(0, 0, 0); + double rad = 0.51; + Vector3d xi = x.row(idx); + Vector3d dir = (xi - cent); + double l = dir.norm(); + if (l > rad) { + return 0; + } + if (grad.rows() == x.rows() && l > 1e-12) { + dir /= l; + dir *= (l - rad) * k; + ClothAssert(dir.allFinite()); + grad.row(idx) += dir; + } + double energy = 0.5 * k * (l - rad) * (l - rad); // (k/2)||l-l0||^2 + ClothAssert(std::isfinite(energy)); + return energy; + } - double energy = dot / (2.0*dt*dt); - if (grad.rows() == x.rows()) - { - M_dx *= (1.0 / (dt*dt)); - grad += M_dx; - } - ClothAssert(std::isfinite(energy)); - return energy; - } + // Momentum potential: ||x-xbar||^2_M / (2dt^2) + // Like springs above, if xrows == grad.rows, gradient is added. + // Unlike springs, the momentum is a global term. + inline double momentum_gradient(const Eigen::MatrixXd& x, + const Eigen::MatrixXd& x_start, + const Eigen::MatrixXd& v, + const Eigen::VectorXd& masses, + double dt, + Eigen::MatrixXd& grad) const + { + using namespace Eigen; + // Some of this can be precomputed to speed up evaluations + MatrixXd x_bar = x_start + dt * v; + MatrixXd dx = x - x_bar; + MatrixXd M_dx = masses.asDiagonal() * dx; + double dot = 0; + for (int i = 0; i < 3; ++i) { + dot += dx.col(i).dot(M_dx.col(i)); + } + double energy = dot / (2.0 * dt * dt); + if (grad.rows() == x.rows()) { + M_dx *= (1.0 / (dt * dt)); + grad += M_dx; + } + ClothAssert(std::isfinite(energy)); + return energy; + } }; } diff --git a/src/Solver.hpp b/src/Solver.hpp index 1b065a3..a959512 100644 --- a/src/Solver.hpp +++ b/src/Solver.hpp @@ -4,129 +4,119 @@ #ifndef SOLVER_HPP #define SOLVER_HPP 1 -#include "Objective.hpp" #include "ClothMesh.hpp" +#include "Objective.hpp" -namespace cloth -{ +namespace cloth { class Solver { -public: - - // State variables - Eigen::MatrixXd x; // positions [nx3] - Eigen::MatrixXd v; // velocities [nx3] - Eigen::MatrixXd x_start; // start of timestep [nx3] - double spring_k; - double bend_k; - double collision_k; - int max_solver_iter; - - Solver() : - spring_k(100), - bend_k(100), - collision_k(300), - max_solver_iter(400) - {} - - // Initalizes state variables, returns true if success - bool initialize(const ClothMesh &cloth) - { - if (cloth.V.rows() == 0 || cloth.V.cols() != 3) - { - return false; - } - - x = cloth.V; - v = Eigen::MatrixXd::Zero(x.rows(), 3); - x_start = x; - return true; - } - - // Computes the energy and gradient (if desired) of the objective - double gradient(const ClothMesh &cloth, const Objective &objective, - double dt, const Eigen::MatrixXd &xk, Eigen::MatrixXd &grad) const - { - double tot_energy = 0; - - // Momentum potential - tot_energy += objective.momentum_gradient( - xk, x_start, v, cloth.masses, dt, grad); - - // Stretch springs - int num_edges = cloth.E.rows(); - for (int i=0; i energy_k) - { - alpha *= 0.5; - x = xk - alpha * grad; - } - } - } - + public: + // State variables + Eigen::MatrixXd x; // positions [nx3] + Eigen::MatrixXd v; // velocities [nx3] + Eigen::MatrixXd x_start; // start of timestep [nx3] + double spring_k; + double bend_k; + double collision_k; + int max_solver_iter; + + Solver() + : spring_k(100) + , bend_k(100) + , collision_k(300) + , max_solver_iter(400) + { + } + + // Initalizes state variables, returns true if success + bool initialize(const ClothMesh& cloth) + { + if (cloth.V.rows() == 0 || cloth.V.cols() != 3) { + return false; + } + + x = cloth.V; + v = Eigen::MatrixXd::Zero(x.rows(), 3); + x_start = x; + return true; + } + + // Computes the energy and gradient (if desired) of the objective + double gradient(const ClothMesh& cloth, + const Objective& objective, + double dt, + const Eigen::MatrixXd& xk, + Eigen::MatrixXd& grad) const + { + double tot_energy = 0; + + // Momentum potential + tot_energy += objective.momentum_gradient(xk, x_start, v, cloth.masses, dt, grad); + + // Stretch springs + int num_edges = cloth.E.rows(); + for (int i = 0; i < num_edges; ++i) { + Eigen::RowVector2i e = cloth.E.row(i); + tot_energy += objective.spring_gradient(xk, e[0], e[1], cloth.E0[i], spring_k, grad); + } + + // Linear bend springs + int num_bend = cloth.B.rows(); + for (int i = 0; i < num_bend; ++i) { + Eigen::RowVector2i b = cloth.B.row(i); + tot_energy += objective.spring_gradient(xk, b[0], b[1], cloth.B0[i], bend_k, grad); + } + + // Collision springs + int num_verts = xk.rows(); + for (int i = 0; i < num_verts; ++i) { + tot_energy += objective.collision_gradient(xk, i, collision_k, grad); + } + + return tot_energy; + } + + double energy(const ClothMesh& cloth, const Objective& objective, double dt, const Eigen::MatrixXd& xk) const + { + Eigen::MatrixXd dummy; // grad computation ignored if grad.size != x.size + return gradient(cloth, objective, dt, xk, dummy); + } + + void solve(const ClothMesh& cloth, const Objective& objective, double dt) + { + using namespace Eigen; + ClothAssert(x.rows() == cloth.V.rows() && x.cols() == 3); + ClothAssert(v.rows() == x.rows() && v.cols() == x.cols()); + + x_start = x; // cache start state + v.col(1).array() += dt * -9.8; // gravity + x = x_start + dt * v; // initial guess + MatrixXd grad = MatrixXd::Zero(x.rows(), 3); + MatrixXd xk = MatrixXd::Zero(x.rows(), 3); + + for (int iter = 0; iter < max_solver_iter; ++iter) { + // Compute current energy and gradient + xk = x; // store state at current iteration + grad.setZero(); + double energy_k = gradient(cloth, objective, dt, x, grad); + + // Check for convergence + if (grad.norm() < 0.01) { + break; + } + + // Take a full step, i.e. alpha=1 + x = xk - grad; + + // Line search until objective decreased + double alpha = 1; + while (energy(cloth, objective, dt, x) > energy_k) { + alpha *= 0.5; + x = xk - alpha * grad; + } + } + } }; } diff --git a/test/main.cpp b/test/main.cpp index 2b04ccd..ac960ec 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -1,82 +1,77 @@ // Copyright Matt Overby 2021. // Distributed under the MIT License. +#include "Solver.hpp" +#include #include #include #include -#include -#include "Solver.hpp" using namespace Eigen; using namespace cloth; -int main(int argc, char *argv[]) +int +main(int argc, char* argv[]) { - (void)(argc); - (void)(argv); + (void)(argc); + (void)(argv); - std::string plane = CLOTH_ROOT_DIR "/data/plane.obj"; - std::string sphere = CLOTH_ROOT_DIR "/data/sphere.obj"; - MatrixXd V, sphereV; - MatrixXi F, sphereF; - if (!igl::readOBJ(plane, V, F) || !igl::readOBJ(sphere, sphereV, sphereF)) - { - return EXIT_FAILURE; - } + std::string plane = CLOTH_ROOT_DIR "/data/plane.obj"; + std::string sphere = CLOTH_ROOT_DIR "/data/sphere.obj"; + MatrixXd V, sphereV; + MatrixXi F, sphereF; + if (!igl::readOBJ(plane, V, F) || !igl::readOBJ(sphere, sphereV, sphereF)) { + return EXIT_FAILURE; + } - // Move cloth above the sphere - V.col(1).array() += 0.6; + // Move cloth above the sphere + V.col(1).array() += 0.6; - // Create cloth object - ClothMesh cloth; - cloth.add_mesh(V, F); + // Create cloth object + ClothMesh cloth; + cloth.add_mesh(V, F); - // Objective function - Objective objective; + // Objective function + Objective objective; - // Create solver - double dt = 1.0 / 24.0; - Solver solver; - if (!solver.initialize(cloth)) - { - return EXIT_FAILURE; - } + // Create solver + double dt = 1.0 / 24.0; + Solver solver; + if (!solver.initialize(cloth)) { + return EXIT_FAILURE; + } - std::cout << "Press A to simulate, R to reset" << std::endl; + std::cout << "Press A to simulate, R to reset" << std::endl; - // Create viewer - igl::opengl::glfw::Viewer viewer; - viewer.data().set_face_based(true); - viewer.data().set_mesh(V, F); - viewer.data().set_colors(RowVector3d(0.659, 0.847, 1)); - viewer.append_mesh(); - viewer.data(1).set_mesh(sphereV, sphereF); - viewer.data(1).set_colors(RowVector3d(0.7, 0.7, 0.7)); - viewer.core().is_animating = false; + // Create viewer + igl::opengl::glfw::Viewer viewer; + viewer.data().set_face_based(true); + viewer.data().set_mesh(V, F); + viewer.data().set_colors(RowVector3d(0.659, 0.847, 1)); + viewer.append_mesh(); + viewer.data(1).set_mesh(sphereV, sphereF); + viewer.data(1).set_colors(RowVector3d(0.7, 0.7, 0.7)); + viewer.core().is_animating = false; - viewer.callback_key_pressed = [&](igl::opengl::glfw::Viewer&, unsigned int key, int)->bool - { - char c = char(key); - if (c == 'r') - { - solver.initialize(cloth); - } - return false; - }; + viewer.callback_key_pressed = [&](igl::opengl::glfw::Viewer&, unsigned int key, int) -> bool { + char c = char(key); + if (c == 'r') { + solver.initialize(cloth); + } + return false; + }; - viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer&)->bool - { - if (viewer.core().is_animating) - { - solver.solve(cloth, objective, dt); // solve the time step - V = solver.x; // update vertex positions - viewer.data(0).set_mesh(V, F); - viewer.data(0).compute_normals(); // update normals after defo - } - return false; - }; + viewer.callback_pre_draw = [&](igl::opengl::glfw::Viewer&) -> bool { + if (viewer.core().is_animating) { + solver.solve(cloth, objective, dt); // solve the time step + V = solver.x; // update vertex positions + viewer.data(0).set_mesh(V, F); + viewer.data(0).compute_normals(); // update normals after defo + } + return false; + }; - viewer.launch(); + viewer.launch(); - return EXIT_SUCCESS; + return EXIT_SUCCESS; } diff --git a/test/test.cpp b/test/test.cpp index 5a6bd9a..d07197b 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -1,56 +1,54 @@ // Copyright Matt Overby 2021. // Distributed under the MIT License. +#include "Solver.hpp" +#include #include #include -#include -#include "Solver.hpp" using namespace Eigen; using namespace cloth; -int main(int argc, char *argv[]) +int +main(int argc, char* argv[]) { - (void)(argc); - (void)(argv); - - std::string plane = CLOTH_ROOT_DIR "/data/plane.obj"; - std::string sphere = CLOTH_ROOT_DIR "/data/sphere.obj"; - MatrixXd V, cV; - MatrixXi F, cF; - if (!igl::readOBJ(plane, V, F) || !igl::readOBJ(sphere, cV, cF)) - { - return EXIT_FAILURE; - } - - // Move cloth above the sphere - V.col(1).array() += 0.5; - - // Create cloth object - ClothMesh cloth; - cloth.add_mesh(V, F); - - // Objective function - Objective objective; - - // Create solver - double dt = 1.0 / 24.0; - Solver solver; - - // 20 iters won't converge but it will test - // all the functions are operational. - solver.max_solver_iter = 20; - - if (!solver.initialize(cloth)) - { - return EXIT_FAILURE; - } + (void)(argc); + (void)(argv); + + std::string plane = CLOTH_ROOT_DIR "/data/plane.obj"; + std::string sphere = CLOTH_ROOT_DIR "/data/sphere.obj"; + MatrixXd V, cV; + MatrixXi F, cF; + if (!igl::readOBJ(plane, V, F) || !igl::readOBJ(sphere, cV, cF)) { + return EXIT_FAILURE; + } + + // Move cloth above the sphere + V.col(1).array() += 0.5; + + // Create cloth object + ClothMesh cloth; + cloth.add_mesh(V, F); + + // Objective function + Objective objective; + + // Create solver + double dt = 1.0 / 24.0; + Solver solver; + + // 20 iters won't converge but it will test + // all the functions are operational. + solver.max_solver_iter = 20; + + if (!solver.initialize(cloth)) { + return EXIT_FAILURE; + } int timesteps = 10; - for (int step=0; step < timesteps; ++step) - { + for (int step = 0; step < timesteps; ++step) { solver.solve(cloth, objective, dt); // solve the time step } - return EXIT_SUCCESS; + return EXIT_SUCCESS; }