diff --git a/src/compas_fd/fd/fd_numerical_data.py b/src/compas_fd/fd/fd_numerical_data.py
deleted file mode 100644
index fb061df5..00000000
--- a/src/compas_fd/fd/fd_numerical_data.py
+++ /dev/null
@@ -1,74 +0,0 @@
-from typing import Any
-from typing import Tuple
-from typing import List
-from typing import Union
-from typing import Sequence
-from typing import Optional
-from typing_extensions import Annotated
-from nptyping import NDArray
-
-from dataclasses import dataclass
-from dataclasses import astuple
-
-from numpy import asarray
-from numpy import float64
-from numpy import zeros_like
-from scipy.sparse import diags
-
-from compas.numerical import connectivity_matrix
-
-from .result import Result
-
-
-@dataclass
-class FDNumericalData:
-    """Stores numerical data used by the force density algorithms."""
-    free: int
-    fixed: int
-    xyz: NDArray[(Any, 3), float64]
-    C: NDArray[(Any, Any), int]
-    q: NDArray[(Any, 1), float64]
-    Q: NDArray[(Any, Any), float64]
-    p: NDArray[(Any, 1), float64]
-    A: NDArray[(Any, Any), float64]
-    Ai: NDArray[(Any, Any), float64]
-    Af: NDArray[(Any, Any), float64]
-    forces: NDArray[(Any, 1), float64] = None
-    lengths: NDArray[(Any, 1), float64] = None
-    residuals: NDArray[(Any, 3), float64] = None
-    tangent_residuals: NDArray[(Any, 3), float64] = None
-    normal_residuals: NDArray[(Any, 1), float64] = None
-
-    def __iter__(self):
-        return iter(astuple(self))
-
-    @classmethod
-    def from_params(cls,
-                    vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]],
-                    fixed: List[int],
-                    edges: List[Tuple[int, int]],
-                    forcedensities: List[float],
-                    loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None):
-        """Construct numerical arrays from force density solver input parameters."""
-        free = list(set(range(len(vertices))) - set(fixed))
-        xyz = asarray(vertices, dtype=float64).reshape((-1, 3))
-        C = connectivity_matrix(edges, 'csr')
-        Ci = C[:, free]
-        Cf = C[:, fixed]
-        q = asarray(forcedensities, dtype=float64).reshape((-1, 1))
-        Q = diags([q.flatten()], [0])
-        p = (zeros_like(xyz) if loads is None else
-             asarray(loads, dtype=float64).reshape((-1, 3)))
-        A = C.T.dot(Q).dot(C)
-        Ai = Ci.T.dot(Q).dot(Ci)
-        Af = Ci.T.dot(Q).dot(Cf)
-        return cls(free, fixed, xyz, C, q, Q, p, A, Ai, Af)
-
-    @classmethod
-    def from_mesh(cls, mesh):
-        """Construct numerical arrays from input mesh."""
-        raise NotImplementedError
-
-    def to_result(self) -> Result:
-        """Parse relevant numerical data into a Result object."""
-        return Result(self.xyz, self.residuals, self.forces, self.lengths)
diff --git a/src/compas_fd/numdata/__init__.py b/src/compas_fd/numdata/__init__.py
new file mode 100644
index 00000000..036c517b
--- /dev/null
+++ b/src/compas_fd/numdata/__init__.py
@@ -0,0 +1,33 @@
+"""
+******************
+compas_fd.numdata
+******************
+
+.. currentmodule:: compas_fd.numdata
+
+
+Classes
+=======
+
+.. autosummary::
+    :toctree: generated/
+    :nosignatures:
+
+    NumericalData
+    FDNumericalData
+
+"""
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+
+import compas
+from .numerical_data import NumericalData
+
+if not compas.IPY:
+    from .fd_numerical_data import FDNumericalData
+
+__all__ = ['NumericalData']
+
+if not compas.IPY:
+    __all__ += ['FDNumericalData']
diff --git a/src/compas_fd/numdata/face_data.py b/src/compas_fd/numdata/face_data.py
new file mode 100644
index 00000000..d5135cd5
--- /dev/null
+++ b/src/compas_fd/numdata/face_data.py
@@ -0,0 +1,116 @@
+from typing import Any
+from typing import List
+from typing import Tuple
+from nptyping import NDArray
+
+from numpy import add
+from numpy import cross
+from numpy import float64
+from numpy import roll
+from numpy import zeros
+from numpy.linalg import norm
+from scipy.sparse import coo_matrix
+
+from compas.datastructures import Mesh
+from compas.numerical import face_matrix
+from .numerical_data import lazy_eval
+
+
+class FaceDataMixin:
+    """Stores numerical mesh face data.
+    Use as mixin for NumericalData classes.
+    """
+
+    def __init__(self, mesh: Mesh) -> None:
+        self.faces_vertices = mesh
+        self.reset_face_data()
+
+    def reset_face_data(self) -> None:
+        self._tributary_area_matrix = None
+        self._face_areas = None
+        self._face_normals = None
+        self._face_centroids = None
+        self._face_cross_products = None
+
+    @property
+    def faces_vertices(self) -> List[Tuple[int]]:
+        return self._faces_vertices
+
+    @faces_vertices.setter
+    def faces_vertices(self, mesh: Mesh) -> None:
+        v_i = mesh.vertex_index()
+        f_i = {fkey: index for index, fkey in enumerate(mesh.faces())}
+        faces_vertices = [None] * mesh.number_of_faces()
+        for f in mesh.faces():
+            faces_vertices[f_i[f]] = [v_i[v] for v in mesh.face_vertices(f)]
+        self._faces_vertices = faces_vertices
+        self._face_matrix = None
+
+    @property
+    @lazy_eval
+    def face_matrix(self) -> NDArray[(Any, Any), float64]:
+        return face_matrix(self.faces_vertices, rtype='csr', normalize=True)
+
+    @property
+    @lazy_eval
+    def tributary_area_matrix(self) -> NDArray[(Any, Any), float64]:
+        """Tributary areas matrix as a sparse (vertices x faces) array.
+        Entry a_ij holds tributary area of face j for vertex i.
+        """
+        v_count = self.xyz.shape[0]
+        f_count = len(self.faces_vertices)
+        data = []
+        rows = []
+        cols = []
+
+        for face, fcp in enumerate(self.face_cross_products):
+            face_vertices = self.faces_vertices[face]
+            partial_areas = norm(fcp, axis=1)
+            tributary_areas = partial_areas + roll(partial_areas, -1)
+            data.extend(tributary_areas)
+            rows.extend(face_vertices)
+            cols.extend([face] * len(face_vertices))
+
+        return coo_matrix((data, (rows, cols)), (v_count, f_count)).tocsr() * 0.25
+
+    @property
+    @lazy_eval
+    def face_areas(self) -> NDArray[(Any, 1), float64]:
+        return self.tributary_area_matrix.sum(axis=0).reshape((-1, 1))
+
+    @property
+    @lazy_eval
+    def face_normals(self) -> NDArray[(Any, 1), float64]:
+        f_count = len(self.faces_vertices)
+        face_normals = zeros((f_count, 3), dtype=float)
+        for face, fcp in enumerate(self.face_cross_products):
+            face_normals[face, :] = add.reduce(fcp)
+        face_normals /= norm(face_normals, axis=1)[:, None]
+        return face_normals
+
+    @property
+    @lazy_eval
+    def face_centroids(self) -> NDArray[(Any, 1), float64]:
+        return self.F.dot(self.xyz)
+
+    @property
+    @lazy_eval
+    def face_cross_products(self) -> List[NDArray[(Any, 3), float64]]:
+        """Cross products of the consecutive (centroid -> vertex)
+        vectors for each of the vertices of the face.
+        """
+        face_cross_products = []
+        centroids = self.face_centroids
+        for face, face_vertices in enumerate(self.faces_vertices):
+            vecs = self.xyz[face_vertices] - centroids[face]
+            vecs_shift = roll(vecs, -1, axis=0)
+            face_cross_products.append(cross(vecs, vecs_shift))
+        return face_cross_products
+
+    # aliases
+    F = face_matrix
+    TA = tributary_area_matrix
+    fa = face_areas
+    fn = face_normals
+    fc = face_centroids
+    fcp = face_cross_products
diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py
new file mode 100644
index 00000000..58c4d63f
--- /dev/null
+++ b/src/compas_fd/numdata/fd_numerical_data.py
@@ -0,0 +1,210 @@
+from typing import Any
+from typing import Tuple
+from typing import List
+from typing import Union
+from typing import Sequence
+from typing import Optional
+from typing_extensions import Annotated
+from nptyping import NDArray
+
+from numpy import asarray
+from numpy import float64
+from numpy import zeros_like
+from scipy.sparse import diags
+from scipy.sparse import vstack as svstack
+
+from compas.datastructures import Mesh
+from compas.numerical import connectivity_matrix
+from compas.numerical import normrow
+
+from compas_fd.result import Result
+from .numerical_data import NumericalData
+from .numerical_data import lazy_eval
+from .face_data import FaceDataMixin
+
+
+class FDNumericalData(FaceDataMixin, NumericalData):
+    """Stores numerical data used by the force density algorithms."""
+
+    def __init__(self,
+                 vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]],
+                 fixed: List[int],
+                 edges: List[Tuple[int, int]],
+                 force_densities: List[float],
+                 loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None
+                 ):
+        self.xyz = vertices
+        self.fixed = fixed
+        self.edges = edges
+        self.force_densities = force_densities
+        self.loads = loads
+
+    @classmethod
+    def from_mesh(cls, mesh: Mesh):
+        """Construct numerical arrays from input mesh."""
+        v_i = mesh.vertex_index()
+        vertices = asarray(mesh.vertices_attributes('xyz'))
+        fixed = [v_i[v] for v in mesh.vertices_where({'is_anchor': True})]
+        edges = [(v_i[u], v_i[v]) for u, v in mesh.edges_where({'_is_edge': True})]
+        force_densities = asarray([attr['q'] for key, attr in
+                                   mesh.edges_where({'_is_edge': True}, True)])
+        loads = asarray(mesh.vertices_attributes(('px', 'py', 'pz')))
+
+        numdata = cls(vertices, fixed, edges, force_densities, loads)
+        super(FDNumericalData, numdata).__init__(mesh)
+        numdata.has_face_data = True
+        return numdata
+
+    def to_result(self) -> Result:
+        """Parse relevant numerical data into a Result object."""
+        return Result(self.xyz, self.residuals, self.forces, self.lengths)
+
+    @property
+    def xyz(self) -> NDArray[(Any, 3), float64]:
+        return self._xyz
+
+    @xyz.setter
+    def xyz(self, vertices) -> None:
+        self._xyz = asarray(vertices, dtype=float64).reshape((-1, 3))
+        self.reset_xyz()
+
+    def reset_xyz(self) -> None:
+        self._lengths = None
+        self._residuals = None
+        self._coordinate_difference_matrices = None
+        self._equilibrium_matrix = None
+        if getattr(self, 'has_face_data', False):
+            FaceDataMixin.reset_face_data(self)
+
+    @property
+    def fixed(self) -> List[int]:
+        return self._fixed
+
+    @fixed.setter
+    def fixed(self, fixed) -> None:
+        self._fixed = fixed
+        self._free = list(set(range(len(self.xyz))) - set(fixed))
+
+    @property
+    def free(self) -> List[int]:
+        return self._free
+
+    @property
+    def edges(self) -> List[Tuple[int, int]]:
+        return self._edges
+
+    @edges.setter
+    def edges(self, edges: List[List[int]]) -> None:
+        self._edges = edges
+        self._connectivity_matrix = None
+        self._connectivity_matrix_free = None
+        self._connectivity_matrix_fixed = None
+
+    @property
+    @lazy_eval
+    def connectivity_matrix(self) -> NDArray[(Any, Any), float64]:
+        return connectivity_matrix(self.edges, 'csr')
+
+    @property
+    @lazy_eval
+    def connectivity_matrix_free(self) -> NDArray[(Any, Any), float64]:
+        return self.connectivity_matrix[:, self.free]
+
+    @property
+    @lazy_eval
+    def connectivity_matrix_fixed(self) -> NDArray[(Any, Any), float64]:
+        return self.connectivity_matrix[:, self.fixed]
+
+    @property
+    @lazy_eval
+    def coordinate_difference_matrices(self) -> Tuple[NDArray[(Any, Any), float64]]:
+        uvw = self.C.dot(self.xyz)
+        U = diags([uvw[:, 0].flatten()], [0])
+        V = diags([uvw[:, 1].flatten()], [0])
+        W = diags([uvw[:, 2].flatten()], [0])
+        return (U, V, W)
+
+    @property
+    @lazy_eval
+    def equilibrium_matrix(self) -> NDArray[(Any, Any), float64]:
+        U, V, W = self.coordinate_difference_matrices
+        return svstack((self.C.T.dot(U), self.C.T.dot(V), self.C.T.dot(W)))
+
+    @property
+    def force_densities(self) -> NDArray[(Any, 1), float64]:
+        return self._force_densities
+
+    @force_densities.setter
+    def force_densities(self,
+                        force_densities: Union[List[float], NDArray[(Any, 3), float64]]
+                        ) -> None:
+        self._force_densities = asarray(force_densities, dtype=float64).reshape((-1, 1))
+        self._force_densities_matrix = None
+        self._stifness_matrix = None
+        self._stifness_matrix_free = None
+        self._stifness_matrix_fixed = None
+        self._forces = None
+        self._residuals = None
+
+    @property
+    @lazy_eval
+    def force_densities_matrix(self) -> NDArray[(Any, Any), float64]:
+        return diags([self.force_densities.flatten()], [0])
+
+    @property
+    def loads(self) -> NDArray[(Any, 3), float64]:
+        return self._loads
+
+    @loads.setter
+    def loads(self,
+              loads: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]
+              ) -> None:
+        self._loads = (zeros_like(self.xyz) if loads is None else
+                       asarray(loads, dtype=float64).reshape((-1, 3)))
+        self._residuals = None
+
+    @property
+    @lazy_eval
+    def stifness_matrix(self) -> NDArray[(Any, Any), float64]:
+        return self.C.T.dot(self.Q).dot(self.C)
+
+    @property
+    @lazy_eval
+    def stifness_matrix_free(self) -> NDArray[(Any, Any), float64]:
+        return self.Ci.T.dot(self.Q).dot(self.Ci)
+
+    @property
+    @lazy_eval
+    def stifness_matrix_fixed(self) -> NDArray[(Any, Any), float64]:
+        return self.Ci.T.dot(self.Q).dot(self.Cf)
+
+    @property
+    @lazy_eval
+    def forces(self) -> NDArray[(Any, 1), float64]:
+        return self.q * self.ls
+
+    @property
+    @lazy_eval
+    def lengths(self) -> NDArray[(Any, 1), float64]:
+        return normrow(self.C.dot(self.xyz))
+
+    @property
+    @lazy_eval
+    def residuals(self) -> NDArray[(Any, 3), float64]:
+        return self.p - self.D.dot(self.xyz)
+
+    # aliases
+    C = connectivity_matrix
+    Ci = connectivity_matrix_free
+    Cf = connectivity_matrix_fixed
+    UVW = coordinate_difference_matrices
+    E = equilibrium_matrix
+    q = force_densities
+    Q = force_densities_matrix
+    p = loads
+    D = stifness_matrix
+    Di = stifness_matrix_free
+    Df = stifness_matrix_fixed
+    f = forces
+    ls = lengths
+    r = residuals
diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py
new file mode 100644
index 00000000..71499790
--- /dev/null
+++ b/src/compas_fd/numdata/numerical_data.py
@@ -0,0 +1,39 @@
+from typing import Callable
+from functools import wraps
+
+from compas_fd.result import Result
+
+
+class NumericalData:
+    """Storage class for numerical arrays used by solver algorithms."""
+
+    @classmethod
+    def from_mesh(cls, mesh):
+        """Construct numerical arrays from input mesh."""
+        raise NotImplementedError
+
+    def to_result(self) -> Result:
+        """Parse relevant numerical data into a Result object."""
+        raise NotImplementedError
+
+
+def lazy_eval(getter: Callable) -> Callable:
+    """Decorator for lazy evaluation on property getters.
+
+    A private attribute name is assumed as ('_' + getter name).
+    The computation in the getter function is done only if:
+        - the private attribute does not exist in the instance dictionary.
+        - or the private attribute is None.
+    After computation by the getter function, the private attribute is stored
+    in the instance attribute dictionary.
+    """
+    private_attr = '_' + getter.__name__
+
+    @wraps(getter)
+    def wrapper(self, *args, **kwargs):
+        attr = getattr(self, private_attr, None)
+        if attr is None:
+            attr = getter(self, *args, **kwargs)
+            self.__dict__[private_attr] = attr
+        return attr
+    return wrapper
diff --git a/src/compas_fd/result/__init__.py b/src/compas_fd/result/__init__.py
new file mode 100644
index 00000000..fa376bd3
--- /dev/null
+++ b/src/compas_fd/result/__init__.py
@@ -0,0 +1,25 @@
+"""
+******************
+compas_fd.result
+******************
+
+.. currentmodule:: compas_fd.result
+
+
+Classes
+=======
+
+.. autosummary::
+    :toctree: generated/
+    :nosignatures:
+
+    Result
+
+"""
+from __future__ import print_function
+from __future__ import absolute_import
+from __future__ import division
+
+from .result import Result
+
+__all__ = ['Result']
diff --git a/src/compas_fd/fd/result.py b/src/compas_fd/result/result.py
similarity index 100%
rename from src/compas_fd/fd/result.py
rename to src/compas_fd/result/result.py