Function to set additional loads

CHANGELOG.md

@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 * Added instructions for fixing problems on Windows.
+* Added gravity variable and function to set additional vertical loads.
 
 ### Changed
 

docs/api/compas_cra.equilibrium.rst

@@ -31,6 +31,7 @@ Equilibrium Helper Functions
     friction_setup
     external_force_setup
     density_setup
+    load_setup
     make_aeq
     make_afr
     unit_basis

src/compas_cra/equilibrium/__init__.py

@@ -6,6 +6,7 @@
     friction_setup,
     external_force_setup,
     density_setup,
+    load_setup,
     make_aeq,
     make_afr,
     unit_basis,
@@ -31,6 +32,7 @@
     "friction_setup",
     "external_force_setup",
     "density_setup",
+    "load_setup",
     "make_aeq",
     "make_afr",
     "unit_basis",

src/compas_cra/equilibrium/cra_helper.py

@@ -58,7 +58,7 @@ def friction_setup(assembly, mu, penalty=False, friction_net=False):
     return afr
 
 
-def external_force_setup(assembly, density):
+def external_force_setup(assembly, density, gravity):
     """Set up external force vector.
 
     Parameters
@@ -68,6 +68,8 @@ def external_force_setup(assembly, density):
     density : float
         Density of the material.
         If density attribute is not set, optimisation will use this density value.
+    gravity: float
+        Gravitational acceleration.
 
     Returns
     -------
@@ -80,12 +82,21 @@ def external_force_setup(assembly, density):
     num_nodes = assembly.graph.number_of_nodes()
     key_index = {key: index for index, key in enumerate(assembly.graph.nodes())}
 
+    print("\nblock   ext. force p")
+    print("-" * 20)
+
     p = [[0, 0, 0, 0, 0, 0] for i in range(num_nodes)]
     for node in assembly.graph.nodes():
         block = assembly.node_block(node)
         index = key_index[node]
-        p[index][2] = -block.volume() * (block.attributes["density"] if "density" in block.attributes else density)
-        print((block.attributes["density"] if "density" in block.attributes else density))
+        # determine weight and load
+        weight = block.volume() * (block.attributes["density"] if "density" in block.attributes else density)
+        load = (block.attributes["load"] if "load" in block.attributes else 0) / gravity  # normalized load
+        # set external force vector
+        p[index][2] = -(weight + load)
+        print(f"{node:>4}    {round(p[index][2] * gravity, 3):>10}")
+
+    print()
 
     p = np.array(p, dtype=float)
     p = p[free, :].reshape((-1, 1), order="C")
@@ -115,6 +126,28 @@ def density_setup(assembly, density):
             block.attributes["density"] = density[node]
 
 
+def load_setup(assembly, load):
+    """Set up additional load in vertical direction.
+
+    Parameters
+    ----------
+    assembly : :class:`~compas_assembly.datastructures.Assembly`
+        The rigid block assembly.
+    load : dict of float
+        load values, the dict key should match with assembly.graph.nodes()
+
+    Returns
+    -------
+    None
+
+    """
+
+    for node in assembly.graph.nodes():
+        block = assembly.graph.node_attribute(node, "block")
+        if node in load:
+            block.attributes["load"] = load[node]
+
+
 def num_free(assembly):
     """Return number of free blocks.
 

src/compas_cra/equilibrium/cra_penalty_pyomo.py

@@ -24,6 +24,7 @@ def cra_penalty_solve(
     assembly: Assembly,
     mu: float = 0.84,
     density: float = 1.0,
+    gravity: float = 9.81,
     d_bnd: float = 1e-3,
     eps: float = 1e-4,
     verbose: bool = False,
@@ -39,6 +40,8 @@ def cra_penalty_solve(
         Friction coefficient value.
     density : float, optional
         Density of the block material.
+    gravity : float, optional
+        Gravitational acceleration.
     d_bnd : float, optional
         The bound of virtual displacement d.
     eps : float, optional
@@ -107,7 +110,7 @@ def cra_penalty_solve(
     aeq = equilibrium_setup(assembly, penalty=False)
     aeq_b = equilibrium_setup(assembly, penalty=True)
     afr_b = friction_setup(assembly, mu, penalty=True)
-    p = external_force_setup(assembly, density)
+    p = external_force_setup(assembly, density, gravity)
 
     model.d = aeq.toarray().T @ model.array_q
     model.forces = f_basis * model.array_f[:, np.newaxis]  # force x in global coordinate

src/compas_cra/equilibrium/cra_pyomo.py

@@ -24,6 +24,7 @@ def cra_solve(
     assembly: Assembly,
     mu: float = 0.84,
     density: float = 1.0,
+    gravity: float = 9.81,
     d_bnd: float = 1e-3,
     eps: float = 1e-4,
     verbose: bool = False,
@@ -39,6 +40,8 @@ def cra_solve(
         Friction coefficient value.
     density : float, optional
         Density of the block material.
+    gravity : float, optional
+        Gravitational acceleration.
     d_bnd : float, optional
         The bound of virtual displacement d.
     eps : float, optional
@@ -101,7 +104,7 @@ def cra_solve(
 
     aeq = equilibrium_setup(assembly)
     afr = friction_setup(assembly, mu)
-    p = external_force_setup(assembly, density)
+    p = external_force_setup(assembly, density, gravity)
 
     model.d = aeq.toarray().T @ model.array_q
     model.forces = basis * model.array_f[:, np.newaxis]  # force x in global coordinate

src/compas_cra/equilibrium/rbe_pyomo.py

@@ -21,6 +21,7 @@ def rbe_solve(
     assembly: Assembly,
     mu: float = 0.84,
     density: float = 1.0,
+    gravity: float = 9.81,
     verbose: bool = False,
     timer: bool = False,
 ) -> Assembly:
@@ -34,6 +35,8 @@ def rbe_solve(
         Friction coefficient value.
     density : float, optional
         Density of the block material.
+    gravity : float, optional
+        Gravitational acceleration.
     verbose : bool, optional
         Print information during the execution of the algorithm.
     timer : bool, optional
@@ -78,7 +81,7 @@ def rbe_solve(
 
     aeq_b = equilibrium_setup(assembly, penalty=True)
     afr_b = friction_setup(assembly, mu, penalty=True)
-    p = external_force_setup(assembly, density)
+    p = external_force_setup(assembly, density, gravity)
 
     obj_rbe = objectives("rbe", (0, 1e0, 1e6, 1e0))
     eq_con, fr_con = static_equilibrium_constraints(model, aeq_b, afr_b, p)