Solver

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

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

scripts/test_constraints.py → 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

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

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

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

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()):