From 397885fcdfe0cb7ab5200b7f87299836aa7e19a5 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Wed, 30 Oct 2024 19:57:06 +0100 Subject: [PATCH 01/10] add variable for gravitational acceleration --- src/compas_cra/equilibrium/cra_helper.py | 6 ++++-- src/compas_cra/equilibrium/cra_penalty_pyomo.py | 5 ++++- src/compas_cra/equilibrium/cra_pyomo.py | 5 ++++- src/compas_cra/equilibrium/rbe_pyomo.py | 5 ++++- 4 files changed, 16 insertions(+), 5 deletions(-) diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index 24e91dee..3047084b 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 ------- @@ -84,7 +86,7 @@ def external_force_setup(assembly, density): 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) + p[index][2] = -block.volume() * (block.attributes["density"] if "density" in block.attributes else density) * gravity print((block.attributes["density"] if "density" in block.attributes else density)) p = np.array(p, dtype=float) 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) From cbe006048d8700b0d8b413a1d55e7b6d08ab5c1a Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Thu, 31 Oct 2024 08:51:51 +0100 Subject: [PATCH 02/10] add function to setup additional loads --- docs/api/compas_cra.equilibrium.rst | 1 + src/compas_cra/equilibrium/__init__.py | 2 ++ src/compas_cra/equilibrium/cra_helper.py | 22 ++++++++++++++++++++++ 3 files changed, 25 insertions(+) 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 3047084b..62fb7a03 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/src/compas_cra/equilibrium/cra_helper.py @@ -117,6 +117,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. From 22f7d3c256796c843d52c78374ece2045bc30965 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Thu, 31 Oct 2024 11:37:06 +0100 Subject: [PATCH 03/10] add additional loads to external_force_setup --- src/compas_cra/equilibrium/cra_helper.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index 62fb7a03..c3480019 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/src/compas_cra/equilibrium/cra_helper.py @@ -87,6 +87,7 @@ def external_force_setup(assembly, density, gravity): block = assembly.node_block(node) index = key_index[node] p[index][2] = -block.volume() * (block.attributes["density"] if "density" in block.attributes else density) * gravity + p[index][2] += -block.attributes["load"] if "load" in block.attributes else 0 print((block.attributes["density"] if "density" in block.attributes else density)) p = np.array(p, dtype=float) From c77465c7479d40087ee149126a9a701d3b8f4492 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Thu, 31 Oct 2024 13:42:25 +0100 Subject: [PATCH 04/10] clarify print output of external_force_setup --- src/compas_cra/equilibrium/cra_helper.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index c3480019..43347d4f 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/src/compas_cra/equilibrium/cra_helper.py @@ -82,13 +82,15 @@ def external_force_setup(assembly, density, gravity): num_nodes = assembly.graph.number_of_nodes() key_index = {key: index for index, key in enumerate(assembly.graph.nodes())} + print("resulting vertical force p on blocks:") + 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) * gravity p[index][2] += -block.attributes["load"] if "load" in block.attributes else 0 - print((block.attributes["density"] if "density" in block.attributes else density)) + print(f"{index}: {round(p[index][2],3)}") p = np.array(p, dtype=float) p = p[free, :].reshape((-1, 1), order="C") From b380b43a43c47a41afe9feb30e11d6d90330d58a Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Thu, 31 Oct 2024 14:39:34 +0100 Subject: [PATCH 05/10] adjust value of gravity in tests --- tests/test_cra.py | 2 +- tests/test_cra_penalty.py | 2 +- tests/test_rbe.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_cra.py b/tests/test_cra.py index 190ed62e..e6844074 100644 --- a/tests/test_cra.py +++ b/tests/test_cra.py @@ -25,7 +25,7 @@ def test_cra(): assembly.add_interfaces_from_meshes([interface1], 0, 1) - cra_solve(assembly, density=1) + cra_solve(assembly, density=1, gravity=1) IS_FORCE_CORRECT = True for edge in assembly.graph.edges(): diff --git a/tests/test_cra_penalty.py b/tests/test_cra_penalty.py index fb8ba65f..d2267b3e 100644 --- a/tests/test_cra_penalty.py +++ b/tests/test_cra_penalty.py @@ -18,7 +18,7 @@ def test_cra_penalty(): assembly_interfaces_numpy(assembly, amin=1e-6, tmax=1e-4) - cra_penalty_solve(assembly, density=1) + cra_penalty_solve(assembly, density=1, gravity=1) block = assembly.graph.node_attribute(1, "block") weight = block.volume() diff --git a/tests/test_rbe.py b/tests/test_rbe.py index f02c8d03..00c2abbc 100644 --- a/tests/test_rbe.py +++ b/tests/test_rbe.py @@ -18,7 +18,7 @@ def test_cra_penalty(): assembly_interfaces_numpy(assembly, amin=1e-6, tmax=1e-4) - rbe_solve(assembly, density=1) + rbe_solve(assembly, density=1, gravity=1) block = assembly.graph.node_attribute(1, "block") weight = block.volume() From 12a7b49bf040bc5ec77431d00087bbb9772ee777 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Fri, 1 Nov 2024 10:24:22 +0100 Subject: [PATCH 06/10] refactor code of external_force_setup --- src/compas_cra/equilibrium/cra_helper.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index 43347d4f..5dabd7bb 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/src/compas_cra/equilibrium/cra_helper.py @@ -82,15 +82,21 @@ def external_force_setup(assembly, density, gravity): num_nodes = assembly.graph.number_of_nodes() key_index = {key: index for index, key in enumerate(assembly.graph.nodes())} - print("resulting vertical force p on blocks:") + 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) * gravity - p[index][2] += -block.attributes["load"] if "load" in block.attributes else 0 - print(f"{index}: {round(p[index][2],3)}") + # determine weight and load + weight = block.volume() * (block.attributes["density"] if "density" in block.attributes else density) * gravity + load = block.attributes["load"] if "load" in block.attributes else 0 + # set external force + p[index][2] = -(weight + load) + print(f"{node:>4} {round(p[index][2],3):>10}") + + print() p = np.array(p, dtype=float) p = p[free, :].reshape((-1, 1), order="C") From 918096f119b07cdd263757a4aa7a321e5f187d66 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Sat, 2 Nov 2024 22:24:32 +0100 Subject: [PATCH 07/10] add changes to CHANGELOG.md --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 033f2440..0d3c6ffa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,9 +10,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added * Added instructions for fixing problems on Windows. +* Added function to set additional vertical loads. ### Changed +* The force calculation now works in SI-units. + ### Removed From f57c98b664ad372a6cc676236b4ba83221f5d33a Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Sun, 12 Jan 2025 17:53:12 +0100 Subject: [PATCH 08/10] revert external_force_setup to normalized forces --- src/compas_cra/equilibrium/cra_helper.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/compas_cra/equilibrium/cra_helper.py b/src/compas_cra/equilibrium/cra_helper.py index 5dabd7bb..3c8d0066 100644 --- a/src/compas_cra/equilibrium/cra_helper.py +++ b/src/compas_cra/equilibrium/cra_helper.py @@ -90,11 +90,11 @@ def external_force_setup(assembly, density, gravity): block = assembly.node_block(node) index = key_index[node] # determine weight and load - weight = block.volume() * (block.attributes["density"] if "density" in block.attributes else density) * gravity - load = block.attributes["load"] if "load" in block.attributes else 0 - # set external force + 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],3):>10}") + print(f"{node:>4} {round(p[index][2] * gravity, 3):>10}") print() From 5a8142cf630eb7b6896568616626ba12d3951fa8 Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Sun, 12 Jan 2025 18:33:12 +0100 Subject: [PATCH 09/10] restore tests without gravity --- tests/test_cra.py | 2 +- tests/test_cra_penalty.py | 2 +- tests/test_rbe.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_cra.py b/tests/test_cra.py index e6844074..190ed62e 100644 --- a/tests/test_cra.py +++ b/tests/test_cra.py @@ -25,7 +25,7 @@ def test_cra(): assembly.add_interfaces_from_meshes([interface1], 0, 1) - cra_solve(assembly, density=1, gravity=1) + cra_solve(assembly, density=1) IS_FORCE_CORRECT = True for edge in assembly.graph.edges(): diff --git a/tests/test_cra_penalty.py b/tests/test_cra_penalty.py index d2267b3e..fb8ba65f 100644 --- a/tests/test_cra_penalty.py +++ b/tests/test_cra_penalty.py @@ -18,7 +18,7 @@ def test_cra_penalty(): assembly_interfaces_numpy(assembly, amin=1e-6, tmax=1e-4) - cra_penalty_solve(assembly, density=1, gravity=1) + cra_penalty_solve(assembly, density=1) block = assembly.graph.node_attribute(1, "block") weight = block.volume() diff --git a/tests/test_rbe.py b/tests/test_rbe.py index 00c2abbc..f02c8d03 100644 --- a/tests/test_rbe.py +++ b/tests/test_rbe.py @@ -18,7 +18,7 @@ def test_cra_penalty(): assembly_interfaces_numpy(assembly, amin=1e-6, tmax=1e-4) - rbe_solve(assembly, density=1, gravity=1) + rbe_solve(assembly, density=1) block = assembly.graph.node_attribute(1, "block") weight = block.volume() From d14c559c3a714ffb427c85a183986029f3fc1bab Mon Sep 17 00:00:00 2001 From: Korbinian Reischl Date: Sun, 12 Jan 2025 18:33:54 +0100 Subject: [PATCH 10/10] update changelog --- CHANGELOG.md | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0d3c6ffa..a3328c16 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,12 +10,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added * Added instructions for fixing problems on Windows. -* Added function to set additional vertical loads. +* Added gravity variable and function to set additional vertical loads. ### Changed -* The force calculation now works in SI-units. - ### Removed