BBO:WIP: plots

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):

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:

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):

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]