Add c++ solver

Force density solver implemented in C++ through pybind11. The solver gets a speed bump, especially for medium-sized networks (5-20x).

Exposes fd_solve() in the fdm.pyd library. Current built for win64. A makefile for other platforms still needs to be included.

Author: Sam-Bouten

scripts/test_cpp.py

@@ -0,0 +1,82 @@
+import compas
+from compas.datastructures import mesh_subdivide
+from compas_fd.datastructures import CableMesh
+from compas_fd.fd import fd_cpp, fd_numpy
+from compas_view2.app import App
+from compas_view2.objects import Object, MeshObject
+from timeit import timeit
+
+Object.register(CableMesh, MeshObject)
+
+
+# =================================================
+# input mesh
+# =================================================
+
+mesh = CableMesh.from_obj(compas.get('faces.obj'))
+
+mesh.vertices_attribute('is_anchor', True, mesh.vertices_where({'vertex_degree': 2}))
+dva = {'rx': .0, 'ry': .0, 'rz': .0,
+       'px': .0, 'py': .0, 'pz': -.0,
+       'is_anchor': False}
+
+mesh = mesh_subdivide(mesh, scheme='quad', k=1)
+mesh.update_default_vertex_attributes(dva)
+
+
+# =================================================
+# pre-process mesh data
+# =================================================
+
+vertices = [mesh.vertex_coordinates(v) for v in mesh.vertices()]
+fixed = list(mesh.vertices_where({'is_anchor': True}))
+edges = list(mesh.edges())
+force_densities = mesh.edges_attribute('q')
+loads = mesh.vertices_attributes(['px', 'py', 'pz'])
+
+
+# =================================================
+# solvers
+# =================================================
+
+result_eigen = fd_cpp(vertices=vertices, fixed=fixed, edges=edges,
+                      force_densities=force_densities, loads=loads)
+
+result_np = fd_numpy(vertices=vertices, fixed=fixed, edges=edges,
+                     forcedensities=force_densities, loads=loads)
+
+forces_comparison = ((abs(l1 - l2) < 1E-14) for l1, l2
+                     in zip(result_eigen.forces, result_np.forces))
+print("Check if forces are equal between solvers: ", all(forces_comparison))
+
+
+for vertex, coo in zip(mesh.vertices(), result_eigen.vertices):
+    mesh.vertex_attributes(vertex, 'xyz', coo)
+
+
+# =================================================
+# benchmarks
+# =================================================
+
+def fn_cpp():
+    fd_cpp(vertices=vertices, fixed=fixed, edges=edges,
+           force_densities=force_densities, loads=loads)
+
+
+def fn_numpy():
+    fd_numpy(vertices=vertices, fixed=fixed, edges=edges,
+             forcedensities=force_densities, loads=loads)
+
+
+nr = 30
+print(f"Solver time in cpp:   {round(1000 * timeit(fn_cpp, number=nr) / nr, 2)} ms")
+print(f"Solver time in numpy: {round(1000 * timeit(fn_numpy, number=nr) / nr, 2)} ms")
+
+
+# =================================================
+# viz
+# =================================================
+
+viewer = App()
+viewer.add(mesh)
+viewer.run()

src/compas_fd/fd/__init__.py

@@ -15,6 +15,7 @@
 
     fd_numpy
     mesh_fd_numpy
+    fd_cpp
 
 """
 from __future__ import print_function
@@ -26,11 +27,13 @@
 if not compas.IPY:
     from .fd_numpy import fd_numpy
     from .mesh_fd_numpy import mesh_fd_numpy
+    from .fd_cpp import fd_cpp
 
 __all__ = []
 
 if not compas.IPY:
     __all__ += [
         'fd_numpy',
         'mesh_fd_numpy',
+        'fd_cpp'
     ]

src/compas_fd/fd/fd_cpp.py

@@ -0,0 +1,31 @@
+"""
+******************
+compas_fd.fd.fd_cpp
+******************
+
+.. currentmodule:: compas_fd.fd.fd_cpp
+
+
+Functions
+=========
+
+.. autosummary::
+    :toctree: generated/
+    :nosignatures:
+
+    fd_cpp
+
+"""
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+
+import compas
+
+if not compas.IPY:
+    from .fd_cpp import fd_cpp
+
+__all__ = []
+
+if not compas.IPY:
+    __all__ += ['fd_cpp']

src/compas_fd/fd/fd_cpp/__init__.py

@@ -0,0 +1,17 @@
+from typing import Optional, Sequence
+from typing_extensions import Annotated
+
+from ..result import Result
+from .fdm import fd_solve   # win64 build
+
+
+def fd_cpp(*,
+           vertices: Sequence[Annotated[Sequence[float], 3]],
+           fixed: Sequence[int],
+           edges: Sequence[Annotated[Sequence[int], 2]],
+           force_densities: Sequence[float],
+           loads: Optional[Sequence[Annotated[Sequence[float], 3]]] = None
+           ) -> Result:
+    """Compute the equilibrium conditions by the force density method.
+    The algorithms backend is implemented in C++ through pybind."""
+    return Result(*fd_solve(vertices, fixed, edges, force_densities, loads))

src/compas_fd/fd/fd_cpp/docs/PLACEHOLDER

@@ -0,0 +1,29 @@
+#ifndef COMPAS
+#define COMPAS
+
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <pybind11/stl.h>
+
+#include <vector>
+#include <tuple>
+#include <numeric>
+#include <exception>
+#include <Eigen/Sparse>
+#include <Eigen/SparseCholesky>
+
+
+namespace py = pybind11;
+
+using vi = std::vector<int>;
+using vd = std::vector<double>;
+using vvi = std::vector<std::vector<int>>;
+using vvd = std::vector<std::vector<double>>;
+using Md = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+using Md3 = Eigen::Matrix<double, Eigen::Dynamic, 3>;
+using Md1 = Eigen::Vector<double, Eigen::Dynamic>;
+using DMd = Eigen::DiagonalMatrix<double, Eigen::Dynamic, Eigen::Dynamic>;
+using SMd = Eigen::SparseMatrix<double>;
+
+
+#endif // COMPAS

src/compas_fd/fd/fd_cpp/fd_cpp.py

@@ -0,0 +1,41 @@
+#include "compas.h"
+#include "connectivity_matrices.h"
+
+
+//aliases
+using T = Eigen::Triplet<int>;
+
+void
+set_connectivity_matrices(
+    const vvi& edges,
+    const vi& free_vertices,
+    const vi& fixed_vertices,
+    SMd& C, SMd& Ci, SMd& Cf)
+{
+    std::vector<T> triplets, free_triplets, fixed_triplets;
+    triplets.reserve(2 * edges.size());
+    free_triplets.reserve(free_vertices.size());
+    fixed_triplets.reserve(fixed_vertices.size());
+    SMd Ci_mask(C.cols(), free_vertices.size());
+    SMd Cf_mask(C.cols(), fixed_vertices.size());
+
+    // full connectivity matrix
+    for (unsigned int i = 0; i < edges.size(); ++i)
+    {
+        triplets.emplace_back(i, edges[i][0],  1);
+        triplets.emplace_back(i, edges[i][1], -1);
+    }
+    C.setFromTriplets(triplets.begin(), triplets.end());
+
+    // free vertices connectivity matrix
+    for (unsigned int i = 0; i < free_vertices.size(); ++i)
+        free_triplets.emplace_back(free_vertices[i], i, 1);
+    Ci_mask.setFromTriplets(free_triplets.begin(), free_triplets.end());
+    Ci = C * Ci_mask;
+
+    // fixed vertices connectivity matrix
+    for (unsigned int i = 0; i < fixed_vertices.size(); ++i)
+        fixed_triplets.emplace_back(fixed_vertices[i], i, 1);
+    Cf_mask.setFromTriplets(fixed_triplets.begin(), fixed_triplets.end());
+    Cf = C * Cf_mask;
+}

src/compas_fd/fd/fd_cpp/include/compas.h

@@ -0,0 +1,20 @@
+#ifndef CONNECTIVITY_MATRICES
+#define CONNECTIVITY_MATRICES
+
+#include "compas.h"
+
+
+// Construct the connectivity matrices of a connected graph.
+// Coefficient (i, j) equals 1 if edge i starts at vertex i,
+// -1 if edge i ends at vertex j, and 0 otherwise.
+// Matrices are defined in-place as output parameters:
+// full matrix C, free vertices matrix Ci, fixed vertices matrix Cf.
+void
+set_connectivity_matrices(
+    const vvi& edges,
+    const vi& free_vertices,
+    const vi& fixed_vertices,
+    SMd& C, SMd& Ci, SMd& Cf);
+
+
+#endif  // CONNECTIVITY_MATRICES

src/compas_fd/fd/fd_cpp/src/connectivity_matrices.cpp

@@ -0,0 +1,70 @@
+#include "compas.h"
+#include "fd_solvers.h"
+#include "process_vertices.h"
+#include "connectivity_matrices.h"
+#include "linalg_conversion.h"
+
+
+std::tuple<Md3, Md3, Md1, Md1>
+fd_solve(
+    vvd& vertex_coordinates,
+    vi& fixed_vertices,
+    vvi& edges,
+    vd& force_densities,
+    vvd& loads)
+{
+    // pre-process vertex indices arrays
+    int vertex_count = vertex_coordinates.size();
+    int edge_count = edges.size();
+    vi free_vertices;
+    set_free_vertices(vertex_count, edge_count,
+                      fixed_vertices, free_vertices);
+
+    // set primary data matrices
+    Md3 X = matrixX3_from_vec2d(vertex_coordinates);
+    Md3 P = matrixX3_from_vec2d(loads);
+    auto& q = matrix_from_vec1d(force_densities);
+    DMd Q = q.asDiagonal();
+
+    // set connectivity matrices
+    SMd C(edge_count, vertex_count);
+    SMd Ci(edge_count, free_vertices.size());
+    SMd Cf(edge_count, fixed_vertices.size());
+    set_connectivity_matrices(edges, free_vertices, fixed_vertices, C, Ci, Cf);
+
+    // solve for current force densities
+    // build stiffness matrices
+    SMd Cit = Ci.transpose();
+    SMd D = C.transpose() * Q * C;
+    SMd Di = Cit * Q * Ci;
+    SMd Df = Cit * Q * Cf;
+
+    // solve for free coordinates
+    Md b = P(free_vertices, Eigen::all) - (Df * X(fixed_vertices, Eigen::all));
+    Eigen::SimplicialLDLT<SMd> solver;
+    solver.compute(Di);
+    if (solver.info() != Eigen::Success)
+        throw std::exception("Singular stiffness matrix");
+    Md3 X_free = solver.solve(b);                   // updated free vertex coordinates
+    for (unsigned int i = 0; i < X_free.rows(); ++i)
+        X.row(free_vertices[i]) = X_free.row(i);
+
+    // get dependent variables
+    Md3 R = P - (D * X);                            // residuals and reactions
+    Md1 L = (C * X).rowwise().norm();               // edge lengths
+    Md1 F = Q.diagonal().array() * L.array();       // edge forces
+
+    return {X, R, F, L};
+}
+
+
+void init_fd_solvers(py::module& m)
+{
+    m.def("fd_solve",
+          &fd_solve,
+          py::arg("vertex_coordinates").noconvert(),
+          py::arg("fixed_vertices").noconvert(),
+          py::arg("edges").noconvert(),
+          py::arg("force_densities").noconvert(),
+          py::arg("loads").noconvert());
+};