From 1219c4c336ee80efc5eb65303a0fac1785f76851 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 11:21:54 +0100 Subject: [PATCH 01/26] Add numerical data class Class for storing numerical arrays between iterations. Add iterative fd solver Constraint update at each iteration. Add test for iterative solver Lint remove whitespace after unpacking operator Lint clean whitespace update python req to >= 3.7 Update build versions remove python version 3.6 Add mesh iterative solver Clean mesh iterative solver Add failsafe for residuals fix residuals checking Clean Refactor for general readability. --- scripts/test_iterative_solver.py | 77 ++++++++++++++++++ src/compas_fd/fd/fd_iter_numpy.py | 104 +++++++++++++++++++++++++ src/compas_fd/fd/fd_numerical_data.py | 4 +- src/compas_fd/fd/mesh_fd_iter_numpy.py | 60 ++++++++++++++ 4 files changed, 242 insertions(+), 3 deletions(-) create mode 100644 scripts/test_iterative_solver.py create mode 100644 src/compas_fd/fd/fd_iter_numpy.py create mode 100644 src/compas_fd/fd/mesh_fd_iter_numpy.py diff --git a/scripts/test_iterative_solver.py b/scripts/test_iterative_solver.py new file mode 100644 index 00000000..9638d5b6 --- /dev/null +++ b/scripts/test_iterative_solver.py @@ -0,0 +1,77 @@ +import compas +from compas.geometry import Vector, Plane, Point, Line +from compas_fd.datastructures import CableMesh +from compas_fd.fd import fd_iter_numpy + +from compas_fd.constraints import Constraint + +from compas_view2.app import App +from compas_view2.objects import Object, MeshObject + +Object.register(CableMesh, MeshObject) + +mesh = CableMesh.from_obj(compas.get('faces.obj')) + +mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) +mesh.vertices_attribute('t', 0.0) + +# line constraint +vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] +line = Line(Point(10, 0, 0), Point(5, 10, 5)) +constraint = Constraint(line) +mesh.vertex_attribute(vertex, 'constraint', constraint) + +# plane constraint +vertex = list(mesh.vertices_where({'x': 0, 'y': 0}))[0] +plane = Plane((-1, 0, 0), (2, 1, 1)) +constraint = Constraint(plane) +mesh.vertex_attribute(vertex, 'constraint', constraint) + +# input solver parameters +vertex_index = mesh.vertex_index() +vertices = mesh.vertices_attributes('xyz') +fixed = list(mesh.vertices_where({'is_anchor': True})) +edges = [(vertex_index[u], vertex_index[v]) for + u, v in mesh.edges_where({'_is_edge': True})] +forcedensities = mesh.edges_attribute('q') +loads = mesh.vertices_attributes(['px', 'py', 'pz']) +constraints = list(mesh.vertices_attribute('constraint')) + +# solver +result = fd_iter_numpy(vertices=vertices, + fixed=fixed, + edges=edges, + forcedensities=forcedensities, + loads=loads, + constraints=constraints, + max_iter=100, + tol_res=1E-3, + tol_dxyz=1E-1) + +# update mesh +for index, vertex in enumerate(mesh.vertices()): + mesh.vertex_attributes(vertex, 'xyz', result.vertices[index]) + mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], result.residuals[index]) + + +# ============================================================================== +# Viz +# ============================================================================== + +viewer = App() + +viewer.add(mesh) + +for vertex in fixed: + viewer.add(Point(* mesh.vertex_attributes(vertex, 'xyz')), size=20, color=(0, 0, 0)) + +for constraint in constraints: + if constraint: + viewer.add(constraint.geometry, linewidth=5, linecolor=(0, 1, 1)) + +for vertex in fixed: + a = Point(* mesh.vertex_attributes(vertex, 'xyz')) + b = a - Vector(* mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'])) + viewer.add(Line(a, b), linewidth=3, linecolor=(0, 0, 1)) + +viewer.run() diff --git a/src/compas_fd/fd/fd_iter_numpy.py b/src/compas_fd/fd/fd_iter_numpy.py new file mode 100644 index 00000000..050861e4 --- /dev/null +++ b/src/compas_fd/fd/fd_iter_numpy.py @@ -0,0 +1,104 @@ +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 scipy.linalg import norm +from scipy.sparse.linalg import spsolve + +from compas.numerical import normrow + +from .fd_numerical_data import FDNumericalData +from .result import Result +from ..constraints import Constraint + + +def fd_iter_numpy(*, + 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, + constraints: Sequence[Constraint], + max_iter: int = 100, + tol_res: float = 1E-3, + tol_dxyz: 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) + + for k in range(max_iter): + xyz_prev = numdata.xyz + _solve_fd(numdata) + _update_constraints(numdata, constraints) + if (_is_converged_residuals(numdata.tangent_residuals, tol_res) and + _is_converged_xyz(xyz_prev, numdata.xyz, tol_dxyz)): + break + + _post_process_fd(numdata) + return numdata.to_result() + + +def _solve_fd(numdata: FDNumericalData) -> None: + """Solve a single iteration for the equilibrium coordinates of a system. + All updated numerical arrays are stored in the numdata parameter. + """ + free, fixed, xyz, C, _, Q, p, Ai, Af, *_ = numdata + b = p[free] - Af.dot(xyz[fixed]) + xyz[free] = spsolve(Ai, b) + numdata.residuals = p - C.T.dot(Q).dot(C).dot(xyz) + numdata.xyz = xyz + + +def _post_process_fd(numdata: FDNumericalData) -> None: + """Compute dependent numerical arrays from the numerical data after running solver. + All updated numerical arrays are stored in the numdata parameter. + """ + _, _, xyz, C, q, *_ = numdata + lengths = normrow(C.dot(xyz)) + numdata.forces = q * lengths + numdata.lengths = lengths + + +def _update_constraints(numdata: FDNumericalData, + constraints: Sequence[Constraint]) -> None: + """Update all vertex constraints by the residuals of the current iteration, + and store their updated vertex coordinates in the numdata parameter. + """ + xyz = numdata.xyz + residuals = numdata.residuals + for vertex, constraint in enumerate(constraints): + if not constraint: + continue + constraint.location = xyz[vertex] + constraint.residual = residuals[vertex] + xyz[vertex] = constraint.location + constraint.tangent * 0.5 + numdata.tangent_residuals = asarray([c.tangent for c in constraints if c]) + numdata.xyz = xyz + + +def _is_converged_residuals(residuals: NDArray[(Any, 3), float64], + tol_res: float) -> bool: + """Verify whether the maximum constraint residual is within tolerance.""" + if residuals is None or not residuals.any(): + return True + max_res = max(norm(residuals, axis=1)) + return max_res < tol_res + + +def _is_converged_xyz(old_xyz: NDArray[(Any, 3), float64], + new_xyz: NDArray[(Any, 3), float64], + tol_xyz: float) -> bool: + """Verify whether the maximum coordinate displacement + between consecutive iterations is within tolerance. + """ + max_dxyz = max(norm(new_xyz - old_xyz, axis=1)) + return max_dxyz < tol_xyz diff --git a/src/compas_fd/fd/fd_numerical_data.py b/src/compas_fd/fd/fd_numerical_data.py index fb061df5..a322d80c 100644 --- a/src/compas_fd/fd/fd_numerical_data.py +++ b/src/compas_fd/fd/fd_numerical_data.py @@ -30,7 +30,6 @@ class FDNumericalData: 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 @@ -59,10 +58,9 @@ def from_params(cls, 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) + return cls(free, fixed, xyz, C, q, Q, p, Ai, Af) @classmethod def from_mesh(cls, mesh): diff --git a/src/compas_fd/fd/mesh_fd_iter_numpy.py b/src/compas_fd/fd/mesh_fd_iter_numpy.py new file mode 100644 index 00000000..d97b7339 --- /dev/null +++ b/src/compas_fd/fd/mesh_fd_iter_numpy.py @@ -0,0 +1,60 @@ +from numpy import array, asarray +from numpy import float64 + +import compas_fd +from .fd_iter_numpy import fd_iter_numpy + + +def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': + """Iteratively find the equilibrium shape of a mesh for the given force densities. + + Parameters + ---------- + mesh : :class:`compas_fd.datastructures.CableMesh` + The mesh to equilibriate. + + Returns + ------- + :class:`compas_fd.datastructures.CableMesh` + The function updates the mesh in place, + but returns a reference to the updated mesh as well + for compatibility with RPCs. + + """ + k_i = mesh.key_index() + vertices = array(mesh.vertices_attributes('xyz'), dtype=float64) + fixed = [k_i[v] for v in mesh.vertices_where({'is_anchor': True})] + edges = [(k_i[u], k_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] + forcedensities = asarray([attr['q'] for key, attr in mesh.edges_where({'_is_edge': True}, True)], + dtype=float64).reshape((-1, 1)) + loads = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) + constraints = list(mesh.vertices_attribute('constraint')) + + result = fd_iter_numpy(vertices=vertices, + fixed=fixed, + edges=edges, + forcedensities=forcedensities, + loads=loads, + constraints=constraints, + max_iter=100, + tol_res=1E-3, + tol_xyz=1E-3) + + _update_mesh(mesh, result) + + return mesh + + +def _update_mesh(mesh, result): + for key, attr in mesh.vertices(True): + index = mesh.key_index()[key] + attr['x'] = result.vertices[index, 0] + attr['y'] = result.vertices[index, 1] + attr['z'] = result.vertices[index, 2] + attr['_rx'] = result.residuals[index, 0] + attr['_ry'] = result.residuals[index, 1] + attr['_rz'] = result.residuals[index, 2] + + for index, (key, attr) in enumerate(mesh.edges_where({'_is_edge': True}, True)): + attr['_f'] = result.forces[index, 0] + attr['_l'] = result.lenghts[index, 0] From 0cffdf2495a37571568cd30435f110297945041e Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:16:49 +0100 Subject: [PATCH 02/26] Add abstract Solver class Base solver manager class for running iterative algorithms, storing numerical data, and results --- src/compas_fd/solvers/solver.py | 41 +++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 src/compas_fd/solvers/solver.py diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py new file mode 100644 index 00000000..df5aee16 --- /dev/null +++ b/src/compas_fd/solvers/solver.py @@ -0,0 +1,41 @@ +from ..numdata import NumericalData +from ..result import Result + + +class Solver: + """Manager class for running iterative algorithms, + updating numerical data and storing results. + """ + + def __init__(self, + numdata: NumericalData, + max_iter: int = 100, + **kwargs + ) -> None: + self.numdata = numdata + self.max_iter = max_iter + 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): + self.solve() + if self.is_converged: + break + self.post_process() + self.result = self.numdata.to_result() + return self.result + + def solve(self) -> None: + """Apply solver algorithm for a single iteration.""" + raise NotImplementedError + + @property + def is_converged(self) -> bool: + """Verify if all convergence criteria are met.""" + raise NotImplementedError + + def post_process(self) -> None: + """Compute dependent variables after ending the solver loop.""" + raise NotImplementedError From 12f94683c9f3ac15c7a6a87658a23903c6181316 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:17:13 +0100 Subject: [PATCH 03/26] Add single iteration fd solver --- src/compas_fd/solvers/fd_solver.py | 35 ++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 src/compas_fd/solvers/fd_solver.py diff --git a/src/compas_fd/solvers/fd_solver.py b/src/compas_fd/solvers/fd_solver.py new file mode 100644 index 00000000..0a87e2e9 --- /dev/null +++ b/src/compas_fd/solvers/fd_solver.py @@ -0,0 +1,35 @@ +from scipy.sparse.linalg import spsolve + +from compas.numerical import normrow + +from .solver import Solver +from ..numdata import FDNumericalData + + +class FDSolver(Solver): + """Manager class for running a single step force density algorithm, + updating numerical data and storing results. + """ + + def __init__(self, numdata: FDNumericalData, **kwargs) -> None: + super(FDSolver, self).__init__(numdata, 1, **kwargs) + + 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 + + @property + def is_converged(self) -> bool: + """Verify if all convergence criteria are met.""" + return True + + 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) From a08166b5655c638ba0abc95c96dc666cebb3cc12 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:17:28 +0100 Subject: [PATCH 04/26] Add iterative constraint fd solver --- src/compas_fd/solvers/fd_constraint_solver.py | 86 +++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 src/compas_fd/solvers/fd_constraint_solver.py diff --git a/src/compas_fd/solvers/fd_constraint_solver.py b/src/compas_fd/solvers/fd_constraint_solver.py new file mode 100644 index 00000000..da07d638 --- /dev/null +++ b/src/compas_fd/solvers/fd_constraint_solver.py @@ -0,0 +1,86 @@ +from typing import Sequence + +from numpy import asarray +from numpy.linalg import norm +from scipy.sparse.linalg import spsolve + +from compas.numerical import normrow + +from .solver import Solver +from ..numdata import FDNumericalData +from ..constraints import Constraint + + +class FDConstraintSolver(Solver): + """Manager class for running iterative constraint force density algorithm, + updating numerical data and storing results. + """ + + def __init__(self, + numdata: FDNumericalData, + constraints: Sequence[Constraint] = None, + max_iter: int = 100, + tol_dxyz: float = 1E-3, + tol_res: float = 1E-3, + **kwargs + ) -> None: + super(FDConstraintSolver, self).__init__(numdata, max_iter, **kwargs) + self.constraints = constraints + self.tol_dxyz = tol_dxyz + self.tol_res = tol_res + + 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 + 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 + 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 + + @property + def is_converged(self) -> bool: + """Verify if all convergence criteria are met.""" + return self._is_converged_residuals and self._is_converged_xyz + + @property + def _is_converged_residuals(self) -> bool: + """Verify whether the maximum constraint residual is within tolerance.""" + residuals = self.numdata.tangent_residuals + if residuals is None or not residuals.any(): + return True + max_res = max(norm(residuals, axis=1)) + return max_res < self.tol_res + + @property + def _is_converged_xyz(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 + + 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 From 277822aa680f7defeb39b6c7f1cb1947bc0359de Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:17:46 +0100 Subject: [PATCH 05/26] Add solvers init --- src/compas_fd/solvers/__init__.py | 37 +++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/compas_fd/solvers/__init__.py diff --git a/src/compas_fd/solvers/__init__.py b/src/compas_fd/solvers/__init__.py new file mode 100644 index 00000000..c0f12e8b --- /dev/null +++ b/src/compas_fd/solvers/__init__.py @@ -0,0 +1,37 @@ +""" +****************** +compas_fd.solvers +****************** + +.. currentmodule:: compas_fd.solvers + + +Classes +======= + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + Solver + FDSolver + FDConstraintSolver + +""" +from __future__ import print_function +from __future__ import absolute_import +from __future__ import division + +import compas +from .solver import Solver +from .fd_solver import FDSolver +from .fd_constraint_solver import FDConstraintSolver + +if not compas.IPY: + pass + +__all__ = ['Solver'] + +if not compas.IPY: + __all__ += ['FDSolver', + 'FDConstraintSolver'] From bb3ca26ca44f354618f72b2a5cf69b4a1ac1ee7a Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:18:05 +0100 Subject: [PATCH 06/26] Relocate Result class --- src/compas_fd/result/__init__.py | 25 +++++++++++++++++++++++++ src/compas_fd/{fd => result}/result.py | 0 2 files changed, 25 insertions(+) create mode 100644 src/compas_fd/result/__init__.py rename src/compas_fd/{fd => result}/result.py (100%) 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 From e44b8363c55dc3faa2bd6cb2f501d5ef0cd573b2 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:19:04 +0100 Subject: [PATCH 07/26] Refactor NumericalData class Relocate and refactor with base class and FD specific class --- src/compas_fd/numdata/__init__.py | 33 +++++++++++++++++++ .../{fd => numdata}/fd_numerical_data.py | 12 ++++--- src/compas_fd/numdata/numerical_data.py | 26 +++++++++++++++ 3 files changed, 67 insertions(+), 4 deletions(-) create mode 100644 src/compas_fd/numdata/__init__.py rename src/compas_fd/{fd => numdata}/fd_numerical_data.py (88%) create mode 100644 src/compas_fd/numdata/numerical_data.py 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/fd/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py similarity index 88% rename from src/compas_fd/fd/fd_numerical_data.py rename to src/compas_fd/numdata/fd_numerical_data.py index a322d80c..e3fd4689 100644 --- a/src/compas_fd/fd/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -17,11 +17,12 @@ from compas.numerical import connectivity_matrix -from .result import Result +from .numerical_data import NumericalData +from ..result import Result @dataclass -class FDNumericalData: +class FDNumericalData(NumericalData): """Stores numerical data used by the force density algorithms.""" free: int fixed: int @@ -30,6 +31,7 @@ class FDNumericalData: 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 @@ -47,7 +49,8 @@ def from_params(cls, fixed: List[int], edges: List[Tuple[int, int]], forcedensities: List[float], - loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None): + 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)) @@ -58,9 +61,10 @@ def from_params(cls, 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, Ai, Af) + return cls(free, fixed, xyz, C, q, Q, p, A, Ai, Af) @classmethod def from_mesh(cls, mesh): diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py new file mode 100644 index 00000000..752cac39 --- /dev/null +++ b/src/compas_fd/numdata/numerical_data.py @@ -0,0 +1,26 @@ +from dataclasses import dataclass +from dataclasses import astuple + +from ..result import Result + + +@dataclass +class NumericalData: + """Storage class for numerical arrays used by solver algorithms.""" + + def __iter__(self): + return iter(astuple(self)) + + @classmethod + def from_params(cls, *args, **kwargs): + """Construct nuerical arrays from algorithm input parameters.""" + raise NotImplementedError + + @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 From 2a670e5ddc9b55d91ef54b117e4c8c69d17fca56 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:19:33 +0100 Subject: [PATCH 08/26] Refactor fd_iter_numpy Call new solver backend implementation --- src/compas_fd/fd/fd_iter_numpy.py | 81 +++---------------------------- 1 file changed, 6 insertions(+), 75 deletions(-) diff --git a/src/compas_fd/fd/fd_iter_numpy.py b/src/compas_fd/fd/fd_iter_numpy.py index 050861e4..83b54eb7 100644 --- a/src/compas_fd/fd/fd_iter_numpy.py +++ b/src/compas_fd/fd/fd_iter_numpy.py @@ -7,16 +7,12 @@ from typing_extensions import Annotated from nptyping import NDArray -from numpy import asarray from numpy import float64 -from scipy.linalg import norm -from scipy.sparse.linalg import spsolve -from compas.numerical import normrow - -from .fd_numerical_data import FDNumericalData -from .result import Result +from ..result import Result from ..constraints import Constraint +from ..numdata import FDNumericalData +from ..solvers import FDConstraintSolver def fd_iter_numpy(*, @@ -34,71 +30,6 @@ def fd_iter_numpy(*, Vertex constraints are recomputed at each iteration. """ numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) - - for k in range(max_iter): - xyz_prev = numdata.xyz - _solve_fd(numdata) - _update_constraints(numdata, constraints) - if (_is_converged_residuals(numdata.tangent_residuals, tol_res) and - _is_converged_xyz(xyz_prev, numdata.xyz, tol_dxyz)): - break - - _post_process_fd(numdata) - return numdata.to_result() - - -def _solve_fd(numdata: FDNumericalData) -> None: - """Solve a single iteration for the equilibrium coordinates of a system. - All updated numerical arrays are stored in the numdata parameter. - """ - free, fixed, xyz, C, _, Q, p, Ai, Af, *_ = numdata - b = p[free] - Af.dot(xyz[fixed]) - xyz[free] = spsolve(Ai, b) - numdata.residuals = p - C.T.dot(Q).dot(C).dot(xyz) - numdata.xyz = xyz - - -def _post_process_fd(numdata: FDNumericalData) -> None: - """Compute dependent numerical arrays from the numerical data after running solver. - All updated numerical arrays are stored in the numdata parameter. - """ - _, _, xyz, C, q, *_ = numdata - lengths = normrow(C.dot(xyz)) - numdata.forces = q * lengths - numdata.lengths = lengths - - -def _update_constraints(numdata: FDNumericalData, - constraints: Sequence[Constraint]) -> None: - """Update all vertex constraints by the residuals of the current iteration, - and store their updated vertex coordinates in the numdata parameter. - """ - xyz = numdata.xyz - residuals = numdata.residuals - for vertex, constraint in enumerate(constraints): - if not constraint: - continue - constraint.location = xyz[vertex] - constraint.residual = residuals[vertex] - xyz[vertex] = constraint.location + constraint.tangent * 0.5 - numdata.tangent_residuals = asarray([c.tangent for c in constraints if c]) - numdata.xyz = xyz - - -def _is_converged_residuals(residuals: NDArray[(Any, 3), float64], - tol_res: float) -> bool: - """Verify whether the maximum constraint residual is within tolerance.""" - if residuals is None or not residuals.any(): - return True - max_res = max(norm(residuals, axis=1)) - return max_res < tol_res - - -def _is_converged_xyz(old_xyz: NDArray[(Any, 3), float64], - new_xyz: NDArray[(Any, 3), float64], - tol_xyz: float) -> bool: - """Verify whether the maximum coordinate displacement - between consecutive iterations is within tolerance. - """ - max_dxyz = max(norm(new_xyz - old_xyz, axis=1)) - return max_dxyz < tol_xyz + solver = FDConstraintSolver(numdata, constraints, max_iter, tol_res, tol_dxyz) + result = solver() + return result From 0cef4c944660deaf10038fc30ba27bfb4ea3c0e8 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:20:05 +0100 Subject: [PATCH 09/26] Update result location reference --- src/compas_fd/fd/fd_numpy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/compas_fd/fd/fd_numpy.py b/src/compas_fd/fd/fd_numpy.py index 241d917d..310d15b9 100644 --- a/src/compas_fd/fd/fd_numpy.py +++ b/src/compas_fd/fd/fd_numpy.py @@ -15,7 +15,7 @@ from compas.numerical import connectivity_matrix from compas.numerical import normrow -from .result import Result +from ..result import Result def fd_numpy(*, From 02387b07bd2dc28e75dfd6d3a1fbf14bbf8e7c5b Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:20:48 +0100 Subject: [PATCH 10/26] Update mesh_fd solvers Reference updated solver backend --- src/compas_fd/__init__.py | 3 +++ src/compas_fd/fd/mesh_fd_iter_numpy.py | 17 ++++++----------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/compas_fd/__init__.py b/src/compas_fd/__init__.py index ee078b11..debe0ce6 100644 --- a/src/compas_fd/__init__.py +++ b/src/compas_fd/__init__.py @@ -12,6 +12,9 @@ compas_fd.datastructures compas_fd.fd compas_fd.loads + compas_fd.numdata + compas_fd.results + compas_fd.solvers """ diff --git a/src/compas_fd/fd/mesh_fd_iter_numpy.py b/src/compas_fd/fd/mesh_fd_iter_numpy.py index d97b7339..5ba4ac05 100644 --- a/src/compas_fd/fd/mesh_fd_iter_numpy.py +++ b/src/compas_fd/fd/mesh_fd_iter_numpy.py @@ -2,7 +2,9 @@ from numpy import float64 import compas_fd -from .fd_iter_numpy import fd_iter_numpy + +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDConstraintSolver def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': @@ -30,18 +32,11 @@ def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd loads = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) constraints = list(mesh.vertices_attribute('constraint')) - result = fd_iter_numpy(vertices=vertices, - fixed=fixed, - edges=edges, - forcedensities=forcedensities, - loads=loads, - constraints=constraints, - max_iter=100, - tol_res=1E-3, - tol_xyz=1E-3) + numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) + solver = FDConstraintSolver(numdata, constraints, max_iter=100, tol_res=1E-3, tol_dxyz=1E-3) + result = solver() _update_mesh(mesh, result) - return mesh From 254f8d5094e2dd95e6a52bc65658b6c40c305bb4 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 15:21:36 +0100 Subject: [PATCH 11/26] Update test scripts --- scripts/test_fd_constraint_solver.py | 21 +++++---- ..._iterative_solver.py => test_fd_solver.py} | 43 +++++-------------- 2 files changed, 21 insertions(+), 43 deletions(-) rename scripts/{test_iterative_solver.py => test_fd_solver.py} (56%) diff --git a/scripts/test_fd_constraint_solver.py b/scripts/test_fd_constraint_solver.py index 34a4e80d..3ee17c7a 100644 --- a/scripts/test_fd_constraint_solver.py +++ b/scripts/test_fd_constraint_solver.py @@ -1,9 +1,10 @@ import compas from compas.geometry import Vector, Plane, Point, Line from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_constrained_numpy from compas_fd.constraints import Constraint +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDConstraintSolver from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -37,16 +38,14 @@ loads = mesh.vertices_attributes(['px', 'py', 'pz']) constraints = list(mesh.vertices_attribute('constraint')) -# solver -result = fd_constrained_numpy(vertices=vertices, - fixed=fixed, - edges=edges, - forcedensities=forcedensities, - loads=loads, - constraints=constraints, - kmax=100, - tol_res=1E-3, - tol_disp=1E-3) +# 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) + +# run solver +result = solver() +print(solver.iter_count) # update mesh for index, vertex in enumerate(mesh.vertices()): diff --git a/scripts/test_iterative_solver.py b/scripts/test_fd_solver.py similarity index 56% rename from scripts/test_iterative_solver.py rename to scripts/test_fd_solver.py index 9638d5b6..ad177c1f 100644 --- a/scripts/test_iterative_solver.py +++ b/scripts/test_fd_solver.py @@ -1,9 +1,9 @@ import compas -from compas.geometry import Vector, Plane, Point, Line +from compas.geometry import Point, Line, Vector from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_iter_numpy -from compas_fd.constraints import Constraint +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDSolver from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -15,19 +15,7 @@ mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) mesh.vertices_attribute('t', 0.0) -# line constraint -vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] -line = Line(Point(10, 0, 0), Point(5, 10, 5)) -constraint = Constraint(line) -mesh.vertex_attribute(vertex, 'constraint', constraint) - -# plane constraint -vertex = list(mesh.vertices_where({'x': 0, 'y': 0}))[0] -plane = Plane((-1, 0, 0), (2, 1, 1)) -constraint = Constraint(plane) -mesh.vertex_attribute(vertex, 'constraint', constraint) - -# input solver parameters +# input parameters vertex_index = mesh.vertex_index() vertices = mesh.vertices_attributes('xyz') fixed = list(mesh.vertices_where({'is_anchor': True})) @@ -35,18 +23,13 @@ u, v in mesh.edges_where({'_is_edge': True})] forcedensities = mesh.edges_attribute('q') loads = mesh.vertices_attributes(['px', 'py', 'pz']) -constraints = list(mesh.vertices_attribute('constraint')) - -# solver -result = fd_iter_numpy(vertices=vertices, - fixed=fixed, - edges=edges, - forcedensities=forcedensities, - loads=loads, - constraints=constraints, - max_iter=100, - tol_res=1E-3, - tol_dxyz=1E-1) + +# set up single iteration solver +numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) +solver = FDSolver(numdata) + +# run solver +result = solver() # update mesh for index, vertex in enumerate(mesh.vertices()): @@ -65,10 +48,6 @@ for vertex in fixed: viewer.add(Point(* mesh.vertex_attributes(vertex, 'xyz')), size=20, color=(0, 0, 0)) -for constraint in constraints: - if constraint: - viewer.add(constraint.geometry, linewidth=5, linecolor=(0, 1, 1)) - for vertex in fixed: a = Point(* mesh.vertex_attributes(vertex, 'xyz')) b = a - Vector(* mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'])) From cada522dfd04be5f0988876537615096e73dc2c2 Mon Sep 17 00:00:00 2001 From: Sam Date: Thu, 11 Nov 2021 17:40:14 +0100 Subject: [PATCH 12/26] Cleanup Change naming, refactor unpacking of numdata --- scripts/test_fd_constraint_solver.py | 4 +- scripts/test_fd_solver.py | 1 + src/compas_fd/fd/fd_iter_numpy.py | 6 +-- src/compas_fd/fd/mesh_fd_iter_numpy.py | 2 +- src/compas_fd/numdata/fd_numerical_data.py | 4 -- src/compas_fd/solvers/fd_constraint_solver.py | 50 +++++++++---------- src/compas_fd/solvers/fd_solver.py | 16 +++--- src/compas_fd/solvers/solver.py | 6 +-- 8 files changed, 40 insertions(+), 49 deletions(-) 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 From 3d7aafe31b31169aaecb286d6acc154efb06b0b6 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 17:54:00 +0100 Subject: [PATCH 13/26] Import cleanup --- src/compas_fd/fd/fd_constrained_numpy.py | 4 ++-- src/compas_fd/fd/fd_iter_numpy.py | 8 ++++---- src/compas_fd/fd/fd_numpy.py | 2 +- src/compas_fd/numdata/fd_numerical_data.py | 2 +- src/compas_fd/numdata/numerical_data.py | 2 +- src/compas_fd/solvers/fd_constraint_solver.py | 4 ++-- src/compas_fd/solvers/solver.py | 4 ++-- 7 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/compas_fd/fd/fd_constrained_numpy.py b/src/compas_fd/fd/fd_constrained_numpy.py index d2c3c032..b4272112 100644 --- a/src/compas_fd/fd/fd_constrained_numpy.py +++ b/src/compas_fd/fd/fd_constrained_numpy.py @@ -16,8 +16,8 @@ from compas_fd.constraints import Constraint -from .fd_numerical_data import FDNumericalData -from .result import Result +from compas_fd.numdata import FDNumericalData +from compas_fd.result import Result def fd_constrained_numpy(*, diff --git a/src/compas_fd/fd/fd_iter_numpy.py b/src/compas_fd/fd/fd_iter_numpy.py index 07363bf8..7dda9752 100644 --- a/src/compas_fd/fd/fd_iter_numpy.py +++ b/src/compas_fd/fd/fd_iter_numpy.py @@ -9,10 +9,10 @@ from numpy import float64 -from ..result import Result -from ..constraints import Constraint -from ..numdata import FDNumericalData -from ..solvers import FDConstraintSolver +from compas_fd.result import Result +from compas_fd.constraints import Constraint +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDConstraintSolver def fd_iter_numpy(*, diff --git a/src/compas_fd/fd/fd_numpy.py b/src/compas_fd/fd/fd_numpy.py index 310d15b9..a9da8ba3 100644 --- a/src/compas_fd/fd/fd_numpy.py +++ b/src/compas_fd/fd/fd_numpy.py @@ -15,7 +15,7 @@ from compas.numerical import connectivity_matrix from compas.numerical import normrow -from ..result import Result +from compas_fd.result import Result def fd_numpy(*, diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index be33c9ac..21020d28 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -16,8 +16,8 @@ from compas.numerical import connectivity_matrix +from compas_fd.result import Result from .numerical_data import NumericalData -from ..result import Result @dataclass diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py index 752cac39..573a935c 100644 --- a/src/compas_fd/numdata/numerical_data.py +++ b/src/compas_fd/numdata/numerical_data.py @@ -1,7 +1,7 @@ from dataclasses import dataclass from dataclasses import astuple -from ..result import Result +from compas_fd.result import Result @dataclass diff --git a/src/compas_fd/solvers/fd_constraint_solver.py b/src/compas_fd/solvers/fd_constraint_solver.py index 67dfbad9..09ade246 100644 --- a/src/compas_fd/solvers/fd_constraint_solver.py +++ b/src/compas_fd/solvers/fd_constraint_solver.py @@ -6,9 +6,9 @@ from compas.numerical import normrow +from compas_fd.numdata import FDNumericalData +from compas_fd.constraints import Constraint from .solver import Solver -from ..numdata import FDNumericalData -from ..constraints import Constraint class FDConstraintSolver(Solver): diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py index ced00b2e..3bd82bae 100644 --- a/src/compas_fd/solvers/solver.py +++ b/src/compas_fd/solvers/solver.py @@ -1,5 +1,5 @@ -from ..numdata import NumericalData -from ..result import Result +from compas_fd.numdata import NumericalData +from compas_fd.result import Result class Solver: From 8b8b3a181ef5f63d52d23bfd2ffd322a24c9d2fa Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 17:55:24 +0100 Subject: [PATCH 14/26] Simplify base solver --- src/compas_fd/solvers/solver.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py index 3bd82bae..3fb23e33 100644 --- a/src/compas_fd/solvers/solver.py +++ b/src/compas_fd/solvers/solver.py @@ -14,12 +14,12 @@ def __init__(self, ) -> None: self.numdata = numdata self.kmax = kmax - self.iter_count = 0 + self.kcount = 0 self.result = None def __call__(self) -> Result: - """Iteratively apply the solver algorithm.""" - for self.iter_count in range(1, self.kmax + 1): + """Iteratively apply the solve algorithm.""" + for self.kcount in range(1, self.kmax + 1): self.solve() if self.is_converged: break @@ -37,5 +37,5 @@ def is_converged(self) -> bool: raise NotImplementedError def post_process(self) -> None: - """Compute dependent variables after ending the solver loop.""" - raise NotImplementedError + """Callable after ending the solver loop.""" + pass From 740faec6f507ccdfd7673937d82e02288caa74a3 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 17:56:43 +0100 Subject: [PATCH 15/26] Simplify FDSolvers Change responsibility to post-process from solver to lazy calculation in FDNumericalData --- src/compas_fd/solvers/fd_constraint_solver.py | 21 +++++++------------ src/compas_fd/solvers/fd_solver.py | 15 +++---------- 2 files changed, 10 insertions(+), 26 deletions(-) diff --git a/src/compas_fd/solvers/fd_constraint_solver.py b/src/compas_fd/solvers/fd_constraint_solver.py index 09ade246..7128c976 100644 --- a/src/compas_fd/solvers/fd_constraint_solver.py +++ b/src/compas_fd/solvers/fd_constraint_solver.py @@ -4,8 +4,6 @@ from numpy.linalg import norm from scipy.sparse.linalg import spsolve -from compas.numerical import normrow - from compas_fd.numdata import FDNumericalData from compas_fd.constraints import Constraint from .solver import Solver @@ -32,24 +30,25 @@ def __init__(self, def solve(self) -> None: """Apply force density algorithm for a single iteration.""" 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) + b = nd.p[nd.free] - nd.Df.dot(nd.xyz[nd.fixed]) + nd.xyz[nd.free] = spsolve(nd.Di, b) nd.xyz_prev = nd.xyz.copy() self._update_constraints() - def _update_constraints(self): + def _update_constraints(self) -> None: """Update all vertex constraints by the residuals of the current iteration, and store their updated vertex coordinates in the numdata attribute. """ nd = self.numdata + xyz = nd.xyz for vertex, constraint in enumerate(self.constraints): if not constraint: continue - constraint.location = nd.xyz[vertex] + constraint.location = xyz[vertex] constraint.residual = nd.residuals[vertex] - nd.xyz[vertex] = constraint.location + constraint.tangent * 0.5 + xyz[vertex] = constraint.location + constraint.tangent * 0.5 nd.tangent_residuals = asarray([c.tangent for c in self.constraints if c]) + nd.xyz = xyz @property def is_converged(self) -> bool: @@ -74,9 +73,3 @@ def _is_converged_displacements(self) -> bool: 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.""" - 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 b72aa4b7..7371e89a 100644 --- a/src/compas_fd/solvers/fd_solver.py +++ b/src/compas_fd/solvers/fd_solver.py @@ -1,9 +1,7 @@ from scipy.sparse.linalg import spsolve -from compas.numerical import normrow - +from compas_fd.numdata import FDNumericalData from .solver import Solver -from ..numdata import FDNumericalData class FDSolver(Solver): @@ -17,17 +15,10 @@ def __init__(self, numdata: FDNumericalData, **kwargs) -> None: def solve(self) -> None: """Apply force density algorithm for a single iteration.""" nd = self.numdata - b = nd.p[nd.free] - nd.Af.dot(nd.xyz[nd.fixed]) - nd.xyz[nd.free] = spsolve(nd.Ai, b) + b = nd.p[nd.free] - nd.Df.dot(nd.xyz[nd.fixed]) + nd.xyz[nd.free] = spsolve(nd.Di, b) @property def is_converged(self) -> bool: """Verify if all convergence criteria are met.""" return True - - def post_process(self) -> None: - """Compute dependent variables after ending solver.""" - 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) From 55e575ca680b6679083cfd469916a4bfa72266b3 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 17:57:20 +0100 Subject: [PATCH 16/26] Simplify base NumericalData Change from dataclass to regular class --- src/compas_fd/numdata/numerical_data.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py index 573a935c..33697721 100644 --- a/src/compas_fd/numdata/numerical_data.py +++ b/src/compas_fd/numdata/numerical_data.py @@ -1,21 +1,9 @@ -from dataclasses import dataclass -from dataclasses import astuple - from compas_fd.result import Result -@dataclass class NumericalData: """Storage class for numerical arrays used by solver algorithms.""" - def __iter__(self): - return iter(astuple(self)) - - @classmethod - def from_params(cls, *args, **kwargs): - """Construct nuerical arrays from algorithm input parameters.""" - raise NotImplementedError - @classmethod def from_mesh(cls, mesh): """Construct numerical arrays from input mesh.""" From bcc35acacffac05b5a0eb779126265bed9d8bee7 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:00:56 +0100 Subject: [PATCH 17/26] Refactor FDNumericalData Change FDNumericalData from dataclass into regular class with lazy calculation of property variables. Implicit recalculation of dependent variables. Add shorthand aliases to numerical arrays. --- src/compas_fd/numdata/fd_numerical_data.py | 188 ++++++++++++++++----- 1 file changed, 147 insertions(+), 41 deletions(-) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index 21020d28..b478a38e 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -7,66 +7,172 @@ from typing_extensions import Annotated from nptyping import NDArray -from dataclasses import dataclass - 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 compas.numerical import normrow from compas_fd.result import Result from .numerical_data import NumericalData -@dataclass class FDNumericalData(NumericalData): """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 - @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) + 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.connectivity_matrix = edges + self.force_densities = force_densities + self.loads = loads @classmethod def from_mesh(cls, mesh): """Construct numerical arrays from input mesh.""" - raise NotImplementedError + 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'))) + return cls(vertices, fixed, edges, force_densities, loads) 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): + return self._xyz + + @xyz.setter + def xyz(self, vertices): + self._xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) + self._lengths = None + self._residuals = None + + @property + def fixed(self): + return self._fixed + + @fixed.setter + def fixed(self, fixed): + self._fixed = fixed + self._free = list(set(range(len(self.xyz))) - set(fixed)) + + @property + def free(self): + return self._free + + @property + def connectivity_matrix(self): + return self._connectivity_matrix + + @connectivity_matrix.setter + def connectivity_matrix(self, edges): + self._connectivity_matrix = connectivity_matrix(edges, 'csr') + self._connectivity_matrix_free = None + self._connectivity_matrix_fixed = None + + @property + def connectivity_matrix_free(self): + if self._connectivity_matrix_free is None: + self._connectivity_matrix_free = self.connectivity_matrix[:, self.free] + return self._connectivity_matrix_free + + @property + def connectivity_matrix_fixed(self): + if self._connectivity_matrix_fixed is None: + self._connectivity_matrix_fixed = self.connectivity_matrix[:, self.fixed] + return self._connectivity_matrix_fixed + + @property + def force_densities(self): + return self._force_densities + + @force_densities.setter + def force_densities(self, force_densities): + 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 + def force_densities_matrix(self): + if self._force_densities_matrix is None: + self._force_densities_matrix = diags( + [self.force_densities.flatten()], [0]) + return self._force_densities_matrix + + @property + def loads(self): + return self._loads + + @loads.setter + def loads(self, loads): + self._loads = (zeros_like(self.xyz) if loads is None else + asarray(loads, dtype=float64).reshape((-1, 3))) + self._residuals = None + + @property + def stifness_matrix(self): + if self._stifness_matrix is None: + self._stifness_matrix = self.C.T.dot(self.Q).dot(self.C) + return self._stifness_matrix + + @property + def stifness_matrix_free(self): + if self._stifness_matrix_free is None: + self._stifness_matrix_free = self.Ci.T.dot(self.Q).dot(self.Ci) + return self._stifness_matrix_free + + @property + def stifness_matrix_fixed(self): + if self._stifness_matrix_fixed is None: + self._stifness_matrix_fixed = self.Ci.T.dot(self.Q).dot(self.Cf) + return self._stifness_matrix_fixed + + @property + def forces(self): + if self._forces is None: + self._forces = self.q * self.ls + return self._forces + + @property + def lengths(self): + if self._lengths is None: + self._lengths = normrow(self.C.dot(self.xyz)) + return self._lengths + + @property + def residuals(self): + if self._residuals is None: + self._residuals = self.p - self.D.dot(self.xyz) + return self._residuals + + # aliases + C = connectivity_matrix + Ci = connectivity_matrix_free + Cf = connectivity_matrix_fixed + 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 From ce6bd53ed52dbc5cc7b50e4b9edef1abce91ad9f Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:01:30 +0100 Subject: [PATCH 18/26] Naming --- src/compas_fd/datastructures/cablemesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/compas_fd/datastructures/cablemesh.py b/src/compas_fd/datastructures/cablemesh.py index 7cfbfb2d..2d4d1811 100644 --- a/src/compas_fd/datastructures/cablemesh.py +++ b/src/compas_fd/datastructures/cablemesh.py @@ -13,9 +13,9 @@ class CableMesh(Mesh): if not compas.IPY: from compas_fd.fd import mesh_fd_numpy - from compas_fd.fd import mesh_fd_constrained_numpy + from compas_fd.fd import mesh_fd_constraint_numpy fd_numpy = mesh_fd_numpy - fd_constrained_numpy = mesh_fd_constrained_numpy + fd_constraint_numpy = mesh_fd_constraint_numpy def __init__(self): super(CableMesh, self).__init__() From 9defa6569b958c33e8b54cc0594af2b47e9a6aff Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:01:53 +0100 Subject: [PATCH 19/26] Naming --- tests/test_hypar.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_hypar.py b/tests/test_hypar.py index 29e06962..0238e49c 100644 --- a/tests/test_hypar.py +++ b/tests/test_hypar.py @@ -10,6 +10,6 @@ def test_hypar(): d = [0, 1, 0] vertices = [o, a, b, c, d] edges = [(0, 1), (0, 2), (0, 3), (0, 4)] - forcedensities = [1, 1, 1, 1] - result = fd_numpy(vertices=vertices, fixed=[1, 2, 3, 4], edges=edges, forcedensities=forcedensities) + force_densities = [1, 1, 1, 1] + result = fd_numpy(vertices=vertices, fixed=[1, 2, 3, 4], edges=edges, force_densities=force_densities) assert np.allclose(result.vertices[0], [0.5, 0.5, 0.5]) From 95135b6a6b8dae78ed6b68fa4a6ddb13ab7dcba3 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:03:58 +0100 Subject: [PATCH 20/26] Update fd functions Abstract fd convenience functions, by implementing Solver and FDNumericalData backends. --- src/compas_fd/fd/fd_constraint_numpy.py | 36 +++++++++++ src/compas_fd/fd/fd_iter_numpy.py | 35 ----------- src/compas_fd/fd/fd_numpy.py | 35 +++-------- src/compas_fd/fd/mesh_fd_constrained_numpy.py | 60 ------------------- ...r_numpy.py => mesh_fd_constraint_numpy.py} | 21 ++----- src/compas_fd/fd/mesh_fd_numpy.py | 31 ++++------ 6 files changed, 59 insertions(+), 159 deletions(-) create mode 100644 src/compas_fd/fd/fd_constraint_numpy.py delete mode 100644 src/compas_fd/fd/fd_iter_numpy.py delete mode 100644 src/compas_fd/fd/mesh_fd_constrained_numpy.py rename src/compas_fd/fd/{mesh_fd_iter_numpy.py => mesh_fd_constraint_numpy.py} (56%) diff --git a/src/compas_fd/fd/fd_constraint_numpy.py b/src/compas_fd/fd/fd_constraint_numpy.py new file mode 100644 index 00000000..621d3994 --- /dev/null +++ b/src/compas_fd/fd/fd_constraint_numpy.py @@ -0,0 +1,36 @@ +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 float64 + +from compas_fd.result import Result +from compas_fd.constraints import Constraint +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDConstraintSolver + + +def fd_constraint_numpy(*, + 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, + constraints: Sequence[Constraint], + kmax: int = 100, + tol_res: 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(vertices, fixed, edges, force_densities, loads) + solver = FDConstraintSolver(numdata, constraints, + kmax=100, tol_res=1E-3, tol_disp=1E-3) + result = solver() + return result diff --git a/src/compas_fd/fd/fd_iter_numpy.py b/src/compas_fd/fd/fd_iter_numpy.py deleted file mode 100644 index 7dda9752..00000000 --- a/src/compas_fd/fd/fd_iter_numpy.py +++ /dev/null @@ -1,35 +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 numpy import float64 - -from compas_fd.result import Result -from compas_fd.constraints import Constraint -from compas_fd.numdata import FDNumericalData -from compas_fd.solvers import FDConstraintSolver - - -def fd_iter_numpy(*, - 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, - constraints: Sequence[Constraint], - kmax: int = 100, - tol_res: 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, kmax, tol_res, tol_disp) - result = solver() - return result diff --git a/src/compas_fd/fd/fd_numpy.py b/src/compas_fd/fd/fd_numpy.py index a9da8ba3..fc9f888e 100644 --- a/src/compas_fd/fd/fd_numpy.py +++ b/src/compas_fd/fd/fd_numpy.py @@ -7,45 +7,24 @@ from typing_extensions import Annotated from nptyping import NDArray -from numpy import asarray, zeros_like from numpy import float64 -from scipy.sparse import diags -from scipy.sparse.linalg import spsolve - -from compas.numerical import connectivity_matrix -from compas.numerical import normrow from compas_fd.result import Result +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDSolver def fd_numpy(*, vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], fixed: List[int], edges: List[Tuple[int, int]], - forcedensities: List[float], + force_densities: List[float], loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None, ) -> Result: """ Compute the equilibrium coordinates of a system of vertices connected by edges. """ - v = len(vertices) - free = list(set(range(v)) - set(fixed)) - xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) - q = asarray(forcedensities, dtype=float64).reshape((-1, 1)) - if loads is None: - p = zeros_like(xyz) - else: - p = asarray(loads, dtype=float64).reshape((-1, 3)) - C = connectivity_matrix(edges, 'csr') - Ci = C[:, free] - Cf = C[:, fixed] - Ct = C.transpose() - Cit = Ci.transpose() - Q = diags([q.flatten()], [0]) - A = Cit.dot(Q).dot(Ci) - b = p[free] - Cit.dot(Q).dot(Cf).dot(xyz[fixed]) - xyz[free] = spsolve(A, b) - lengths = normrow(C.dot(xyz)) - forces = q * lengths - residuals = p - Ct.dot(Q).dot(C).dot(xyz) - return Result(xyz, residuals, forces, lengths) + numdata = FDNumericalData(vertices, fixed, edges, force_densities, loads) + solver = FDSolver(numdata) + result = solver() + return result diff --git a/src/compas_fd/fd/mesh_fd_constrained_numpy.py b/src/compas_fd/fd/mesh_fd_constrained_numpy.py deleted file mode 100644 index f98e3aa6..00000000 --- a/src/compas_fd/fd/mesh_fd_constrained_numpy.py +++ /dev/null @@ -1,60 +0,0 @@ -from numpy import array, asarray -from numpy import float64 - -import compas_fd -from .fd_constrained_numpy import fd_constrained_numpy - - -def mesh_fd_constrained_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': - """Iteratively find the equilibrium shape of a mesh for the given force densities. - - Parameters - ---------- - mesh : :class:`compas_fd.datastructures.CableMesh` - The mesh to equilibriate. - - Returns - ------- - :class:`compas_fd.datastructures.CableMesh` - The function updates the mesh in place, - but returns a reference to the updated mesh as well - for compatibility with RPCs. - - """ - v_i = mesh.vertex_index() - vertices = array(mesh.vertices_attributes('xyz'), dtype=float64) - 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})] - forcedensities = asarray([attr['q'] for edge, attr in mesh.edges_where({'_is_edge': True}, True)], dtype=float64).reshape((-1, 1)) - loads = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) - constraints = list(mesh.vertices_attribute('constraint')) - - result = fd_constrained_numpy(vertices=vertices, - fixed=fixed, - edges=edges, - forcedensities=forcedensities, - loads=loads, - constraints=constraints, - kmax=100, - tol_res=1E-3, - tol_disp=1E-3) - - _update_mesh(mesh, result) - - return mesh - - -def _update_mesh(mesh, result): - vertex_index = mesh.vertex_index() - for vertex, attr in mesh.vertices(True): - index = vertex_index[vertex] - attr['x'] = result.vertices[index, 0] - attr['y'] = result.vertices[index, 1] - attr['z'] = result.vertices[index, 2] - attr['_rx'] = result.residuals[index, 0] - attr['_ry'] = result.residuals[index, 1] - attr['_rz'] = result.residuals[index, 2] - - for index, (vertex, attr) in enumerate(mesh.edges_where({'_is_edge': True}, True)): - attr['_f'] = result.forces[index, 0] - attr['_l'] = result.lenghts[index, 0] diff --git a/src/compas_fd/fd/mesh_fd_iter_numpy.py b/src/compas_fd/fd/mesh_fd_constraint_numpy.py similarity index 56% rename from src/compas_fd/fd/mesh_fd_iter_numpy.py rename to src/compas_fd/fd/mesh_fd_constraint_numpy.py index 94113c2e..d450eddd 100644 --- a/src/compas_fd/fd/mesh_fd_iter_numpy.py +++ b/src/compas_fd/fd/mesh_fd_constraint_numpy.py @@ -1,13 +1,10 @@ -from numpy import array, asarray -from numpy import float64 - import compas_fd from compas_fd.numdata import FDNumericalData from compas_fd.solvers import FDConstraintSolver -def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': +def mesh_fd_constraint_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': """Iteratively find the equilibrium shape of a mesh for the given force densities. Parameters @@ -23,19 +20,11 @@ def mesh_fd_iter_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd for compatibility with RPCs. """ - k_i = mesh.key_index() - vertices = array(mesh.vertices_attributes('xyz'), dtype=float64) - fixed = [k_i[v] for v in mesh.vertices_where({'is_anchor': True})] - edges = [(k_i[u], k_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] - forcedensities = asarray([attr['q'] for key, attr in mesh.edges_where({'_is_edge': True}, True)], - dtype=float64).reshape((-1, 1)) - loads = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) - constraints = list(mesh.vertices_attribute('constraint')) - - numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) - solver = FDConstraintSolver(numdata, constraints, kmax=100, tol_res=1E-3, tol_disp=1E-3) + numdata = FDNumericalData.from_mesh(mesh) + constraints = list(c for c in mesh.vertices_attribute('constraint') if c) + solver = FDConstraintSolver(numdata, constraints, + kmax=100, tol_res=1E-3, tol_disp=1E-3) result = solver() - _update_mesh(mesh, result) return mesh diff --git a/src/compas_fd/fd/mesh_fd_numpy.py b/src/compas_fd/fd/mesh_fd_numpy.py index 316fad97..cef94115 100644 --- a/src/compas_fd/fd/mesh_fd_numpy.py +++ b/src/compas_fd/fd/mesh_fd_numpy.py @@ -1,9 +1,7 @@ -from numpy import array -from numpy import float64 - import compas_fd -from compas_fd.loads import SelfweightCalculator -from .fd_numpy import fd_numpy + +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDSolver def mesh_fd_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.datastructures.CableMesh': @@ -22,21 +20,16 @@ def mesh_fd_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.data for compatibility with RPCs. """ - k_i = mesh.key_index() - fixed = mesh.vertices_where({'is_anchor': True}) - fixed = [k_i[key] for key in fixed] - xyz = array(mesh.vertices_attributes('xyz'), dtype=float64) - p = array(mesh.vertices_attributes(('px', 'py', 'pz')), dtype=float64) - edges = [(k_i[u], k_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] - q = array([attr['q'] for key, attr in mesh.edges_where({'_is_edge': True}, True)], dtype=float64).reshape((-1, 1)) - density = mesh.attributes['density'] - calculate_sw = SelfweightCalculator(mesh, density=density) - p[:, 2] -= calculate_sw(xyz)[:, 0] - - result = fd_numpy(vertices=xyz, fixed=fixed, edges=edges, forcedensities=q, loads=p) + numdata = FDNumericalData.from_mesh(mesh) + solver = FDSolver(numdata) + result = solver() + _update_mesh(mesh, result) + return mesh + +def _update_mesh(mesh, result): for key, attr in mesh.vertices(True): - index = k_i[key] + index = mesh.key_index()[key] attr['x'] = result.vertices[index, 0] attr['y'] = result.vertices[index, 1] attr['z'] = result.vertices[index, 2] @@ -47,5 +40,3 @@ def mesh_fd_numpy(mesh: 'compas_fd.datastructures.CableMesh') -> 'compas_fd.data for index, (key, attr) in enumerate(mesh.edges_where({'_is_edge': True}, True)): attr['_f'] = result.forces[index, 0] attr['_l'] = result.lenghts[index, 0] - - return mesh From cf638b0d83f41713307e812ddc529bb482bcc567 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:04:12 +0100 Subject: [PATCH 21/26] Naming again --- src/compas_fd/fd/__init__.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/compas_fd/fd/__init__.py b/src/compas_fd/fd/__init__.py index d011d008..e6832bf7 100644 --- a/src/compas_fd/fd/__init__.py +++ b/src/compas_fd/fd/__init__.py @@ -25,16 +25,16 @@ if not compas.IPY: from .fd_numpy import fd_numpy - from .fd_constrained_numpy import fd_constrained_numpy + from .fd_constraint_numpy import fd_constraint_numpy from .mesh_fd_numpy import mesh_fd_numpy - from .mesh_fd_constrained_numpy import mesh_fd_constrained_numpy + from .mesh_fd_constraint_numpy import mesh_fd_constraint_numpy __all__ = [] if not compas.IPY: __all__ += [ 'fd_numpy', - 'fd_constrained_numpy', + 'fd_constraint_numpy', 'mesh_fd_numpy', - 'mesh_fd_constrained_numpy' + 'mesh_fd_constraint_numpy' ] From 89e4b5a995811d8d7de04576a7eb5255c962bcb3 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:04:48 +0100 Subject: [PATCH 22/26] Update test scripts Change scripts to call new implementations through Solvers and FDNumericalData. --- scripts/PLACEHOLDER | 0 scripts/test_constraints_curve.py | 33 ++++------------- scripts/test_constraints_frame.py | 31 ++++------------ ...onstraints.py => test_constraints_line.py} | 34 ++++-------------- scripts/test_constraints_plane.py | 36 +++++-------------- scripts/test_constraints_surface.py | 31 ++++------------ scripts/test_constraints_vector.py | 31 +++------------- scripts/test_fd_constraint_solver.py | 10 +++--- scripts/test_fd_solver.py | 6 ++-- 9 files changed, 48 insertions(+), 164 deletions(-) delete mode 100644 scripts/PLACEHOLDER rename scripts/{test_constraints.py => test_constraints_line.py} (56%) diff --git a/scripts/PLACEHOLDER b/scripts/PLACEHOLDER deleted file mode 100644 index e69de29b..00000000 diff --git a/scripts/test_constraints_curve.py b/scripts/test_constraints_curve.py index 5b3b4bed..94a13f5f 100644 --- a/scripts/test_constraints_curve.py +++ b/scripts/test_constraints_curve.py @@ -1,13 +1,11 @@ -from functools import partial import compas from compas.geometry import Vector, Point, Line from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy - from compas.geometry import Polyline, Bezier from compas.geometry import NurbsCurve from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -19,6 +17,7 @@ mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) mesh.vertices_attribute('t', 0.0) +# curve constraint vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] points = [Point(10, 0, 0), Point(6, 13, -4), Point(-3, 16, 3), Point(0, 20, 0)] bezier = Bezier(points) @@ -27,30 +26,12 @@ constraint = Constraint(curve) mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) - -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 - -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) +constraints = mesh.vertices_attribute('constraint') + +# solve +mesh = mesh_fd_constraint_numpy(mesh) + # ============================================================================== # Viz diff --git a/scripts/test_constraints_frame.py b/scripts/test_constraints_frame.py index afa0bdce..128f2717 100644 --- a/scripts/test_constraints_frame.py +++ b/scripts/test_constraints_frame.py @@ -1,10 +1,9 @@ -from functools import partial import compas from compas.geometry import Vector, Point, Line, Frame from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -22,30 +21,12 @@ constraint = Constraint(frame) mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) - -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 - -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) +constraints = mesh.vertices_attribute('constraint') + +# solve +mesh = mesh_fd_constraint_numpy(mesh) + # ============================================================================== # Viz diff --git a/scripts/test_constraints.py b/scripts/test_constraints_line.py similarity index 56% rename from scripts/test_constraints.py rename to scripts/test_constraints_line.py index 2e50763a..de4662a3 100644 --- a/scripts/test_constraints.py +++ b/scripts/test_constraints_line.py @@ -1,10 +1,9 @@ -from functools import partial import compas from compas.geometry import Vector, Point, Line from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -16,37 +15,18 @@ mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) mesh.vertices_attribute('t', 0.0) +# line constraint vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] line = Line(Point(10, 0, 0), Point(14, 10, 0)) constraint = Constraint(line) - mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() - -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) - -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 - -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) +constraints = mesh.vertices_attribute('constraint') + +# solve +mesh = mesh_fd_constraint_numpy(mesh) + # ============================================================================== # Viz diff --git a/scripts/test_constraints_plane.py b/scripts/test_constraints_plane.py index 3cce0c67..439a3bc9 100644 --- a/scripts/test_constraints_plane.py +++ b/scripts/test_constraints_plane.py @@ -1,10 +1,9 @@ -from functools import partial import compas from compas.geometry import Vector, Point, Line, Plane from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -16,35 +15,18 @@ mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) mesh.vertices_attribute('t', 0.0) -vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] -plane = Plane((15, 15, 0), (1, 0, 0)) +# plane constraint +vertex = list(mesh.vertices_where({'x': 0, 'y': 0}))[0] +plane = Plane((0, 0, 0), (1, 1, -1)) constraint = Constraint(plane) mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) - -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 - -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) +constraints = mesh.vertices_attribute('constraint') + +# solve +mesh = mesh_fd_constraint_numpy(mesh) + # ============================================================================== # Viz diff --git a/scripts/test_constraints_surface.py b/scripts/test_constraints_surface.py index 49bcbf41..06ccfd46 100644 --- a/scripts/test_constraints_surface.py +++ b/scripts/test_constraints_surface.py @@ -1,12 +1,11 @@ -from functools import partial import compas from compas.geometry import Vector, Point, Line from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy from compas.geometry import NurbsSurface from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -18,6 +17,7 @@ mesh.vertices_attribute('is_anchor', True, keys=list(mesh.vertices_where({'vertex_degree': 2}))) mesh.vertices_attribute('t', 0.0) +# surface constraint vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] points = [ [Point(0, 0, 10), Point(3, 0, 10), Point(6, 0, 10), Point(9, 0, 10)], @@ -29,30 +29,11 @@ constraint = Constraint(surface) mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) - -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 - -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) +constraints = mesh.vertices_attribute('constraint') + +# solve +mesh = mesh_fd_constraint_numpy(mesh) # ============================================================================== # Viz diff --git a/scripts/test_constraints_vector.py b/scripts/test_constraints_vector.py index 469e57db..07e743c1 100644 --- a/scripts/test_constraints_vector.py +++ b/scripts/test_constraints_vector.py @@ -1,10 +1,9 @@ -from functools import partial import compas -from compas.geometry import Vector, Point, Line, Vector +from compas.geometry import Point, Line, Vector from compas_fd.datastructures import CableMesh -from compas_fd.fd import fd_numpy from compas_fd.constraints import Constraint +from compas_fd.fd import mesh_fd_constraint_numpy from compas_view2.app import App from compas_view2.objects import Object, MeshObject @@ -20,34 +19,14 @@ vertex = list(mesh.vertices_where({'x': 10, 'y': 10}))[0] vector = Vector(0, -1, 1) constraint = Constraint(vector) - mesh.vertex_attribute(vertex, 'constraint', constraint) -vertex_index = mesh.vertex_index() - -vertices = mesh.vertices_attributes('xyz') -edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -loads = mesh.vertices_attributes(['px', 'py', 'pz']) -forcedensities = mesh.edges_attribute('q') fixed = list(mesh.vertices_where({'is_anchor': True})) -constraints = mesh.vertices_attribute('constraint', keys=fixed) - -fd = partial(fd_numpy, edges=edges, loads=loads, forcedensities=forcedensities, fixed=fixed) +constraints = mesh.vertices_attribute('constraint') -for k in range(100): - result = fd(vertices=vertices) - residuals = result.residuals - vertices = result.vertices - for vertex, constraint in zip(fixed, constraints): - if not constraint: - continue - constraint.location = vertices[vertex] - constraint.residual = residuals[vertex] - vertices[vertex] = constraint.location + constraint.tangent * 0.5 +# solve +mesh = mesh_fd_constraint_numpy(mesh) -for index, vertex in enumerate(mesh.vertices()): - mesh.vertex_attributes(vertex, 'xyz', vertices[index]) - mesh.vertex_attributes(vertex, ['_rx', '_ry', '_rz'], residuals[index]) # ============================================================================== # Viz diff --git a/scripts/test_fd_constraint_solver.py b/scripts/test_fd_constraint_solver.py index 692dbf8f..6b7cc21e 100644 --- a/scripts/test_fd_constraint_solver.py +++ b/scripts/test_fd_constraint_solver.py @@ -28,24 +28,24 @@ constraint = Constraint(plane) mesh.vertex_attribute(vertex, 'constraint', constraint) -# input solver parameters +# solver input parameters vertex_index = mesh.vertex_index() vertices = mesh.vertices_attributes('xyz') fixed = list(mesh.vertices_where({'is_anchor': True})) edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -forcedensities = mesh.edges_attribute('q') +force_densities = mesh.edges_attribute('q') loads = mesh.vertices_attributes(['px', 'py', 'pz']) constraints = list(mesh.vertices_attribute('constraint')) # set up iterative constraint solver -numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) -solver = FDConstraintSolver(numdata, constraints, kmax=30, +numdata = FDNumericalData(vertices, fixed, edges, force_densities, loads) +solver = FDConstraintSolver(numdata, constraints, kmax=5, tol_res=1E-3, tol_disp=1E-3) # run solver result = solver() -print(solver.iter_count) +print(solver.kcount) # update mesh for index, vertex in enumerate(mesh.vertices()): diff --git a/scripts/test_fd_solver.py b/scripts/test_fd_solver.py index e6e12860..c359c3ff 100644 --- a/scripts/test_fd_solver.py +++ b/scripts/test_fd_solver.py @@ -21,16 +21,16 @@ fixed = list(mesh.vertices_where({'is_anchor': True})) edges = [(vertex_index[u], vertex_index[v]) for u, v in mesh.edges_where({'_is_edge': True})] -forcedensities = mesh.edges_attribute('q') +force_densities = mesh.edges_attribute('q') loads = mesh.vertices_attributes(['px', 'py', 'pz']) # set up single iteration solver -numdata = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) +numdata = FDNumericalData(vertices, fixed, edges, force_densities, loads) solver = FDSolver(numdata) # run solver result = solver() -print(solver.iter_count) +print(solver.kcount) # update mesh for index, vertex in enumerate(mesh.vertices()): From ed695171a069ec825bdb18a8ac62a7d4fb47b103 Mon Sep 17 00:00:00 2001 From: Sam Date: Sun, 14 Nov 2021 18:08:45 +0100 Subject: [PATCH 23/26] Cleanup --- src/compas_fd/fd/fd_constrained_numpy.py | 101 ----------------------- 1 file changed, 101 deletions(-) delete mode 100644 src/compas_fd/fd/fd_constrained_numpy.py diff --git a/src/compas_fd/fd/fd_constrained_numpy.py b/src/compas_fd/fd/fd_constrained_numpy.py deleted file mode 100644 index b4272112..00000000 --- a/src/compas_fd/fd/fd_constrained_numpy.py +++ /dev/null @@ -1,101 +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 numpy import asarray -from numpy import float64 -from scipy.linalg import norm -from scipy.sparse.linalg import spsolve - -from compas.numerical import normrow - -from compas_fd.constraints import Constraint - -from compas_fd.numdata import FDNumericalData -from compas_fd.result import Result - - -def fd_constrained_numpy(*, - 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, - constraints: Sequence[Constraint], - kmax: int = 100, - tol_res: 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. - """ - nd = FDNumericalData.from_params(vertices, fixed, edges, forcedensities, loads) - - for k in range(kmax): - xyz_prev = nd.xyz - _solve_fd(nd) - _update_constraints(nd, constraints) - if (_is_converged_residuals(nd.tangent_residuals, tol_res) and - _is_converged_disp(xyz_prev, nd.xyz, tol_disp)): - break - - _post_process_fd(nd) - return nd.to_result() - - -def _solve_fd(numdata: FDNumericalData) -> None: - """Solve a single iteration for the equilibrium coordinates of a system. - All updated numerical arrays are stored in the numdata parameter. - """ - nd = numdata - b = nd.p[nd.free] - nd.Af.dot(nd.xyz[nd.fixed]) - nd.xyz[nd.free] = spsolve(nd.Ai, b) - numdata.residuals = nd.p - nd.A.dot(nd.xyz) - - -def _post_process_fd(numdata: FDNumericalData) -> None: - """Compute dependent numerical arrays from the numerical data after running solver. - All updated numerical arrays are stored in the numdata parameter. - """ - nd = numdata - nd.lengths = normrow(nd.C.dot(nd.xyz)) - numdata.forces = nd.q * nd.lengths - - -def _update_constraints(numdata: FDNumericalData, - constraints: Sequence[Constraint]) -> None: - """Update all vertex constraints by the residuals of the current iteration, - and store their updated vertex coordinates in the numdata parameter. - """ - nd = numdata - for vertex, constraint in enumerate(constraints): - if not constraint: - continue - 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 constraints if c]) - - -def _is_converged_residuals(residuals: NDArray[(Any, 3), float64], - tol_res: float) -> bool: - """Verify whether the maximum constraint residual is within tolerance.""" - if residuals is None or not residuals.any(): - return True - max_res = max(norm(residuals, axis=1)) - return max_res < tol_res - - -def _is_converged_disp(old_xyz: NDArray[(Any, 3), float64], - new_xyz: NDArray[(Any, 3), float64], - tol_disp: float) -> bool: - """Verify whether the maximum coordinate displacement - between consecutive iterations is within tolerance. - """ - max_dxyz = max(norm(new_xyz - old_xyz, axis=1)) - return max_dxyz < tol_disp From 48af8aa1e6c43be2c4e06c120133e038a74f2a56 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 09:41:41 +0100 Subject: [PATCH 24/26] Docstrings --- src/compas_fd/solvers/solver.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py index 3fb23e33..52e340c4 100644 --- a/src/compas_fd/solvers/solver.py +++ b/src/compas_fd/solvers/solver.py @@ -18,7 +18,7 @@ def __init__(self, self.result = None def __call__(self) -> Result: - """Iteratively apply the solve algorithm.""" + """Iteratively apply the solver algorithm.""" for self.kcount in range(1, self.kmax + 1): self.solve() if self.is_converged: @@ -28,12 +28,12 @@ def __call__(self) -> Result: return self.result def solve(self) -> None: - """Apply solver algorithm for a single iteration.""" + """Apply the solver algorithm for a single iteration.""" raise NotImplementedError @property def is_converged(self) -> bool: - """Verify if all convergence criteria are met.""" + """Verify if convergence criteria are met.""" raise NotImplementedError def post_process(self) -> None: From 4ff216f98441250ac903c6b35e55ac68b1a8aaaf Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 09:43:56 +0100 Subject: [PATCH 25/26] Store edge topology in numdata --- src/compas_fd/numdata/fd_numerical_data.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py index b478a38e..19e7234d 100644 --- a/src/compas_fd/numdata/fd_numerical_data.py +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -31,7 +31,7 @@ def __init__(self, ): self.xyz = vertices self.fixed = fixed - self.connectivity_matrix = edges + self.edges = edges self.force_densities = force_densities self.loads = loads @@ -75,15 +75,22 @@ def free(self): return self._free @property - def connectivity_matrix(self): - return self._connectivity_matrix + def edges(self): + return self._edges - @connectivity_matrix.setter - def connectivity_matrix(self, edges): - self._connectivity_matrix = connectivity_matrix(edges, 'csr') + @edges.setter + def edges(self, edges): + self._edges = edges + self._connectivity_matrix = None self._connectivity_matrix_free = None self._connectivity_matrix_fixed = None + @property + def connectivity_matrix(self): + if self._connectivity_matrix is None: + self._connectivity_matrix = connectivity_matrix(self.edges, 'csr') + return self._connectivity_matrix + @property def connectivity_matrix_free(self): if self._connectivity_matrix_free is None: From 206a278cd57809787d9ec87c5b1c6dd2db18fc55 Mon Sep 17 00:00:00 2001 From: Sam Date: Mon, 15 Nov 2021 14:08:47 +0100 Subject: [PATCH 26/26] Update fdsolver xyz explicit reset --- src/compas_fd/solvers/fd_solver.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/compas_fd/solvers/fd_solver.py b/src/compas_fd/solvers/fd_solver.py index 7371e89a..01229b60 100644 --- a/src/compas_fd/solvers/fd_solver.py +++ b/src/compas_fd/solvers/fd_solver.py @@ -17,6 +17,7 @@ def solve(self) -> None: nd = self.numdata b = nd.p[nd.free] - nd.Df.dot(nd.xyz[nd.fixed]) nd.xyz[nd.free] = spsolve(nd.Di, b) + nd.reset_xyz() @property def is_converged(self) -> bool: