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 34a4e80d..6b7cc21e 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 @@ -27,26 +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')) -# 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(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.kcount) # update mesh for index, vertex in enumerate(mesh.vertices()): diff --git a/scripts/test_fd_solver.py b/scripts/test_fd_solver.py new file mode 100644 index 00000000..c359c3ff --- /dev/null +++ b/scripts/test_fd_solver.py @@ -0,0 +1,57 @@ +import compas +from compas.geometry import Point, Line, Vector +from compas_fd.datastructures import CableMesh + +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 + +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) + +# 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})] +force_densities = mesh.edges_attribute('q') +loads = mesh.vertices_attributes(['px', 'py', 'pz']) + +# set up single iteration solver +numdata = FDNumericalData(vertices, fixed, edges, force_densities, loads) +solver = FDSolver(numdata) + +# run solver +result = solver() +print(solver.kcount) + +# 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 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/__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/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__() 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' ] 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 d2c3c032..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 .fd_numerical_data import FDNumericalData -from .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 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_numerical_data.py b/src/compas_fd/fd/fd_numerical_data.py deleted file mode 100644 index fb061df5..00000000 --- a/src/compas_fd/fd/fd_numerical_data.py +++ /dev/null @@ -1,74 +0,0 @@ -from typing import Any -from typing import Tuple -from typing import List -from typing import Union -from typing import Sequence -from typing import Optional -from typing_extensions import Annotated -from nptyping import NDArray - -from dataclasses import dataclass -from dataclasses import astuple - -from numpy import asarray -from numpy import float64 -from numpy import zeros_like -from scipy.sparse import diags - -from compas.numerical import connectivity_matrix - -from .result import Result - - -@dataclass -class FDNumericalData: - """Stores numerical data used by the force density algorithms.""" - free: int - fixed: int - xyz: NDArray[(Any, 3), float64] - C: NDArray[(Any, Any), int] - q: NDArray[(Any, 1), float64] - Q: NDArray[(Any, Any), float64] - p: NDArray[(Any, 1), float64] - A: NDArray[(Any, Any), float64] - Ai: NDArray[(Any, Any), float64] - Af: NDArray[(Any, Any), float64] - forces: NDArray[(Any, 1), float64] = None - lengths: NDArray[(Any, 1), float64] = None - residuals: NDArray[(Any, 3), float64] = None - tangent_residuals: NDArray[(Any, 3), float64] = None - normal_residuals: NDArray[(Any, 1), float64] = None - - def __iter__(self): - return iter(astuple(self)) - - @classmethod - def from_params(cls, - vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], - fixed: List[int], - edges: List[Tuple[int, int]], - forcedensities: List[float], - loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None): - """Construct numerical arrays from force density solver input parameters.""" - free = list(set(range(len(vertices))) - set(fixed)) - xyz = asarray(vertices, dtype=float64).reshape((-1, 3)) - C = connectivity_matrix(edges, 'csr') - Ci = C[:, free] - Cf = C[:, fixed] - q = asarray(forcedensities, dtype=float64).reshape((-1, 1)) - Q = diags([q.flatten()], [0]) - p = (zeros_like(xyz) if loads is None else - asarray(loads, dtype=float64).reshape((-1, 3))) - A = C.T.dot(Q).dot(C) - Ai = Ci.T.dot(Q).dot(Ci) - Af = Ci.T.dot(Q).dot(Cf) - return cls(free, fixed, xyz, C, q, Q, p, A, Ai, Af) - - @classmethod - def from_mesh(cls, mesh): - """Construct numerical arrays from input mesh.""" - raise NotImplementedError - - def to_result(self) -> Result: - """Parse relevant numerical data into a Result object.""" - return Result(self.xyz, self.residuals, self.forces, self.lengths) diff --git a/src/compas_fd/fd/fd_numpy.py b/src/compas_fd/fd/fd_numpy.py index 241d917d..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 .result import Result +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_constraint_numpy.py b/src/compas_fd/fd/mesh_fd_constraint_numpy.py new file mode 100644 index 00000000..d450eddd --- /dev/null +++ b/src/compas_fd/fd/mesh_fd_constraint_numpy.py @@ -0,0 +1,44 @@ +import compas_fd + +from compas_fd.numdata import FDNumericalData +from compas_fd.solvers import FDConstraintSolver + + +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 + ---------- + 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. + + """ + 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 + + +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] 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 diff --git a/src/compas_fd/numdata/__init__.py b/src/compas_fd/numdata/__init__.py new file mode 100644 index 00000000..036c517b --- /dev/null +++ b/src/compas_fd/numdata/__init__.py @@ -0,0 +1,33 @@ +""" +****************** +compas_fd.numdata +****************** + +.. currentmodule:: compas_fd.numdata + + +Classes +======= + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + NumericalData + FDNumericalData + +""" +from __future__ import print_function +from __future__ import absolute_import +from __future__ import division + +import compas +from .numerical_data import NumericalData + +if not compas.IPY: + from .fd_numerical_data import FDNumericalData + +__all__ = ['NumericalData'] + +if not compas.IPY: + __all__ += ['FDNumericalData'] diff --git a/src/compas_fd/numdata/fd_numerical_data.py b/src/compas_fd/numdata/fd_numerical_data.py new file mode 100644 index 00000000..19e7234d --- /dev/null +++ b/src/compas_fd/numdata/fd_numerical_data.py @@ -0,0 +1,185 @@ +from typing import Any +from typing import Tuple +from typing import List +from typing import Union +from typing import Sequence +from typing import Optional +from typing_extensions import Annotated +from nptyping import NDArray + +from numpy import asarray +from numpy import float64 +from numpy import zeros_like +from scipy.sparse import diags + +from compas.numerical import connectivity_matrix +from compas.numerical import normrow + +from compas_fd.result import Result +from .numerical_data import NumericalData + + +class FDNumericalData(NumericalData): + """Stores numerical data used by the force density algorithms.""" + + def __init__(self, + vertices: Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]], + fixed: List[int], + edges: List[Tuple[int, int]], + force_densities: List[float], + loads: Optional[Union[Sequence[Annotated[List[float], 3]], NDArray[(Any, 3), float64]]] = None + ): + self.xyz = vertices + self.fixed = fixed + self.edges = edges + self.force_densities = force_densities + self.loads = loads + + @classmethod + def from_mesh(cls, mesh): + """Construct numerical arrays from input mesh.""" + v_i = mesh.vertex_index() + vertices = asarray(mesh.vertices_attributes('xyz')) + fixed = [v_i[v] for v in mesh.vertices_where({'is_anchor': True})] + edges = [(v_i[u], v_i[v]) for u, v in mesh.edges_where({'_is_edge': True})] + force_densities = asarray([attr['q'] for key, attr in + mesh.edges_where({'_is_edge': True}, True)]) + loads = asarray(mesh.vertices_attributes(('px', 'py', 'pz'))) + 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 edges(self): + return self._edges + + @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: + 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 diff --git a/src/compas_fd/numdata/numerical_data.py b/src/compas_fd/numdata/numerical_data.py new file mode 100644 index 00000000..33697721 --- /dev/null +++ b/src/compas_fd/numdata/numerical_data.py @@ -0,0 +1,14 @@ +from compas_fd.result import Result + + +class NumericalData: + """Storage class for numerical arrays used by solver algorithms.""" + + @classmethod + def from_mesh(cls, mesh): + """Construct numerical arrays from input mesh.""" + raise NotImplementedError + + def to_result(self) -> Result: + """Parse relevant numerical data into a Result object.""" + raise NotImplementedError 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 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'] 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..7128c976 --- /dev/null +++ b/src/compas_fd/solvers/fd_constraint_solver.py @@ -0,0 +1,75 @@ +from typing import Sequence + +from numpy import asarray +from numpy.linalg import norm +from scipy.sparse.linalg import spsolve + +from compas_fd.numdata import FDNumericalData +from compas_fd.constraints import Constraint +from .solver import Solver + + +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, + kmax: int = 100, + tol_res: float = 1E-3, + tol_disp: float = 1E-3, + **kwargs + ) -> None: + super(FDConstraintSolver, self).__init__(numdata, kmax, **kwargs) + self.constraints = constraints + self.tol_res = tol_res + self.tol_disp = tol_disp + + def solve(self) -> None: + """Apply force density algorithm for a single iteration.""" + nd = self.numdata + 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) -> 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 = xyz[vertex] + constraint.residual = nd.residuals[vertex] + 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: + """Verify if all convergence criteria are met.""" + return (self._is_converged_residuals and + self._is_converged_displacements) + + @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_displacements(self) -> bool: + """Verify whether the maximum coordinate displacement + between consecutive iterations is within tolerance. + """ + nd = self.numdata + max_dxyz = max(norm(nd.xyz - nd.xyz_prev, axis=1)) + return max_dxyz < self.tol_disp diff --git a/src/compas_fd/solvers/fd_solver.py b/src/compas_fd/solvers/fd_solver.py new file mode 100644 index 00000000..01229b60 --- /dev/null +++ b/src/compas_fd/solvers/fd_solver.py @@ -0,0 +1,25 @@ +from scipy.sparse.linalg import spsolve + +from compas_fd.numdata import FDNumericalData +from .solver import Solver + + +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.""" + 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: + """Verify if all convergence criteria are met.""" + return True diff --git a/src/compas_fd/solvers/solver.py b/src/compas_fd/solvers/solver.py new file mode 100644 index 00000000..52e340c4 --- /dev/null +++ b/src/compas_fd/solvers/solver.py @@ -0,0 +1,41 @@ +from compas_fd.numdata import NumericalData +from compas_fd.result import Result + + +class Solver: + """Manager class for running iterative algorithms, + updating numerical data and storing results. + """ + + def __init__(self, + numdata: NumericalData, + kmax: int = 100, + **kwargs + ) -> None: + self.numdata = numdata + self.kmax = kmax + self.kcount = 0 + self.result = None + + def __call__(self) -> Result: + """Iteratively apply the solver algorithm.""" + for self.kcount in range(1, self.kmax + 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 the solver algorithm for a single iteration.""" + raise NotImplementedError + + @property + def is_converged(self) -> bool: + """Verify if convergence criteria are met.""" + raise NotImplementedError + + def post_process(self) -> None: + """Callable after ending the solver loop.""" + pass 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])