From 4a61d630c724a118d23d5a71034f2130492f89c1 Mon Sep 17 00:00:00 2001
From: Sam <Sam_Bouten@hotmail.com>
Date: Thu, 11 Nov 2021 15:21:36 +0100
Subject: [PATCH] Update test scripts

---
 ...solver.py => test_fd_constraint_solver.py} | 23 ++++----
 scripts/test_fd_solver.py                     | 56 +++++++++++++++++++
 2 files changed, 67 insertions(+), 12 deletions(-)
 rename scripts/{test_iterative_solver.py => test_fd_constraint_solver.py} (82%)
 create mode 100644 scripts/test_fd_solver.py

diff --git a/scripts/test_iterative_solver.py b/scripts/test_fd_constraint_solver.py
similarity index 82%
rename from scripts/test_iterative_solver.py
rename to scripts/test_fd_constraint_solver.py
index 9638d5b6..3ee17c7a 100644
--- a/scripts/test_iterative_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_iter_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
@@ -23,7 +24,7 @@
 
 # plane constraint
 vertex = list(mesh.vertices_where({'x': 0, 'y': 0}))[0]
-plane = Plane((-1, 0, 0), (2, 1, 1))
+plane = Plane((-1, 0, 0), (2, 1, -1))
 constraint = Constraint(plane)
 mesh.vertex_attribute(vertex, 'constraint', constraint)
 
@@ -37,16 +38,14 @@
 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 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_fd_solver.py b/scripts/test_fd_solver.py
new file mode 100644
index 00000000..ad177c1f
--- /dev/null
+++ b/scripts/test_fd_solver.py
@@ -0,0 +1,56 @@
+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})]
+forcedensities = 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)
+solver = FDSolver(numdata)
+
+# run solver
+result = solver()
+
+# 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()