Versions bump

.github/workflows/conmech_tests.yml

@@ -14,8 +14,8 @@ jobs:
   run_tests:
     strategy:
       matrix:
-        platform: [ ubuntu-latest, macos-latest, windows-latest ]
-        python-version: [ "3.8", "3.9", "3.10", "3.11" ]
+        platform: [ ubuntu-latest, macos-13, windows-latest ]
+        python-version: [ "3.9", "3.10", "3.11"]
     runs-on: ${{ matrix.platform }}
     steps:
       - uses: actions/checkout@v2

conmech/dynamics/contact/__init__.py

@@ -0,0 +1,18 @@
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2024  Piotr Bartman-Szwarc <piotr.bartman@uj.edu.pl>
+#
+# 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.

conmech/dynamics/contact/constant_constact_law.py

@@ -0,0 +1,46 @@
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2024  Piotr Bartman-Szwarc <piotr.bartman@uj.edu.pl>
+#
+# 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 Type
+
+from conmech.dynamics.contact.contact_law import (
+    ContactLaw,
+    DirectContactLaw,
+    PotentialOfContactLaw,
+)
+
+
+def make_const_contact_law(resistance: float) -> Type[ContactLaw]:
+    class PSlopeContactLaw(DirectContactLaw, PotentialOfContactLaw):
+        @staticmethod
+        def potential_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+        ) -> float:
+            if var_nu <= 0:
+                return 0 * static_displacement_nu
+            return (resistance * var_nu) * static_displacement_nu
+
+        @staticmethod
+        def subderivative_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+        ) -> float:
+            if var_nu <= 0:
+                return 0
+            return resistance
+
+    return PSlopeContactLaw

conmech/dynamics/contact/contact_law.py

@@ -0,0 +1,126 @@
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2024  Piotr Bartman-Szwarc <piotr.bartman@uj.edu.pl>
+#
+# 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 abc import ABC, abstractmethod
+
+
+class ContactLaw:
+    # pylint: disable=unused-argument)
+    """
+    "Abstract" class for all contact conditions.
+    """
+
+    @staticmethod
+    def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float:
+        """
+        :param var_nu: variable normal vector length
+        :param static_displacement_nu: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns foundation response
+        """
+        return 1.0
+
+    @staticmethod
+    def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float:
+        """
+        Friction bound
+
+        :param var_nu: variable normal vector length
+        :param static_displacement_nu: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns foundation response
+        """
+        return 1.0
+
+
+class DirectContactLaw(ContactLaw, ABC):
+    # pylint: disable=unused-argument)
+    """
+    Abstract class for contact conditions given in a direct form.
+    Since usually contact law is a multifunction it can be treated as
+    a subderivative of *some* function - potential. Hence, from point of view
+    of conmech package we call subderivative form.
+    """
+
+    @staticmethod
+    @abstractmethod
+    def subderivative_normal_direction(
+        var_nu: float, static_displacement_nu: float, dt: float
+    ) -> float:
+        """
+        :param var_nu: variable normal vector length
+        :param static_displacement_nu: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns foundation response
+        """
+        raise NotImplementedError()
+
+    @staticmethod
+    def subderivative_tangential_direction(
+        var_tau: float, static_displacement_tau: float, dt: float
+    ) -> float:
+        """
+        :param var_tau: variable normal vector length
+        :param static_displacement_tau: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns potential of foundation friction response
+        """
+        return 0.0
+
+
+class PotentialOfContactLaw(ContactLaw, ABC):
+    # pylint: disable=unused-argument)
+    """
+    Abstract class for contact conditions given in a potential form.
+    Since usually contact law is a multifunction it can be treated as
+    a subderivative of *some* function - potential. To solve contact problem
+    numerically with optimization approach we need only the potential of
+    contact condition. Hence, from point of view of conmech package,
+    we call potential form.
+    """
+
+    @staticmethod
+    @abstractmethod
+    def potential_normal_direction(
+        var_nu: float, static_displacement_nu: float, dt: float
+    ) -> float:
+        """
+        :param var_nu: variable normal vector length
+        :param static_displacement_nu: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns potential of foundation response
+        """
+        raise NotImplementedError()
+
+    @staticmethod
+    def potential_tangential_direction(
+        var_tau: float, static_displacement_tau: float, dt: float
+    ) -> float:
+        """
+        :param var_tau: variable normal vector length
+        :param static_displacement_tau: normal vector length of displacement from
+         the previous step. For time independent problems this likely be 0.
+        :param dt: time step
+        :returns potential of foundation friction response
+        """
+        return 0

conmech/dynamics/contact/damped_normal_compliance.py

@@ -0,0 +1,55 @@
+# CONMECH @ Jagiellonian University in Kraków
+#
+# Copyright (C) 2024  Piotr Bartman-Szwarc <piotr.bartman@uj.edu.pl>
+#
+# 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 conmech.dynamics.contact.contact_law import PotentialOfContactLaw
+from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw
+
+
+def make_damped_norm_compl(obstacle_level: float, kappa: float, beta: float, interior=False):
+    superclass = InteriorContactLaw if interior else PotentialOfContactLaw
+
+    class DampedNormalCompliance(superclass):
+        @property
+        def kappa(self):
+            return kappa
+
+        @property
+        def beta(self):
+            return beta
+
+        @property
+        def obstacle_level(self):
+            return obstacle_level
+
+        @staticmethod
+        def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float:
+            """
+            Since multiply by var_nu
+            """
+            return 0.5 * var_nu
+
+        @staticmethod
+        def potential_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+        ) -> float:
+            displacement = static_displacement_nu + var_nu * dt
+            if displacement < obstacle_level:
+                return 0
+            return kappa * (displacement - obstacle_level) + beta * var_nu
+
+    return DampedNormalCompliance

conmech/mesh/utils.py → conmech/dynamics/contact/interior_contact_law.py

@@ -1,6 +1,6 @@
 # CONMECH @ Jagiellonian University in Kraków
 #
-# Copyright (C) 2023  Piotr Bartman <piotr.bartman@uj.edu.pl>
+# Copyright (C) 2024  Piotr Bartman-Szwarc <piotr.bartman@uj.edu.pl>
 #
 # This program is free software; you can redistribute it and/or
 # modify it under the terms of the GNU General Public License
@@ -16,16 +16,8 @@
 # along with this program; if not, write to the Free Software
 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,
 # USA.
-import numpy as np
+from conmech.dynamics.contact.contact_law import PotentialOfContactLaw
 
 
-def interpolate_nodes(scaled_nodes, corner_vectors):
-    input_dim = scaled_nodes.shape[-1]
-    output_dim = corner_vectors.shape[-1]
-    values = np.zeros((scaled_nodes.shape[0], output_dim))
-    for i in range(input_dim):
-        coordinate_i = scaled_nodes[..., [i]]
-        values += (
-            coordinate_i * corner_vectors[i] + (1 - coordinate_i) * corner_vectors[i + input_dim]
-        ) / input_dim
-    return values
+class InteriorContactLaw(PotentialOfContactLaw):
+    pass

conmech/dynamics/contact/relu_slope_contact_law.py

@@ -0,0 +1,32 @@
+"""
+Example contact law
+"""
+
+from typing import Type
+
+from conmech.dynamics.contact.contact_law import (
+    ContactLaw,
+    DirectContactLaw,
+    PotentialOfContactLaw,
+)
+
+
+def make_slope_contact_law(slope: float) -> Type[ContactLaw]:
+    class ReLuContactLaw(DirectContactLaw, PotentialOfContactLaw):
+        @staticmethod
+        def potential_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+        ) -> float:
+            if var_nu <= 0:
+                return 0.0
+            return 0.5 * slope * var_nu**2
+
+        @staticmethod
+        def subderivative_normal_direction(
+            var_nu: float, static_displacement_nu: float, dt: float
+        ) -> float:
+            if var_nu <= 0:
+                return 0.0
+            return slope * var_nu * static_displacement_nu
+
+    return ReLuContactLaw

conmech/dynamics/dynamics.py

@@ -23,9 +23,8 @@ def __init__(
         self.body.dynamics = self
 
         self.force = BodyForces(body)
-        self.temperature = BodyForces(
-            body,
-        )
+        self.temperature = BodyForces(body)
+
         self.factory = get_factory(body.mesh.dimension)
         self.element_initial_volume: np.ndarray
         self.volume_at_nodes: np.ndarray

conmech/dynamics/factory/_abstract_dynamics_factory.py

@@ -1,6 +1,8 @@
 from typing import Tuple
 import numpy as np
 
+from conmech.struct.stiffness_matrix import SM1
+
 
 class AbstractDynamicsFactory:
     @property
@@ -34,5 +36,9 @@ def get_permittivity_tensor(self, W: np.ndarray, coeff: np.ndarray) -> np.ndarra
         raise NotImplementedError()
 
     @staticmethod
-    def calculate_poisson_matrix(W: np.ndarray) -> np.ndarray:
+    def calculate_poisson_matrix(W: np.ndarray, propagation: float) -> SM1:
+        return SM1(propagation**2 * np.sum(W.diagonal(), axis=2))
+
+    @staticmethod
+    def calculate_wave_matrix(W: np.ndarray) -> np.ndarray:
         return np.sum(W.diagonal(), axis=2)