From 2e4b95be71a411c39241f8d0eceb006aaa0bf5f8 Mon Sep 17 00:00:00 2001 From: Piotr Bartman-Szwarc Date: Wed, 26 Feb 2025 03:55:45 +0100 Subject: [PATCH] BBO:WIP: plots --- conmech/plotting/drawer.py | 4 + conmech/solvers/optimization/optimization.py | 3 +- examples/Bagirrov_Bartman_Ochal_2024.py | 78 ++++++++++---------- examples/Makela_et_al_1998.py | 66 ++++++++++++++++- 4 files changed, 109 insertions(+), 42 deletions(-) diff --git a/conmech/plotting/drawer.py b/conmech/plotting/drawer.py index c64c6a4a..0b327427 100644 --- a/conmech/plotting/drawer.py +++ b/conmech/plotting/drawer.py @@ -194,6 +194,10 @@ def draw_boundaries(self, axes): edges=self.mesh.neumann_boundary, nodes=nodes, axes=axes, edge_color="g" ) else: + nodes = self.state.displaced_nodes + self.draw_boundary( + edges=self.mesh.boundary_surfaces, nodes=nodes, axes=axes, edge_color="black" + ) self.draw_dirichlet(axes) def draw_dirichlet(self, axes): diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index f4be7511..f58d34f3 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -182,7 +182,8 @@ def _solve_impl( # pylint: disable=import-outside-toplevel,import-error) from kosopt import subgradient - _, nstep = method.rsplit(" ", 1) + # _, nstep = method.rsplit(" ", 1) + nstep = 5 try: nstep = int(nstep) except ValueError: diff --git a/examples/Bagirrov_Bartman_Ochal_2024.py b/examples/Bagirrov_Bartman_Ochal_2024.py index b2984bde..6ce30661 100644 --- a/examples/Bagirrov_Bartman_Ochal_2024.py +++ b/examples/Bagirrov_Bartman_Ochal_2024.py @@ -19,7 +19,7 @@ from examples.Makela_et_al_1998 import loss_value, plot_losses cc = 0.7 -mesh_density = 8 +mesh_density = 20 kN = 1000 mm = 0.001 E = cc * 1.378e8 * kN @@ -126,8 +126,9 @@ def plot_setup(config: Config, mesh_descr, contact): state = State(runner.body) drawer = Drawer(state=state, config=config) drawer.node_size = 0 + drawer.colorful = False drawer.original_mesh_color = None - drawer.deformed_mesh_color = "black" + drawer.deformed_mesh_color = "white" # drawer.normal_stress_scale = 10 drawer.field_name = None drawer.xlabel = "x" @@ -137,7 +138,7 @@ def plot_setup(config: Config, mesh_descr, contact): axes = (axes,) drawer.outer_forces_scale = 0.35 drawer.outer_forces_size = 5 - plt.title("Reference configuration") + # plt.title("Reference configuration") # to have nonzero force interface on Neumann boundary. # state.time = 4 drawer.x_min = 0 @@ -167,15 +168,14 @@ def plot_setup(config: Config, mesh_descr, contact): drawer.outer_forces_arrows = np.column_stack((f_x, f_y, d_x, d_y)) drawer.draw( fig_axes=(fig, axes[0]), - show=True, + show=False, # field_min=f_limits[0], # field_max=f_limits[1], - save=False, - + save="reference.pdf", ) -def main(config: Config, methods, forces, contact, prefix=""): +def main(config: Config, methods, forces, contact, prefix="", layers_num=None): """ Entrypoint to example. @@ -191,6 +191,7 @@ def main(config: Config, methods, forces, contact, prefix=""): try: path = f"{config.outputs_path}/{PREFIX}_mtd_{m}_frc_{f:.2e}" with open(path, "rb") as output: + print(path) _ = pickle.load(output) except IOError as ioe: print(ioe) @@ -259,7 +260,7 @@ def outer_forces(x, t=None): path = f"{config.outputs_path}/{PREFIX}_losses" # if config.show: - plot_losses(path) + plot_losses(path, slopes=layers_num) # for m in methods[:]: # for f in forces[3:4]: @@ -285,7 +286,7 @@ def outer_forces(x, t=None): # plt.show() -def composite_problem(config, layers_limits, thickness, methods, forces): +def composite_problem(config, layers_limits, thickness, methods, forces, layers_num): def bottom(x): if x >= thickness: return 0.0 @@ -298,57 +299,54 @@ def top(x): contact = make_composite_contact_law(layers_limits, alpha=bottom, beta=top) - # X = np.linspace(-thickness / 4, 5 / 4 * thickness, 1000) - # Y = np.empty(1000) - # for i in range(1000): - # Y[i] = contact.potential_normal_direction(X[i], 0.0, 0.0) - # plt.plot(X, Y) + X = np.linspace(-thickness / 4, 5 / 4 * thickness, 1000) + Y = np.empty(1000) + for i in range(1000): + Y[i] = contact.potential_normal_direction(X[i], 0.0, 0.0) + plt.plot(X, Y) + plt.show() + for i in range(1000): + Y[i] = contact.subderivative_normal_direction(X[i], 0.0, 0.0) + plt.plot(X, Y) # plt.show() - # for i in range(1000): - # Y[i] = contact.subderivative_normal_direction(X[i], 0.0, 0.0) - # plt.plot(X, Y) - # # plt.show() - # for i in range(1000): - # Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) - # plt.plot(X, Y) - # plt.show() - main(config, methods, forces, contact, prefix=str(len(layers_limits))) + for i in range(1000): + Y[i] = contact.sub2derivative_normal_direction(X[i], 0.0, 0.0) + plt.plot(X, Y) + plt.show() + main(config, methods, forces, contact, prefix=str(len(layers_limits)), layers_num=layers_num) def survey(config): methods = ( - # "gradiented BFGS", - # "gradiented CG", + "gradiented BFGS", + "gradiented CG", "BFGS", - # "CG", + "CG", "Powell", "qsm", - "globqsm 5", - "globqsm 8", - # "globqsm 10", - # "dc qsm", - # "dc globqsm" + "globqsm", )[:] forces = np.asarray(( # 15e3 * kN, - 14e3 * kN, - # 15e3 * kN, - # 15.5e3 * kN, 16e3 * kN, - # 16.5e3 * kN, - # 17e3 * kN, + 17e3 * kN, 18e3 * kN, - # 19e3 * kN, + 19e3 * kN, 20e3 * kN, - 22e3 * kN, - 24e3 * kN, + # 21e3 * kN, + 22.5e3 * kN, + # 23e3 * kN, + 25e3 * kN, + 30e3 * kN, + 35e3 * kN, + 40e3 * kN, )) for i in [0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12][:]: #range(1, 10 + 1, 2): thickness = 3 * mm layers_num = i # layers_limits = np.logspace(0, thickness, layers_num + 1) layers_limits = partition(0, thickness, layers_num + 1, p=1.25) - composite_problem(config, layers_limits, thickness, methods, forces) + composite_problem(config, layers_limits, thickness, methods, forces, layers_num) def partition(start, stop, num, p=1.0): diff --git a/examples/Makela_et_al_1998.py b/examples/Makela_et_al_1998.py index cdeb57ad..d07723a7 100644 --- a/examples/Makela_et_al_1998.py +++ b/examples/Makela_et_al_1998.py @@ -242,11 +242,75 @@ def outer_forces(x, t=None): plt.show() -def plot_losses(path): +def plot_losses(path, slopes=None): print("Plotting...") with open(path, "rb") as output: losses = pickle.load(output) + data = { + "BFGS": (losses["BFGS"], "0.7"), + "BFGS (subgrad)": (losses["gradiented BFGS"], "0.6"), + "CG": (losses["CG"], "0.4"), + "CG (subgrad)": (losses["gradiented CG"], "0.3"), + "Powell": (losses["Powell"], "blue"), + "subgradient": (losses["qsm"], "pink"), + "global subgradient": (losses["globqsm"], "red"), + } + + # Grid of plots + fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=True) + styles = ['-', '--', '-.', ':', '-.', '--', '-',] + + # loss + for i, (method, (results, color)) in enumerate(data.items()): + forces = np.array(list(results.keys())) / 1e6 # N/m2 to MPa + losses = -np.array([val[0] for val in results.values()]) + axes[0].plot(forces, losses, linestyle=styles[i % len(styles)], marker='o', label=method, + linewidth=2, color=color) + + axes[0].set_ylabel("$-\mathcal{L}(u)$", fontsize=14) + axes[0].set_xlabel("Load [MPa]", fontsize=14) + axes[0].grid(True, linestyle='--', alpha=0.6) + axes[0].legend(fontsize=10, loc="lower right") + # axes[0].set_xscale("log") + axes[0].set_yscale("log") + if slopes is not None: + fig.suptitle(f"Num. of layers: {slopes + 2}") + + for i, (method, (results, color)) in enumerate(data.items()): + forces = np.array(list(results.keys())) / 1e6 + times = np.array([val[1] for val in results.values()]) + + best_losses = {force: min(data[m][0][force][0] for m in data if force in data[m][0]) for force in + results.keys()} + + markers = [] + for force, (loss, _) in results.items(): + best_loss = best_losses[force] + if loss == best_loss: # best solution + marker = 's' + elif loss <= best_loss * 0.9 + 0.01: # near to solution + marker = 'o' + else: # fail + marker = 'x' + markers.append(marker) + + line, = axes[1].plot(forces, times, linestyle=styles[i % len(styles)], linewidth=2, + label=method, alpha=0.7, color=color) + color = line.get_color() # line color + for j, force in enumerate(forces): + axes[1].scatter(force, times[j], marker=markers[j], color=color, s=70) + + axes[1].set_ylabel("Computation Time [s]", fontsize=14) + axes[1].set_xlabel("Load [MPa]", fontsize=14) + axes[1].grid(True, linestyle='--', alpha=0.6) + axes[1].set_yscale("log") + # axes[1].legend(fontsize=10, loc="upper left") + + plt.tight_layout() + plt.savefig(path + ".pdf") + plt.show() + return for i, (mtd, values) in enumerate(losses.items()): forces_ = np.asarray(list(values.keys())) values_ = np.asarray(list(values.values()))[:, 0]