BOST 2024

conmech/mesh/boundaries.py

@@ -102,7 +102,7 @@ def contact_boundary(self):
 
     @property
     def contact_normals(self):
-        return self.surface_normals[: self.contact_nodes_count]
+        return self.surface_normals[: len(self.boundaries["contact"].surfaces)]
 
     @property
     def neumann_boundary(self):

conmech/mesh/zoo/barochsoftar_2024.py

@@ -0,0 +1,45 @@
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2023  Piotr Bartman <[email protected]>
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 3
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,
+# USA.
+import dmsh
+
+from conmech.mesh.zoo.raw_mesh import RawMesh
+from conmech.properties.mesh_description import JOB2023MeshDescription
+
+
+class BOST2023(RawMesh):
+    def __init__(self, mesh_descr: JOB2023MeshDescription):
+        # pylint: disable=no-member  # for dmsh
+        geo = dmsh.Polygon(
+            [
+                [0.0, 0.2],
+                [0.2, 0.4],
+                [0.6, 0.4],
+                [0.8, 0.2],
+                [0.8, 0.0],
+                [1.2, 0.0],
+                [1.2, 1.2],
+                [0.8, 1.2],
+                [0.8, 1.0],
+                [0.6, 0.8],
+                [0.2, 0.8],
+                [0.0, 1.0],
+            ]
+        )
+        nodes, elements = dmsh.generate(geo, mesh_descr.max_element_perimeter)
+        super().__init__(nodes, elements)

conmech/plotting/drawer.py

@@ -84,28 +84,32 @@ def draw(
             fig.tight_layout()
             plt.show()
         if save:
-            self.save_plot(save_format)
+            self.save_plot(save_format, name=save)
 
     def set_axes_limits(self, axes, foundation):
         # pylint: disable=nested-min-max
         if self.x_min is None:
             self.x_min = min(
-                min(self.state.body.mesh.nodes[:, 0]), min(self.state.displaced_nodes[:, 0])
+                min(self.state.body.mesh.nodes[:, 0]),
+                min(self.state.displaced_nodes[:, 0]),
             )
         if self.x_max is None:
             self.x_max = max(
-                max(self.state.body.mesh.nodes[:, 0]), max(self.state.displaced_nodes[:, 0])
+                max(self.state.body.mesh.nodes[:, 0]),
+                max(self.state.displaced_nodes[:, 0]),
             )
         dx = self.x_max - self.x_min
         x_margin = dx * 0.2
         xlim = (self.x_min - x_margin, self.x_max + x_margin)
         if self.y_min is None:
             self.y_min = min(
-                min(self.state.body.mesh.nodes[:, 1]), min(self.state.displaced_nodes[:, 1])
+                min(self.state.body.mesh.nodes[:, 1]),
+                min(self.state.displaced_nodes[:, 1]),
             )
         if self.y_max is None:
             self.y_max = max(
-                max(self.state.body.mesh.nodes[:, 1]), max(self.state.displaced_nodes[:, 1])
+                max(self.state.body.mesh.nodes[:, 1]),
+                max(self.state.displaced_nodes[:, 1]),
             )
         dy = self.y_max - self.y_min
         y_margin = dy * 0.2
@@ -163,7 +167,10 @@ def draw_boundaries(self, axes):
                 edges=self.mesh.contact_boundary, nodes=nodes, axes=axes, edge_color="b"
             )
             self.draw_boundary(
-                edges=self.mesh.dirichlet_boundary, nodes=nodes, axes=axes, edge_color="r"
+                edges=self.mesh.dirichlet_boundary,
+                nodes=nodes,
+                axes=axes,
+                edge_color="r",
             )
             self.draw_boundary(
                 edges=self.mesh.neumann_boundary, nodes=nodes, axes=axes, edge_color="g"
@@ -175,7 +182,7 @@ def draw_dirichlet(self, axes):
         dirichlet_nodes = self.state.body.mesh.dirichlet_boundary
         dirichlet_nodes = list(set(dirichlet_nodes.flatten()))
         x = self.state.displaced_nodes[dirichlet_nodes][[0, 1]]
-        v = np.zeros_like(x) + np.asarray([0, 1])
+        v = np.zeros_like(x) + np.asarray([1, 0])
         if any(v[:, 0]) or any(v[:, 1]):  # to avoid warning
             axes.quiver(
                 x[:, 0],
@@ -184,7 +191,7 @@ def draw_dirichlet(self, axes):
                 v[:, 1],
                 angles="xy",
                 scale_units="xy",
-                scale=2.5,
+                scale=8.0,
                 headlength=10,
                 headaxislength=10,
                 headwidth=10,
@@ -245,7 +252,7 @@ def draw_stress(self, axes):
     def get_output_path(config, format_, name):
         output_dir = config.output_dir or str(config.timestamp)
         directory = f"{config.outputs_path}/{output_dir}"
-        name = name if name else config.timestamp
+        name = name if isinstance(name, str) else config.timestamp
         path = f"{directory}/{name}.{format_}"
         return path
 
@@ -261,7 +268,9 @@ def save_plot(self, format_, name=None):
         )
         plt.close()
 
-    def draw_boundary(self, edges, nodes, axes, label="", node_color="k", edge_color="k"):
+    def draw_boundary(
+        self, edges, nodes, axes, label="", node_color="k", edge_color="k"
+    ):
         graph = nx.Graph()
         for edge in edges:
             graph.add_edge(edge[0], edge[1])
@@ -303,6 +312,8 @@ def draw_field(self, field, v_min, v_max, axes, fig):
         # from mpl_toolkits.axes_grid1 import make_axes_locatable
         # divider = make_axes_locatable(axes)
         # cax = divider.append_axes("bottom", size="5%", pad=0.15)
-        sm = plt.cm.ScalarMappable(cmap=self.cmap, norm=plt.Normalize(vmin=v_min, vmax=v_max))
+        sm = plt.cm.ScalarMappable(
+            cmap=self.cmap, norm=plt.Normalize(vmin=v_min, vmax=v_max)
+        )
         sm.set_array([])
         fig.colorbar(sm, orientation="horizontal", label=self.field_label, ax=axes)

conmech/properties/mesh_description.py

@@ -89,6 +89,14 @@ def build(self):
         return Barboteu2008(self)
 
 
+@dataclass
+class BOST2023MeshDescription(GeneratedMeshDescription):
+    def build(self):
+        from conmech.mesh.zoo.barochsoftar_2024 import BOST2023
+
+        return BOST2023(self)
+
+
 @dataclass
 class JOB2023MeshDescription(GeneratedMeshDescription):
     def build(self):

conmech/solvers/optimization/optimization.py

@@ -93,8 +93,8 @@ def _solve_impl(
         displacement = np.squeeze(displacement.copy().reshape(1, -1))
         old_solution = np.squeeze(initial_guess.copy().reshape(1, -1))
         disp = kwargs.get("disp", False)
-        maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9))
-        tol = kwargs.get("tol", 1e-12)
+        maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e12))
+        tol = kwargs.get("tol", 1e-18)
         args = (
             self.body.mesh.nodes,
             self.body.mesh.contact_boundary,
@@ -105,6 +105,13 @@ def _solve_impl(
             self.time_step,
         )
 
+        loss = []
+        sols = []
+        sols.append(solution)
+        loss.append(self.loss(solution, *args)[0])
+        print("pre", method, self.loss(solution, *args))
+
+
         while norm >= fixed_point_abs_tol:
             if method.lower() in (  # TODO
                 "quasi secant method",
@@ -117,6 +124,34 @@ def _solve_impl(
                 from kosopt import qsmlm
 
                 solution = qsmlm.minimize(self.loss, solution, args=args, maxiter=maxiter)
+                sols[self.loss(solution, *args)] = solution
+            elif method.lower() == 'constrained':
+                import numba
+                contact_nodes_count = self.body.mesh.boundaries.contact_nodes_count
+                GAP = 0.0
+                @numba.njit()
+                def constr(x):
+                    offset = len(x) // 2
+                    t = x[offset:offset+contact_nodes_count]
+                    return np.min(t + GAP)
+
+                maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9))
+                tol = kwargs.get("tol", 1e-12)
+                result = scipy.optimize.minimize(
+                    self.loss,
+                    solution,
+                    args=args,
+                    options={"disp": disp, "maxiter": maxiter},
+                    tol=tol,
+                    constraints=({'type': 'ineq',
+                                 'fun': constr
+                                 })
+                )
+                solution = result.x
+                sols.append(solution)
+                loss.append(self.loss(solution, *args)[0])
+                print(method, self.loss(solution, *args))
+                break
             else:
                 result = scipy.optimize.minimize(
                     self.loss,
@@ -127,6 +162,14 @@ def _solve_impl(
                     tol=tol,
                 )
                 solution = result.x
+                sols.append(solution.copy())
+                loss.append(self.loss(solution, *args)[0])
+                print(method, self.loss(solution, *args))
+                break
+
             norm = np.linalg.norm(np.subtract(solution, old_solution))
             old_solution = solution.copy()
+        min_index = loss.index(np.min(loss))
+        solution = sols[min_index]
+        print(loss, "final:", self.loss(solution, *args))
         return solution