diff --git a/CHANGELOG.md b/CHANGELOG.md index 033f2440..a3328c16 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/api/compas_cra.equilibrium.rst b/docs/api/compas_cra.equilibrium.rst index b16f0a99..4afb5330 100644 --- a/docs/api/compas_cra.equilibrium.rst +++ b/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 diff --git a/src/compas_cra/equilibrium/__init__.py b/src/compas_cra/equilibrium/__init__.py index 0d00ca3f..904260cd 100644 --- a/src/compas_cra/equilibrium/__init__.py +++ b/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", diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index 24e91dee..3c8d0066 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/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. diff --git a/src/compas_cra/equilibrium/cra_penalty_pyomo.py b/src/compas_cra/equilibrium/cra_penalty_pyomo.py index 205972ce..9cb41755 100644 --- a/src/compas_cra/equilibrium/cra_penalty_pyomo.py +++ b/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 diff --git a/src/compas_cra/equilibrium/cra_pyomo.py b/src/compas_cra/equilibrium/cra_pyomo.py index c54924e8..b44b7076 100644 --- a/src/compas_cra/equilibrium/cra_pyomo.py +++ b/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 diff --git a/src/compas_cra/equilibrium/rbe_pyomo.py b/src/compas_cra/equilibrium/rbe_pyomo.py index 8660bfdc..47e10778 100644 --- a/src/compas_cra/equilibrium/rbe_pyomo.py +++ b/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)