Cleanup

Change naming, refactor unpacking of numdata

scripts/test_fd_constraint_solver.py

@@ -40,8 +40,8 @@
 
 # set up iterative constraint solver
 numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads)
-solver = FDConstraintSolver(numdata, constraints, max_iter=30,
-                            tol_res=1E-3, tol_dxyz=1E-3)
+solver = FDConstraintSolver(numdata, constraints, kmax=30,
+                            tol_res=1E-3, tol_disp=1E-3)
 
 # run solver
 result = solver()

scripts/test_fd_solver.py

@@ -30,6 +30,7 @@
 
 # run solver
 result = solver()
+print(solver.iter_count)
 
 # update mesh
 for index, vertex in enumerate(mesh.vertices()):

src/compas_fd/fd/fd_iter_numpy.py

@@ -22,14 +22,14 @@ def fd_iter_numpy(*,
                   forcedensities: List[float],
                   loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None,
                   constraints: Sequence[Constraint],
-                  max_iter: int = 100,
+                  kmax: int = 100,
                   tol_res: float = 1E-3,
-                  tol_dxyz: float = 1E-3
+                  tol_disp: float = 1E-3
                   ) -> Result:
     """Iteratively compute the equilibrium coordinates of a system of vertices connected by edges.
     Vertex constraints are recomputed at each iteration.
     """
     numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads)
-    solver = FDConstraintSolver(numdata, constraints, max_iter, tol_res, tol_dxyz)
+    solver = FDConstraintSolver(numdata, constraints, kmax, tol_res, tol_disp)
     result = solver()
     return result

src/compas_fd/fd/mesh_fd_iter_numpy.py

@@ -33,7 +33,7 @@ def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd
     constraints = list(mesh.vertices_attribute('constraint'))
 
     numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads)
-    solver = FDConstraintSolver(numdata, constraints, max_iter=100, tol_res=1E-3, tol_dxyz=1E-3)
+    solver = FDConstraintSolver(numdata, constraints, kmax=100, tol_res=1E-3, tol_disp=1E-3)
     result = solver()
 
     _update_mesh(mesh, result)

src/compas_fd/numdata/fd_numerical_data.py

@@ -8,7 +8,6 @@
 from nptyping import NDArray
 
 from dataclasses import dataclass
-from dataclasses import astuple
 
 from numpy import asarray
 from numpy import float64
@@ -40,9 +39,6 @@ class FDNumericalData(NumericalData):
     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]],

src/compas_fd/solvers/fd_constraint_solver.py

@@ -19,45 +19,43 @@ class FDConstraintSolver(Solver):
     def __init__(self,
                  numdata: FDNumericalData,
                  constraints: Sequence[Constraint] = None,
-                 max_iter: int = 100,
-                 tol_dxyz: float = 1E-3,
+                 kmax: int = 100,
                  tol_res: float = 1E-3,
+                 tol_disp: float = 1E-3,
                  **kwargs
                  ) -> None:
-        super(FDConstraintSolver, self).__init__(numdata, max_iter, **kwargs)
+        super(FDConstraintSolver, self).__init__(numdata, kmax, **kwargs)
         self.constraints = constraints
-        self.tol_dxyz = tol_dxyz
         self.tol_res = tol_res
+        self.tol_disp = tol_disp
 
     def solve(self) -> None:
         """Apply force density algorithm for a single iteration."""
-        free, fixed, xyz, C, q, Q, p, A, Ai, Af, *_ = self.numdata
-        b = p[free] - Af.dot(xyz[fixed])
-        xyz[free] = spsolve(Ai, b)
-        self.numdata.residuals = p - A.dot(xyz)
-        self.numdata.xyz_prev = xyz.copy()
-        self.numdata.xyz = xyz
+        nd = self.numdata
+        b = nd.p[nd.free] - nd.Af.dot(nd.xyz[nd.fixed])
+        nd.xyz[nd.free] = spsolve(nd.Ai, b)
+        nd.residuals = nd.p - nd.A.dot(nd.xyz)
+        nd.xyz_prev = nd.xyz.copy()
         self._update_constraints()
 
     def _update_constraints(self):
         """Update all vertex constraints by the residuals of the current iteration,
         and store their updated vertex coordinates in the numdata attribute.
         """
-        xyz = self.numdata.xyz
-        residuals = self.numdata.residuals
+        nd = self.numdata
         for vertex, constraint in enumerate(self.constraints):
             if not constraint:
                 continue
-            constraint.location = xyz[vertex]
-            constraint.residual = residuals[vertex]
-            xyz[vertex] = constraint.location + constraint.tangent * 0.5
-        self.numdata.tangent_residuals = asarray([c.tangent for c in self.constraints if c])
-        self.numdata.xyz = xyz
+            constraint.location = nd.xyz[vertex]
+            constraint.residual = nd.residuals[vertex]
+            nd.xyz[vertex] = constraint.location + constraint.tangent * 0.5
+        nd.tangent_residuals = asarray([c.tangent for c in self.constraints if c])
 
     @property
     def is_converged(self) -> bool:
         """Verify if all convergence criteria are met."""
-        return self._is_converged_residuals and self._is_converged_xyz
+        return (self._is_converged_residuals and
+                self._is_converged_displacements)
 
     @property
     def _is_converged_residuals(self) -> bool:
@@ -69,18 +67,16 @@ def _is_converged_residuals(self) -> bool:
         return max_res < self.tol_res
 
     @property
-    def _is_converged_xyz(self) -> bool:
+    def _is_converged_displacements(self) -> bool:
         """Verify whether the maximum coordinate displacement
         between consecutive iterations is within tolerance.
         """
-        new_xyz = self.numdata.xyz
-        old_xyz = self.numdata.xyz_prev
-        max_dxyz = max(norm(new_xyz - old_xyz, axis=1))
-        return max_dxyz < self.tol_dxyz
+        nd = self.numdata
+        max_dxyz = max(norm(nd.xyz - nd.xyz_prev, axis=1))
+        return max_dxyz < self.tol_disp
 
     def post_process(self) -> None:
         """Compute dependent variables after ending the solver loop."""
-        _, _, xyz, C, q, *_ = self.numdata
-        lengths = normrow(C.dot(xyz))
-        self.numdata.forces = q * lengths
-        self.numdata.lengths = lengths
+        nd = self.numdata
+        nd.lengths = normrow(nd.C.dot(nd.xyz))
+        nd.forces = nd.q * nd.lengths

src/compas_fd/solvers/fd_solver.py

@@ -16,10 +16,9 @@ def __init__(self, numdata: FDNumericalData, **kwargs) -> None:
 
     def solve(self) -> None:
         """Apply force density algorithm for a single iteration."""
-        free, fixed, xyz, C, q, Q, p, _, Ai, Af, *_ = self.numdata
-        b = p[free] - Af.dot(xyz[fixed])
-        xyz[free] = spsolve(Ai, b)
-        self.numdata.xyz = xyz
+        nd = self.numdata
+        b = nd.p[nd.free] - nd.Af.dot(nd.xyz[nd.fixed])
+        nd.xyz[nd.free] = spsolve(nd.Ai, b)
 
     @property
     def is_converged(self) -> bool:
@@ -28,8 +27,7 @@ def is_converged(self) -> bool:
 
     def post_process(self) -> None:
         """Compute dependent variables after ending solver."""
-        _, _, xyz, C, q, Q, p, A, *_ = self.numdata
-        lengths = normrow(C.dot(xyz))
-        self.numdata.lengths = lengths
-        self.numdata.forces = q * lengths
-        self.numdata.residuals = p - A.dot(xyz)
+        nd = self.numdata
+        nd.lengths = normrow(nd.C.dot(nd.xyz))
+        nd.forces = nd.q * nd.lengths
+        nd.residuals = nd.p - nd.A.dot(nd.xyz)

src/compas_fd/solvers/solver.py

@@ -9,17 +9,17 @@ class Solver:
 
     def __init__(self,
                  numdata: NumericalData,
-                 max_iter: int = 100,
+                 kmax: int = 100,
                  **kwargs
                  ) -> None:
         self.numdata = numdata
-        self.max_iter = max_iter
+        self.kmax = kmax
         self.iter_count = 0
         self.result = None
 
     def __call__(self) -> Result:
         """Iteratively apply the solver algorithm."""
-        for self.iter_count in range(1, self.max_iter + 1):
+        for self.iter_count in range(1, self.kmax + 1):
             self.solve()
             if self.is_converged:
                 break