diff --git a/scripts/test_cpp.py b/scripts/test_cpp.py
new file mode 100644
index 00000000..e4f13344
--- /dev/null
+++ b/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()
diff --git a/src/compas_fd/fd/__init__.py b/src/compas_fd/fd/__init__.py
index d9b2165d..01b206d1 100644
--- a/src/compas_fd/fd/__init__.py
+++ b/src/compas_fd/fd/__init__.py
@@ -15,6 +15,7 @@
 
     fd_numpy
     mesh_fd_numpy
+    fd_cpp
 
 """
 from __future__ import print_function
@@ -26,6 +27,7 @@
 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__ = []
 
@@ -33,4 +35,5 @@
     __all__ += [
         'fd_numpy',
         'mesh_fd_numpy',
+        'fd_cpp'
     ]
diff --git a/src/compas_fd/fd/fd_cpp.py b/src/compas_fd/fd/fd_cpp.py
deleted file mode 100644
index 19fe1556..00000000
--- a/src/compas_fd/fd/fd_cpp.py
+++ /dev/null
@@ -1 +0,0 @@
-# should be compiled with pybind11
diff --git a/src/compas_fd/fd/fd_cpp/__init__.py b/src/compas_fd/fd/fd_cpp/__init__.py
new file mode 100644
index 00000000..bd2f9de4
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/__init__.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']
diff --git a/src/compas_fd/fd/fd_cpp/docs/PLACEHOLDER b/src/compas_fd/fd/fd_cpp/docs/PLACEHOLDER
new file mode 100644
index 00000000..e69de29b
diff --git a/src/compas_fd/fd/fd_cpp/fd_cpp.py b/src/compas_fd/fd/fd_cpp/fd_cpp.py
new file mode 100644
index 00000000..89d09c8d
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/fd_cpp.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))
diff --git a/src/compas_fd/fd/fd_cpp/include/compas.h b/src/compas_fd/fd/fd_cpp/include/compas.h
new file mode 100644
index 00000000..141032ec
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/include/compas.h
@@ -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
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.cpp b/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.cpp
new file mode 100644
index 00000000..5fd1b5eb
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.cpp
@@ -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;
+}
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.h b/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.h
new file mode 100644
index 00000000..845231e2
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/connectivity_matrices.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
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/fd_solvers.cpp b/src/compas_fd/fd/fd_cpp/src/fd_solvers.cpp
new file mode 100644
index 00000000..5920d337
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/fd_solvers.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());
+};
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/fd_solvers.h b/src/compas_fd/fd/fd_cpp/src/fd_solvers.h
new file mode 100644
index 00000000..1d75a1da
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/fd_solvers.h
@@ -0,0 +1,14 @@
+#ifndef FD_SOLVERS
+#define FD_SOLVERS
+
+
+// One-shot equilibrium calculation by the force density method.
+std::tuple<Md3, Md3, Md1, Md1>
+fd_solve(
+	vvd& vertex_coordinates,
+	vi& fixed_vertices,
+	vvi& edges,
+	vd& force_densities,
+	vvd& loads);
+
+#endif  // FD_SOLVERS
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/fdm.cpp b/src/compas_fd/fd/fd_cpp/src/fdm.cpp
new file mode 100644
index 00000000..0aa4bdbe
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/fdm.cpp
@@ -0,0 +1,11 @@
+#include "compas.h"
+#include "fd_solvers.h"
+
+
+void init_fd_solvers(py::module&);
+
+PYBIND11_MODULE(fdm, m)
+{
+    m.doc() = "force density algorithms using Eigen.";
+    init_fd_solvers(m);
+}
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/linalg_conversion.h b/src/compas_fd/fd/fd_cpp/src/linalg_conversion.h
new file mode 100644
index 00000000..836a8148
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/linalg_conversion.h
@@ -0,0 +1,56 @@
+#ifndef LINALG_CONVERSION
+#define LINALG_CONVERSION
+
+#include "compas.h"
+
+
+// Construct a generic Eigen::Vector from an std::vector.
+// The same data memory address is referenced, no copy is made.
+template <typename T> Eigen::Map<Eigen::Vector<double, Eigen::Dynamic>>&
+matrix_from_vec1d(
+    std::vector<T>& vec)
+{
+    Eigen::Map<Eigen::Vector<T, Eigen::Dynamic>> out_vec(vec.data(), vec.size());
+    return out_vec;
+}
+
+
+// Construct a generic row-major Eigen::MatrixX3 from a nested std::vector.
+// A new contiguous slot of memory is allocated for the Eigen::Matrix.
+template <typename T> Eigen::Matrix<T, Eigen::Dynamic, 3>
+matrixX3_from_vec2d(
+    const std::vector<std::vector<T>>& vec)
+{
+    int row_count = vec.size();
+    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> out_mat(vec.size(), 3);
+    for (unsigned int i = 0; i < vec.size(); ++i)
+        out_mat.row(i) = Eigen::Vector<T, Eigen::Dynamic>::Map(vec[i].data(), 3);
+    return out_mat;
+}
+
+
+// Construct a generic std::vector from an Eigen::Vector.
+// A new slot of memory is allocated for the std::vector.
+template <typename T> std::vector<T>
+matrix_to_vec1d(
+    const Eigen::Vector<T, Eigen::Dynamic>& mat)
+{
+    return std::vector<T>(&mat[0], mat.data() + mat.cols() * mat.rows());
+}
+
+
+// Construct a nested std::vector from an Eigen::Matrix of 3 columns.
+// A new slot of memory is allocated for the std::vector.
+template <typename T> std::vector<std::array<T, 3>>
+matrixX3_to_vec2d(
+    const Eigen::Matrix<T, Eigen::Dynamic, 3>& mat)
+{
+    std::vector<std::array<T, 3>> out_vec;
+    out_vec.resize(mat.rows());
+    Eigen::Matrix<T, Eigen::Dynamic, 3, Eigen::RowMajor>::Map(
+        out_vec[0].data(), mat.rows(), mat.cols()) = mat;
+    return out_vec;
+}
+
+
+#endif  // LINALG_CONVERSION
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/process_vertices.cpp b/src/compas_fd/fd/fd_cpp/src/process_vertices.cpp
new file mode 100644
index 00000000..f8b76ee5
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/process_vertices.cpp
@@ -0,0 +1,18 @@
+#include "compas.h"
+#include "process_vertices.h"
+
+
+void
+set_free_vertices(
+    const int& vertex_count,
+    const int& edge_count,
+    const vi& fixed_vertices,
+    vi& free_vertices)
+{
+    vi all_vertices(vertex_count);
+    std::iota(all_vertices.begin(), all_vertices.end(), 0);
+    free_vertices.reserve(vertex_count - fixed_vertices.size());
+    std::set_difference(all_vertices.begin(), all_vertices.end(),
+        fixed_vertices.begin(), fixed_vertices.end(),
+        std::inserter(free_vertices, free_vertices.begin()));
+}
\ No newline at end of file
diff --git a/src/compas_fd/fd/fd_cpp/src/process_vertices.h b/src/compas_fd/fd/fd_cpp/src/process_vertices.h
new file mode 100644
index 00000000..dff4a35a
--- /dev/null
+++ b/src/compas_fd/fd/fd_cpp/src/process_vertices.h
@@ -0,0 +1,17 @@
+#ifndef PROCESS_VERTICES
+#define PROCESS_VERTICES
+
+#include "compas.h"
+
+
+// Set the indices of all free vertices.
+// free vertices are defined in-place as output parameters.
+void
+set_free_vertices(
+    const int& vertex_count,
+    const int& edge_count,
+    const vi& fixed_vertices,
+    vi& free_vertices);
+
+
+#endif  // PROCESS_VERTICES
\ No newline at end of file