cleanup

conmech/plotting/drawer.py

@@ -4,6 +4,7 @@
 @author: Michał Jureczka
 @author: Piotr Bartman
 """
+
 import time
 
 import matplotlib.pyplot as plt

conmech/solvers/optimization/optimization.py

@@ -1,6 +1,21 @@
-"""
-Created at 18.02.2021
-"""
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2021-2024  Piotr Bartman-Szwarc <[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 math
 from typing import Optional
 
@@ -17,7 +32,10 @@
 from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw
 from conmech.solvers.solver import Solver
 from conmech.solvers.solver_methods import (
-    make_cost_functional, make_equation, make_cost_functional_subgradient)
+    make_cost_functional,
+    make_equation,
+    make_cost_functional_subgradient,
+)
 
 
 class Optimization(Solver):
@@ -121,22 +139,21 @@ def _solve_impl(
             ):
                 # pylint: disable=import-outside-toplevel,import-error)
                 from kosopt import qsmlm
+
                 solution = qsmlm.minimize(
-                    self.loss, solution, args=args, maxiter=maxiter,
-                    subgradient=self.subgradient
+                    self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient
                 )
-            elif method.lower() in (
-                "subgradient"
-            ):
+            elif method.lower() in ("subgradient",):
+                # pylint: disable=import-outside-toplevel,import-error)
                 from kosopt import subgradient
+
                 solution = subgradient.minimize(
-                    self.loss, solution, args=args, maxiter=maxiter,
-                    subgradient=self.subgradient
+                    self.loss, solution, args=args, maxiter=maxiter, subgradient=self.subgradient
                 )
             elif method.lower() in (  # TODO
-                    "discontinuous gradient",
-                    "discontinuous gradient method",
-                    "dg",
+                "discontinuous gradient",
+                "discontinuous gradient method",
+                "dg",
             ):
                 # pylint: disable=import-outside-toplevel,import-error)
                 from kosopt import qsmlmi

conmech/solvers/solver.py

@@ -1,11 +1,25 @@
-"""
-Created at 18.02.2021
-"""
-
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2021-2024  Piotr Bartman-Szwarc <[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 time
 from typing import Optional
 
 import numpy as np
-import time
 
 from conmech.dynamics.statement import Statement, Variables
 from conmech.dynamics.contact.contact_law import ContactLaw

conmech/solvers/solver_methods.py

@@ -1,7 +1,21 @@
-"""
-Created at 21.08.2019
-"""
-
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2019-2024  Piotr Bartman-Szwarc <[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.
 from typing import Callable, Optional, Any
 
 import numba
@@ -211,6 +225,7 @@ def cost_functional(
 
     return cost_functional
 
+
 def make_cost_functional_subgradient(
     djn: Callable, djt: Optional[Callable] = None, dh_functional: Optional[Callable] = None
 ):
@@ -220,8 +235,7 @@ def make_cost_functional_subgradient(
     DIMENSION = 2
 
     @numba.njit()
-    def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary,
-                            contact_normals):
+    def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary, contact_normals):
         cost = np.zeros_like(u_vector)
         offset = len(u_vector) // DIMENSION
 
@@ -232,11 +246,11 @@ def contact_subgradient(u_vector, u_vector_old, nodes, contact_boundary,
             n_1 = nodes[n_id_1]
             if n_id_0 < offset:
                 um_normal_0 = -n_0[0]  # TODO
-                cost[n_id_0] = djn(um_normal_0, 0., 0.)
+                cost[n_id_0] = djn(um_normal_0, 0.0, 0.0)
                 cost[n_id_0 + offset] = cost[n_id_0]
             if n_id_1 < offset:
                 um_normal_1 = -n_1[0]  # TODO
-                cost[n_id_1] = djn(um_normal_1, 0., 0.)
+                cost[n_id_1] = djn(um_normal_1, 0.0, 0.0)
                 cost[n_id_1 + offset] = cost[n_id_1]
         return cost
 
@@ -254,9 +268,7 @@ def subgradient(
         base_integrals,
         dt,
     ):
-        dj = contact_subgradient(
-            var, var_old, nodes, contact_boundary, contact_normals
-        )
+        dj = contact_subgradient(var, var_old, nodes, contact_boundary, contact_normals)
         result = np.dot(lhs, var) - rhs + dj
         result = result.ravel()
         return result

examples/Bagirrov_Bartman_Ochal_2024.py

@@ -1,6 +1,7 @@
 """
 Created at 21.08.2019
 """
+
 import pickle
 from dataclasses import dataclass
 from matplotlib import pyplot as plt
@@ -37,7 +38,7 @@ def normal_direction(u_nu: float) -> float:
         return k10 * (u_nu * 2) + k11
     if u_nu < 2 * mm:
         return k20 * (u_nu * 2) + k21
-    return u_nu ** 3 * 4 * k30
+    return u_nu**3 * 4 * k30
 
 
 class MMLV99(ContactLaw):
@@ -46,12 +47,12 @@ def potential_normal_direction(u_nu: float) -> float:
         if u_nu <= 0:
             return 0.0
         if u_nu < 0.5 * mm:
-            return k0 * u_nu ** 2
+            return k0 * u_nu**2
         if u_nu < 1 * mm:
-            return k10 * u_nu ** 2 + k11 * u_nu
+            return k10 * u_nu**2 + k11 * u_nu
         if u_nu < 2 * mm:
-            return k20 * u_nu ** 2 + k21 * u_nu + 4
-        return u_nu ** 4 * k30
+            return k20 * u_nu**2 + k21 * u_nu + 4
+        return u_nu**4 * k30
 
     @staticmethod
     def potential_tangential_direction(u_tau: np.ndarray) -> float:
@@ -63,7 +64,7 @@ def subderivative_normal_direction(u_nu: float, v_nu: float) -> float:
 
     @staticmethod
     def regularized_subderivative_tangential_direction(
-            u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7
+        u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7
     ) -> float:
         """
         Coulomb regularization
@@ -167,12 +168,11 @@ def outer_forces(x, t=None):
             drawer.draw(
                 show=config.show,
                 save=config.save,
-                title=f"{m}: {f}, "
+                title=f"{m}: {f}, ",
                 # f"time: {runner.step_solver.last_timing}"
             )
-            x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1,
-                0]
-            u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1]
+            x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0]
+            u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1]
             y1 = [normal_direction(-u_) for u_ in u]
             print(f)
             plt.plot(x, y1, label=f"{f:.2e}")
@@ -182,7 +182,7 @@ def outer_forces(x, t=None):
 
 
 if __name__ == "__main__":
-    X = np.linspace((2-2) * mm, (2+2) * mm, 1000)
+    X = np.linspace((2 - 2) * mm, (2 + 2) * mm, 1000)
     Y = np.empty(1000)
     for i in range(1000):
         Y[i] = MMLV99.potential_normal_direction(X[i])
@@ -236,6 +236,15 @@ def outer_forces(x, t=None):
     # # plt.loglog()
     # plt.show()
     methods = ("BFGS", "CG", "qsm", "Powell", "subgradient")
-    forces = (23e3 * kN, 25e3 * kN, 25.6e3 * kN, 25.9e3 * kN, 26e3 * kN,
-              26.1e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN)[-1::4]
+    forces = (
+        23e3 * kN,
+        25e3 * kN,
+        25.6e3 * kN,
+        25.9e3 * kN,
+        26e3 * kN,
+        26.1e3 * kN,
+        26.2e3 * kN,
+        27e3 * kN,
+        30e3 * kN,
+    )[-1::4]
     main(Config(save=False, show=True, force=False).init(), methods, forces)

examples/Makela_et_al_1998.py

@@ -1,6 +1,21 @@
-"""
-Created at 21.08.2019
-"""
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2023-2024  Piotr Bartman-Szwarc <[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 pickle
 from dataclasses import dataclass
 
@@ -30,17 +45,19 @@
 
 class MMLV99(PotentialOfContactLaw, DirectContactLaw):
     @staticmethod
-    def potential_normal_direction(var_nu: float, static_displacement_nu: float, dt: float) -> float:
+    def potential_normal_direction(
+        var_nu: float, static_displacement_nu: float, dt: float
+    ) -> float:
         u_nu = -var_nu
-        coef = 1.
+        coef = 1.0
         if u_nu <= 0:
             return 0.0
         if u_nu < 0.5 * mm:
-            return k0 * u_nu ** 2 * coef
+            return k0 * u_nu**2 * coef
         if u_nu < 1 * mm:
-            return (k10 * u_nu ** 2 + k11 * u_nu) * coef
+            return (k10 * u_nu**2 + k11 * u_nu) * coef
         if u_nu < 2 * mm:
-            return (k20 * u_nu ** 2 + k21 * u_nu + 4) * coef
+            return (k20 * u_nu**2 + k21 * u_nu + 4) * coef
         return 16 * coef
 
     @staticmethod
@@ -62,7 +79,7 @@ def subderivative_normal_direction(
     def potential_tangential_direction(
         var_tau: float, static_displacement_tau: float, dt: float
     ) -> float:
-        return np.log(np.sum(var_tau ** 2) ** 0.5 + 1)
+        return np.log(np.sum(var_tau**2) ** 0.5 + 1)
 
 
 @dataclass()
@@ -164,9 +181,8 @@ def outer_forces(x, t=None):
                 # title=f"{m}: {f}, "
                 # f"time: {runner.step_solver.last_timing}"
             )
-            x = state.body.mesh.nodes[:state.body.mesh.contact_nodes_count - 1,
-                0]
-            u = state.displacement[:state.body.mesh.contact_nodes_count - 1, 1]
+            x = state.body.mesh.nodes[: state.body.mesh.contact_nodes_count - 1, 0]
+            u = state.displacement[: state.body.mesh.contact_nodes_count - 1, 1]
             y1 = [MMLV99().normal_direction(-u_) for u_ in u]
             plt.plot(x, y1, label=f"{f:.2e}")
         plt.ylabel("Interlaminar binding force [kN/m$^2$]")
@@ -228,8 +244,19 @@ def outer_forces(x, t=None):
     #                     -0.8555741662642888,
     #                     -1.2663638426265278, ],
     # }
-    forces = np.asarray((20e3 * kN, 21e3 * kN, 21e3 * kN, 23e3 * kN,
-              25e3 * kN, 26e3 * kN, 26.2e3 * kN, 27e3 * kN, 30e3 * kN))[::2]
+    forces = np.asarray(
+        (
+            20e3 * kN,
+            21e3 * kN,
+            21e3 * kN,
+            23e3 * kN,
+            25e3 * kN,
+            26e3 * kN,
+            26.2e3 * kN,
+            27e3 * kN,
+            30e3 * kN,
+        )
+    )[::2]
     # # for m, losses in results.items():
     # #     plt.plot(forces/1e3, -1 * np.asarray(losses[:]), "-o", label=m)
     # plt.legend()