Merge pull request #9 from Baraa-Elmoussa/main

Remove MatrixConstraint function and fix dependencies

Author: tomvanmele

CHANGELOG.md

@@ -13,6 +13,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+* Replaced kernel-layer `MatrixConstraint` with standard `pyo.Constraint` rules in `static_equilibrium_constraints` for compatibility with recent Pyomo versions.
+* Pinned Python to `< 3.14` for viewer compatibility.
+
 ### Removed
 
 

environment.yml

@@ -2,11 +2,12 @@ name: cra-dev
 channels:
   - conda-forge
 dependencies:
-  - python
+  - python < 3.14
   - pip
   - compas >=2.4
-  - ipopt ==3.14.9
-  - pyomo ==6.4.2
+  - ipopt
+  - pyomo 
   - pip:
       - compas_assembly >=0.7.0
+      - matplotlib
       - -e .[dev]

scripts/tutorial_cubes.py

@@ -1,22 +1,24 @@
 """This is the complete code for the tutorial"""
 
-from compas.geometry import Box, Frame, Translation, Rotation
+from compas.geometry import Box
+from compas.geometry import Frame
+from compas.geometry import Rotation
+from compas.geometry import Translation
 from compas_assembly.datastructures import Block
-from compas_cra.datastructures import CRA_Assembly
+
 from compas_cra.algorithms import assembly_interfaces_numpy
+from compas_cra.datastructures import CRA_Assembly
 from compas_cra.equilibrium import cra_solve
 from compas_cra.viewers import cra_view
 
-
-support = Box(Frame.worldXY(), 4, 2, 1)  # supporting block
+support = Box(frame=Frame.worldXY(), xsize=4, ysize=2, zsize=1)  # supporting block
 free1 = Box(
-    Frame.worldXY().transformed(
-        Translation.from_vector([0, 0, 1])
-        * Rotation.from_axis_and_angle([0, 0, 1], 0.2)
+    frame=Frame.worldXY().transformed(
+        Translation.from_vector([0, 0, 1]) * Rotation.from_axis_and_angle([0, 0, 1], 0.2)
     ),
-    1,
-    3,
-    1,
+    xsize=1,
+    ysize=3,
+    zsize=1,
 )  # block to analyse
 
 assembly = CRA_Assembly()  # create empty assembly
@@ -28,6 +30,6 @@
 
 assembly_interfaces_numpy(assembly)  # identify interface between two blocks
 
-cra_solve(assembly, verbose=True, timer=True)  # solve equilibrium using cra_solve
+cra_solve(assembly, verbose=False, timer=True)  # solve equilibrium using cra_solve
 
-cra_view(assembly, resultant=False, nodal=True, grid=True)  # visiualisation
+# cra_view(assembly, resultant=False, nodal=True, grid=True)  # visiualisation

src/compas_cra/algorithms/interfaces_numpy.py

@@ -2,6 +2,7 @@
 
 from math import fabs
 
+import compas.geometry as cg
 from compas.geometry import Frame
 from compas.geometry import local_to_world_coordinates_numpy
 from numpy import array
@@ -127,8 +128,8 @@ def assembly_interfaces_numpy(assembly, nmax=10, tmax=1e-6, amin=1e-1):
 
                         if area >= amin:
                             coords = [[x, y, 0.0] for x, y, z in intersection.exterior.coords]
-
-                            coords = local_to_world_coordinates_numpy(Frame(o, A[0], A[1]), coords)
+                            origin = cg.Point(*o.flatten())
+                            coords = local_to_world_coordinates_numpy(Frame(origin, A[0], A[1]), coords)
 
                             assembly.add_to_interfaces(
                                 node,

src/compas_cra/equilibrium/pyomo_helper.py

@@ -4,8 +4,6 @@
 from typing import Literal
 
 import pyomo.environ as pyo
-from numpy import zeros
-from pyomo.core.base.matrix_constraint import MatrixConstraint
 
 
 def initialisations(
@@ -244,18 +242,46 @@ def static_equilibrium_constraints(model, aeq, afr, p) -> Callable:
 
     """
 
-    equilibrium_constraints = MatrixConstraint(
-        aeq.data, aeq.indices, aeq.indptr, -p.flatten(), -p.flatten(), model.array_f
-    )
-
-    friction_constraint = MatrixConstraint(
-        afr.data,
-        afr.indices,
-        afr.indptr,
-        [None for i in range(afr.shape[0])],
-        zeros(afr.shape[0]),
-        model.array_f,
-    )
+    # The previous MatrixConstraint function was causing problems,
+    # it's now replaced by
+
+    # aeq = aeq.tocsr()  # tocsr gives data, indicies and indptr instead of decomposing them one by one
+    # afr = afr.tocsr()
+    rhs = (-p).flatten()
+    f_var = model.array_f
+
+    # Deprecated code for reference:
+    # equilibrium_constraints = matrix_constraint(
+    # aeq.data, aeq.indices, aeq.indptr, -p.flatten(),
+    # -p.flatten(), model.array_f )
+    # A_eq is the tangent matrix, equal to -p,
+    # A container for constraints of the form lb <= Ax <= ub.
+    # here, lp and ub are both equal to -p,
+    # so the constraint is Ax = -p, which is the equilibrium equation.
+    # https://pyomo.readthedocs.io/en/6.9.1/api/pyomo.core.kernel.matrix_constraint.matrix_constraint.html
+
+    def eq_rule(m, i):
+        s, e = aeq.indptr[i], aeq.indptr[i + 1]
+        if s == e:
+            return pyo.Constraint.Skip
+        expr = pyo.quicksum(float(aeq.data[k]) * f_var[int(aeq.indices[k])] for k in range(s, e))
+        return expr == float(rhs[i])
+
+    # Expression equivalent to:
+    # -p <= A_eq * f <= -p (Friction)
+
+    def fr_rule(m, j):
+        s, e = afr.indptr[j], afr.indptr[j + 1]
+        if s == e:
+            return pyo.Constraint.Skip
+        expr = pyo.quicksum(float(afr.data[k]) * f_var[int(afr.indices[k])] for k in range(s, e))
+        return expr <= 0.0
+
+    # Expression equivalent to:
+    # A_fr * f <= 0 (Unilateral constraint)
+
+    equilibrium_constraints = pyo.Constraint(range(aeq.shape[0]), rule=eq_rule)
+    friction_constraint = pyo.Constraint(range(afr.shape[0]), rule=fr_rule)
     return equilibrium_constraints, friction_constraint