From bd07e5b7075a7ca1c0bbac26045dd80ffb2b60c4 Mon Sep 17 00:00:00 2001
From: tomvanmele <tom.v.mele@gmail.com>
Date: Wed, 28 Feb 2024 00:00:46 +0100
Subject: [PATCH] test script for the numpy version

---
 scripts/test_dr_numpy.py | 62 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 scripts/test_dr_numpy.py

diff --git a/scripts/test_dr_numpy.py b/scripts/test_dr_numpy.py
new file mode 100644
index 0000000..23b46a6
--- /dev/null
+++ b/scripts/test_dr_numpy.py
@@ -0,0 +1,62 @@
+from compas.colors import Color
+from compas.datastructures import Mesh
+from compas.geometry import Line
+from compas.geometry import Point
+from compas.geometry import Sphere
+from compas.geometry import Vector
+from compas_dr.dr_numpy import dr_numpy
+from compas_dr.numdata import InputData
+from compas_view2.app import App
+
+# =============================================================================
+# Input
+# =============================================================================
+
+mesh = Mesh.from_meshgrid(dx=10, nx=10)
+
+fixed = list(mesh.vertices_where(vertex_degree=2))
+loads = [[0, 0, 0]] * mesh.number_of_vertices()
+qpre = [1.0] * mesh.number_of_edges()
+
+for index, edge in enumerate(mesh.edges()):
+    if mesh.is_edge_on_boundary(edge):
+        qpre[index] = 10
+
+inputdata = InputData.from_mesh(mesh, fixed, loads, qpre)
+
+# =============================================================================
+# Solve
+# =============================================================================
+
+result = dr_numpy(inputdata, kmax=1000, rk_steps=2)
+
+# =============================================================================
+# Update
+# =============================================================================
+
+for vertex in mesh.vertices():
+    mesh.vertex_attributes(vertex, "xyz", result.xyz[vertex])
+
+# =============================================================================
+# Visualization
+# =============================================================================
+
+viewer = App()
+viewer.add(mesh)
+
+for vertex in mesh.vertices():
+    point = Point(*mesh.vertex_coordinates(vertex))
+    residual = Vector(*result.residuals[vertex])
+
+    if vertex in fixed:
+        ball = Sphere(radius=0.1, point=point)
+        viewer.add(ball.to_brep(), facecolor=Color.red())
+
+        viewer.add(
+            Line(point, point - residual * 0.1),
+            linecolor=Color.green().darkened(50),
+            linewidth=3,
+        )
+
+viewer.view.camera.zoom_extents()
+viewer.run()