diff --git a/scripts/test_fd_constraint_solver.py b/scripts/test_fd_constraint_solver.py index 3ee17c7a..692dbf8f 100644 --- a/scripts/test_fd_constraint_solver.py +++ b/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() diff --git a/scripts/test_fd_solver.py b/scripts/test_fd_solver.py index ad177c1f..e6e12860 100644 --- a/scripts/test_fd_solver.py +++ b/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()): diff --git a/src/compas_fd/fd/fd_iter_numpy.py b/src/compas_fd/fd/fd_iter_numpy.py index 83b54eb7..07363bf8 100644 --- a/src/compas_fd/fd/fd_iter_numpy.py +++ b/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 diff --git a/src/compas_fd/fd/mesh_fd_iter_numpy.py b/src/compas_fd/fd/mesh_fd_iter_numpy.py index 5ba4ac05..94113c2e 100644 --- a/src/compas_fd/fd/mesh_fd_iter_numpy.py +++ b/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) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index e3fd4689..be33c9ac 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/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]], diff --git a/src/compas_fd/solvers/fd_constraint_solver.py b/src/compas_fd/solvers/fd_constraint_solver.py index da07d638..67dfbad9 100644 --- a/src/compas_fd/solvers/fd_constraint_solver.py +++ b/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 diff --git a/src/compas_fd/solvers/fd_solver.py b/src/compas_fd/solvers/fd_solver.py index 0a87e2e9..b72aa4b7 100644 --- a/src/compas_fd/solvers/fd_solver.py +++ b/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) diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py index df5aee16..ced00b2e 100644 --- a/src/compas_fd/solvers/solver.py +++ b/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