BBO: DC functions

conmech/solvers/optimization/optimization.py

@@ -84,6 +84,13 @@ def __init__(
             )
         else:
             self.subgradient = None
+        if hasattr(contact_law, "sub2derivative_normal_direction"):  # TODO
+            self.sub2gradient = make_subgradient(
+                normal_condition=contact_law.sub2derivative_normal_direction,
+                only_boundary=True,
+            )
+        else:
+            self.sub2gradient = None
         if isinstance(statement, WaveStatement):
             if isinstance(contact_law, InteriorContactLaw):
                 self.loss = make_equation(  # TODO!
@@ -150,7 +157,7 @@ def _solve_impl(
             # pylint: disable=import-outside-toplevel,import-error)
             from kosopt.qsmlm import make_minimizer
 
-            self.minimizer = make_minimizer(self.loss, self.subgradient)
+            self.minimizer = make_minimizer(self.loss, self.subgradient, self.sub2gradient)
 
         while norm >= fixed_point_abs_tol:
             if method.lower() in QSMLM_NAMES:

conmech/solvers/solver_methods.py

@@ -228,6 +228,7 @@ def make_subgradient(
     tangential_condition_bound: Optional[Callable] = None,
     problem_dimension=2,
     variable_dimension=2,
+    only_boundary=False,
 ):
     normal_condition = njit(normal_condition)
     normal_condition_bound = njit(normal_condition_bound, value=1)
@@ -304,7 +305,10 @@ def cost_functional(
             var, var_old, u_vector, nodes, contact_boundary, contact_normals, dt
         )
         ind = lhs.shape[0]
-        result_ = np.dot(lhs, var[:ind]) - rhs + dj
+        if only_boundary:
+            result_ = dj
+        else:
+            result_ = np.dot(lhs, var[:ind]) - rhs + dj
         result[:ind] = result_.ravel()
         return result
 

examples/Bagirrov_Bartman_Ochal_2024.py

@@ -65,6 +65,27 @@ def subderivative_normal_direction(
             return k20 * (var_nu * 2) + k21
         return var_nu**3 * 4 * k30
 
+    @staticmethod
+    def sub2derivative_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+    ) -> float:
+        acc = 0.0
+        p1 = 0.5 * mm
+        if var_nu <= p1:
+            return acc
+
+        acc += (k0 * p1 * 2) - (k10 * (p1 * 2) + k11)
+        p2 = 1 * mm
+        if var_nu < p2:
+            return acc
+
+        acc += (k10 * (p2 * 2) + k11) - (k20 * (p2 * 2) + k21)
+        p3 = 2 * mm
+        if var_nu < p3:
+            return acc
+
+        return acc
+
 
 @dataclass()
 class StaticSetup(StaticDisplacementProblem):
@@ -200,6 +221,10 @@ def outer_forces(x, t=None):
         Y[i] = MMLV99().subderivative_normal_direction(X[i], 0.0, 0.0)
     plt.plot(X, Y)
     plt.show()
+    for i in range(1000):
+        Y[i] = MMLV99().sub2derivative_normal_direction(X[i], 0.0, 0.0)
+    plt.plot(X, Y)
+    plt.show()
 
     methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "qsm", "Powell", "globqsm")[:]
     forces = (

examples/Makela_et_al_1998.py

@@ -76,6 +76,30 @@ def subderivative_normal_direction(
             return k20 * (u_nu * 2) + k21
         return 0.0
 
+    @staticmethod
+    def sub2derivative_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+    ) -> float:
+        u_nu = -var_nu
+
+        acc = 0.0
+        p1 = 0.5 * mm
+        if u_nu <= p1:
+            return acc
+
+        acc += (k0 * p1 * 2) - (k10 * (p1 * 2) + k11)
+        p2 = 1.0 * mm
+        if u_nu < p2:
+            return acc
+
+        acc += (k10 * (p2 * 2) + k11) - (k20 * (p2 * 2) + k21)
+        p3 = 2.0 * mm
+        if u_nu < p3:
+            return acc
+
+        acc += (k20 * (p3 * 2) + k21) - 0.0
+        return acc
+
     @staticmethod
     def potential_tangential_direction(
         var_tau: float, static_displacement_tau: float, dt: float
@@ -188,8 +212,8 @@ def outer_forces(x, t=None):
     print("Plotting...")\
 
     path = f"{config.outputs_path}/{PREFIX}_losses"
-    if config.show:
-        plot_losses(path)
+    # if config.show:
+    plot_losses(path)
 
     for m in methods:
         for f in forces:
@@ -265,14 +289,22 @@ def loss_value(state, runner) -> float:
 
 
 if __name__ == "__main__":
-    # X = np.linspace(0, -3 * mm, 1000)
-    # Y = np.empty(1000)
+    X = np.linspace(0, -3 * mm, 1000)
+    Y = np.empty(1000)
+    for i in range(1000):
+        Y[i] = MMLV99.potential_normal_direction(X[i], 0, 0)
+    plt.plot(X, Y)
+    plt.show()
+    for i in range(1000):
+        Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0)
+    plt.plot(X, Y)
+    plt.show()
     # for i in range(1000):
-    #     Y[i] = MMLV99.potential_normal_direction(X[i])
+    #     Y[i] = MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
     # plt.plot(X, Y)
     # plt.show()
     # for i in range(1000):
-    #     Y[i] = MMLV99.normal_direction(X[i])
+    #     Y[i] = MMLV99.subderivative_normal_direction(X[i], 0, 0) + MMLV99.sub2derivative_normal_direction(X[i], 0, 0)
     # plt.plot(X, Y)
     # plt.show()
 
@@ -290,5 +322,5 @@ def loss_value(state, runner) -> float:
         )
     )[:]
 
-    methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[:]
+    methods = ("gradiented BFGS", "gradiented CG", "BFGS", "CG", "Powell", "globqsm", "qsm")[-1:]
     main(Config(save=False, show=False, force=True).init(), methods, forces)