diff --git a/.github/workflows/conmech_tests.yml b/.github/workflows/conmech_tests.yml index ece6637a..8a9bb02a 100644 --- a/.github/workflows/conmech_tests.yml +++ b/.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 diff --git a/conmech/dynamics/contact/__init__.py b/conmech/dynamics/contact/__init__.py new file mode 100644 index 00000000..8961ee91 --- /dev/null +++ b/conmech/dynamics/contact/__init__.py @@ -0,0 +1,18 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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. diff --git a/conmech/dynamics/contact/constant_constact_law.py b/conmech/dynamics/contact/constant_constact_law.py new file mode 100644 index 00000000..2dc56170 --- /dev/null +++ b/conmech/dynamics/contact/constant_constact_law.py @@ -0,0 +1,46 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 diff --git a/conmech/dynamics/contact/contact_law.py b/conmech/dynamics/contact/contact_law.py new file mode 100644 index 00000000..3202e96d --- /dev/null +++ b/conmech/dynamics/contact/contact_law.py @@ -0,0 +1,126 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 diff --git a/conmech/dynamics/contact/damped_normal_compliance.py b/conmech/dynamics/contact/damped_normal_compliance.py new file mode 100644 index 00000000..de25c502 --- /dev/null +++ b/conmech/dynamics/contact/damped_normal_compliance.py @@ -0,0 +1,55 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 diff --git a/conmech/mesh/utils.py b/conmech/dynamics/contact/interior_contact_law.py similarity index 59% rename from conmech/mesh/utils.py rename to conmech/dynamics/contact/interior_contact_law.py index 7e4fe316..dba39fce 100644 --- a/conmech/mesh/utils.py +++ b/conmech/dynamics/contact/interior_contact_law.py @@ -1,6 +1,6 @@ # CONMECH @ Jagiellonian University in Kraków # -# Copyright (C) 2023 Piotr Bartman +# Copyright (C) 2024 Piotr Bartman-Szwarc # # 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 diff --git a/conmech/dynamics/contact/relu_slope_contact_law.py b/conmech/dynamics/contact/relu_slope_contact_law.py new file mode 100644 index 00000000..7e2e8513 --- /dev/null +++ b/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 diff --git a/conmech/dynamics/dynamics.py b/conmech/dynamics/dynamics.py index 9dcd96b2..13933a3b 100644 --- a/conmech/dynamics/dynamics.py +++ b/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 diff --git a/conmech/dynamics/factory/_abstract_dynamics_factory.py b/conmech/dynamics/factory/_abstract_dynamics_factory.py index 69534853..c9e61e65 100644 --- a/conmech/dynamics/factory/_abstract_dynamics_factory.py +++ b/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) diff --git a/conmech/dynamics/factory/_dynamics_factory_2d.py b/conmech/dynamics/factory/_dynamics_factory_2d.py index 9c4172d6..b1767439 100644 --- a/conmech/dynamics/factory/_dynamics_factory_2d.py +++ b/conmech/dynamics/factory/_dynamics_factory_2d.py @@ -3,6 +3,7 @@ import numpy as np from conmech.dynamics.factory._abstract_dynamics_factory import AbstractDynamicsFactory +from conmech.struct.stiffness_matrix import SM2, SM1, SM1to2 DIMENSION = 2 ELEMENT_NODES_COUNT = 3 @@ -159,33 +160,31 @@ def get_edges_features_matrix(self, elements, nodes): def dimension(self) -> int: return DIMENSION - def calculate_constitutive_matrices( - self, W: np.ndarray, mu: float, lambda_: float - ) -> np.ndarray: + def calculate_constitutive_matrices(self, W: np.ndarray, mu: float, lambda_: float) -> SM2: A_11 = (2 * mu + lambda_) * W[0, 0] + mu * W[1, 1] A_12 = mu * W[1, 0] + lambda_ * W[0, 1] A_21 = lambda_ * W[1, 0] + mu * W[0, 1] A_22 = mu * W[0, 0] + (2 * mu + lambda_) * W[1, 1] - return np.block([[A_11, A_12], [A_21, A_22]]) + return SM2(np.block([[A_11, A_12], [A_21, A_22]])) def get_relaxation_tensor(self, W, coeff): A_11 = coeff[0][0][0] * W[0, 0] + coeff[0][1][1] * W[1, 1] A_12 = coeff[0][1][0] * W[1, 0] + coeff[0][0][1] * W[0, 1] A_21 = coeff[1][1][0] * W[1, 0] + coeff[1][0][1] * W[0, 1] A_22 = coeff[1][0][0] * W[0, 0] + coeff[1][1][1] * W[1, 1] - return np.block([[A_11, A_12], [A_21, A_22]]) + return SM2(np.block([[A_11, A_12], [A_21, A_22]])) def calculate_acceleration(self, U, density): Z = np.zeros_like(U) - return density * np.block([[U, Z], [Z, U]]) + return SM2(density * np.block([[U, Z], [Z, U]])) def calculate_thermal_expansion(self, V, coeff): A_11 = coeff[0][0] * V[0] + coeff[0][1] * V[1] A_22 = coeff[1][0] * V[0] + coeff[1][1] * V[1] - return np.block([A_11, A_22]) + return SM1to2(np.block([A_11, A_22])) def calculate_thermal_conductivity(self, W, coeff): - return ( + return SM1( coeff[0][0] * W[0, 0] + coeff[0][1] * W[0, 1] + coeff[1][0] * W[1, 0] @@ -197,10 +196,10 @@ def get_piezoelectric_tensor(self, W, coeff): A_12 = coeff[0][1][0] * W[1, 0] + coeff[0][0][1] * W[0, 1] A_21 = coeff[1][1][0] * W[1, 0] + coeff[1][0][1] * W[0, 1] A_22 = coeff[1][0][0] * W[0, 0] + coeff[1][1][1] * W[1, 1] - return np.block([[A_11 + A_12, A_22 + A_21]]) + return SM2(np.block([[A_11 + A_12, A_22 + A_21]])) def get_permittivity_tensor(self, W, coeff): - return ( + return SM1( coeff[0][0] * W[0, 0] + coeff[0][1] * W[0, 1] + coeff[1][0] * W[1, 0] diff --git a/conmech/dynamics/factory/_dynamics_factory_3d.py b/conmech/dynamics/factory/_dynamics_factory_3d.py index c6ef16c9..65411b32 100644 --- a/conmech/dynamics/factory/_dynamics_factory_3d.py +++ b/conmech/dynamics/factory/_dynamics_factory_3d.py @@ -2,6 +2,7 @@ import numpy as np from conmech.dynamics.factory._abstract_dynamics_factory import AbstractDynamicsFactory +from conmech.struct.stiffness_matrix import SM3, SM1, SM1to3 DIMENSION = 3 ELEMENT_NODES_COUNT = 4 @@ -163,20 +164,28 @@ def calculate_constitutive_matrices(self, W, mu, lambda_): A_13 = lambda_ * W[2, 0] + mu * W[0, 2] A_23 = lambda_ * W[2, 1] + mu * W[1, 2] - return np.block([[A_11, A_12, A_13], [A_21, A_22, A_23], [A_31, A_32, A_33]]) + return SM3( + np.block( + [ + [A_11, A_12, A_13], + [A_21, A_22, A_23], + [A_31, A_32, A_33], + ] + ) + ) def calculate_acceleration(self, U, density): Z = np.zeros_like(U) - return density * np.block([[U, Z, Z], [Z, U, Z], [Z, Z, U]]) + return SM3(density * np.block([[U, Z, Z], [Z, U, Z], [Z, Z, U]])) def calculate_thermal_expansion(self, V, coeff): A_11 = coeff[0][0] * V[0] + coeff[0][1] * V[1] + coeff[0][2] * V[2] A_22 = coeff[1][0] * V[0] + coeff[1][1] * V[1] + coeff[1][2] * V[2] A_33 = coeff[2][0] * V[0] + coeff[2][1] * V[1] + coeff[2][2] * V[2] - return np.block([A_11, A_22, A_33]) + return SM1to3(np.block([A_11, A_22, A_33])) def calculate_thermal_conductivity(self, W, coeff): - return ( + return SM1( coeff[0][0] * W[0, 0] + coeff[0][1] * W[0, 1] + coeff[0][2] * W[0, 2] diff --git a/conmech/dynamics/factory/dynamics_factory_method.py b/conmech/dynamics/factory/dynamics_factory_method.py index 3d4ffa6f..6247d642 100644 --- a/conmech/dynamics/factory/dynamics_factory_method.py +++ b/conmech/dynamics/factory/dynamics_factory_method.py @@ -7,6 +7,7 @@ TemperatureBodyProperties, PiezoelectricBodyProperties, BodyProperties, + MembraneProperties, ) @@ -51,21 +52,22 @@ def get_dynamics(elements: np.ndarray, body_prop: BodyProperties, U, V, W): dimension = len(elements[0]) - 1 factory = get_factory(dimension) - poisson_operator = factory.calculate_poisson_matrix(W) + acceleration_operator = factory.calculate_acceleration(U, body_prop.mass_density) - elasticity = ( - factory.calculate_constitutive_matrices(W, body_prop.mu, body_prop.lambda_) - if isinstance(body_prop, ElasticProperties) - else None - ) + if isinstance(body_prop, MembraneProperties): + poisson_operator = factory.calculate_poisson_matrix(W, body_prop.propagation) + else: + poisson_operator = None - viscosity = ( - factory.calculate_constitutive_matrices(W, body_prop.theta, body_prop.zeta) - if isinstance(body_prop, ViscoelasticProperties) - else None - ) + if isinstance(body_prop, ElasticProperties): + elasticity = factory.calculate_constitutive_matrices(W, body_prop.mu, body_prop.lambda_) + else: + elasticity = None - acceleration_operator = factory.calculate_acceleration(U, body_prop.mass_density) + if isinstance(body_prop, ViscoelasticProperties): + viscosity = factory.calculate_constitutive_matrices(W, body_prop.theta, body_prop.zeta) + else: + viscosity = None if isinstance(body_prop, TemperatureBodyProperties): thermal_expansion = factory.calculate_thermal_expansion(V, body_prop.thermal_expansion) diff --git a/conmech/dynamics/statement.py b/conmech/dynamics/statement.py index d2dc8e4b..78a11eaf 100644 --- a/conmech/dynamics/statement.py +++ b/conmech/dynamics/statement.py @@ -3,6 +3,8 @@ import numpy as np +from conmech.struct.stiffness_matrix import SM1 + @dataclass class Variables: @@ -16,9 +18,10 @@ class Variables: class Statement: - def __init__(self, body, dimension): + def __init__(self, body, dimensions: tuple): self.body = body - self.dimension = dimension + self.dimension_in = dimensions[0] + self.dimension_out = dimensions[1] self.left_hand_side = None self.right_hand_side = None self.dirichlet_cond_name = "dirichlet" @@ -39,13 +42,13 @@ def apply_dirichlet_condition(self): c = self.body.mesh.boundaries.boundaries[dirichlet_cond].node_condition node_count = self.body.mesh.nodes_count for i, j in self.body.mesh.boundaries.get_all_boundary_indices( - dirichlet_cond, node_count, self.dimension + dirichlet_cond, node_count, self.dimension_in ): self.right_hand_side[:] -= self.left_hand_side[:, i] @ c[j] - self.left_hand_side[:, i] = 0 - self.left_hand_side[i, :] = 0 + self.left_hand_side.data[:, i] = 0 + self.left_hand_side.data[i, :] = 0 # have to be "[i][:, i]" instead of just a "[i, i]" because the i may be ndarray - self.left_hand_side[i][:, i] = np.eye(j.stop - j.start) + self.left_hand_side.data[i][:, i] = np.eye(j.stop - j.start) self.right_hand_side[i] = c[j] def find_dirichlet_conditions(self): @@ -65,7 +68,7 @@ def find_dirichlet_conditions(self): class StaticPoissonStatement(Statement): def __init__(self, dynamics): - super().__init__(dynamics, 1) + super().__init__(dynamics, (1, 1)) def update_left_hand_side(self, var: Variables): self.left_hand_side = self.body.dynamics.poisson_operator.copy() @@ -74,9 +77,36 @@ def update_right_hand_side(self, var: Variables): self.right_hand_side = self.body.dynamics.temperature.integrate(time=0) +class WaveStatement(Statement): + def __init__(self, body): + super().__init__(body, (1, 1)) + + def update_left_hand_side(self, var): + assert var.time_step is not None + + self.left_hand_side = ( + (1 / var.time_step) * self.body.dynamics.acceleration_operator.SM1 + + self.body.dynamics.poisson_operator * var.time_step + ) + + def update_right_hand_side(self, var): + assert var.displacement is not None + assert var.velocity is not None + assert var.time_step is not None + assert var.time is not None + + ind = self.body.mesh.nodes_count # 1 dimensional + + A = -1 * self.body.dynamics.poisson_operator @ var.displacement[:ind] + + A += (1 / var.time_step) * self.body.dynamics.acceleration_operator.SM1 @ var.velocity[:ind] + + self.right_hand_side = self.body.dynamics.force.integrate(time=var.time) + A + + class StaticDisplacementStatement(Statement): def __init__(self, body): - super().__init__(body, body.mesh.dimension) + super().__init__(body, (body.mesh.dimension, body.mesh.dimension)) def update_left_hand_side(self, var: Variables): self.left_hand_side = self.body.dynamics.elasticity.copy() @@ -87,7 +117,7 @@ def update_right_hand_side(self, var: Variables): class QuasistaticRelaxationStatement(Statement): def __init__(self, body): - super().__init__(body, 2) + super().__init__(body, (2, 2)) def update_left_hand_side(self, var: Variables): assert var.time_step is not None @@ -109,13 +139,13 @@ def update_right_hand_side(self, var: Variables): class QuasistaticVelocityStatement(Statement): def __init__(self, body): - super().__init__(body, 2) + super().__init__(body, (2, 2)) def update_left_hand_side(self, var: Variables): assert var.time_step is not None self.left_hand_side = ( - self.body.dynamics.viscosity.copy() + self.body.dynamics.elasticity * var.time_step + self.body.dynamics.viscosity + self.body.dynamics.elasticity * var.time_step ) def update_right_hand_side(self, var: Variables): @@ -130,7 +160,7 @@ def update_right_hand_side(self, var: Variables): class DynamicVelocityStatement(Statement): def __init__(self, body): - super().__init__(body, 2) + super().__init__(body, (2, 2)) def update_left_hand_side(self, var): assert var.time_step is not None @@ -167,16 +197,17 @@ def update_right_hand_side(self, var): class TemperatureStatement(Statement): def __init__(self, body): - super().__init__(body, 1) + super().__init__(body, (1, 1)) def update_left_hand_side(self, var): assert var.time_step is not None ind = self.body.mesh.nodes_count # 1 dimensional - self.left_hand_side = (1 / var.time_step) * self.body.dynamics.acceleration_operator[ + left_hand_side = (1 / var.time_step) * self.body.dynamics.acceleration_operator[ :ind, :ind ] + self.body.dynamics.thermal_conductivity[:ind, :ind] + self.left_hand_side = SM1(left_hand_side) def update_right_hand_side(self, var): assert var.velocity is not None @@ -197,7 +228,7 @@ def update_right_hand_side(self, var): class PiezoelectricStatement(Statement): def __init__(self, body): - super().__init__(body, 1) + super().__init__(body, (1, 1)) self.dirichlet_cond_name = "piezo_" + self.dirichlet_cond_name def update_left_hand_side(self, var): diff --git a/conmech/helpers/config.py b/conmech/helpers/config.py index 203c1dce..aed923be 100644 --- a/conmech/helpers/config.py +++ b/conmech/helpers/config.py @@ -27,10 +27,10 @@ def init(self) -> "Config": """ if self.show and self.save: raise ValueError("Cannot show and save at once!") + self.path.mkdir(parents=True, exist_ok=True) + return self + @property + def path(self): output_dir = self.output_dir or str(self.timestamp) - path_part_1 = pathlib.Path(self.outputs_path) - path_part_2 = pathlib.Path(output_dir) - path = path_part_1 / path_part_2 - path.mkdir(parents=True, exist_ok=True) - return self + return pathlib.Path(self.outputs_path) / pathlib.Path(output_dir) diff --git a/conmech/helpers/nph.py b/conmech/helpers/nph.py index 6c367586..18741eff 100644 --- a/conmech/helpers/nph.py +++ b/conmech/helpers/nph.py @@ -3,105 +3,31 @@ """ from ctypes import ArgumentError -from typing import Optional, Tuple import numba import numpy as np -def stack(data: np.ndarray) -> np.ndarray: - return data.T.flatten() - - def stack_column(data: np.ndarray) -> np.ndarray: return data.T.flatten().reshape(-1, 1) -stack_column_numba = numba.njit(stack_column) - - def unstack(vector: np.ndarray, dim: int) -> np.ndarray: return vector.reshape(-1, dim, order="F") -def unstack_and_sum_columns(data: np.ndarray, dim: int, keepdims: bool = False) -> np.ndarray: - return np.sum(unstack(data, dim), axis=1, keepdims=keepdims) - - def elementwise_dot( matrix_1: np.ndarray, matrix_2: np.ndarray, keepdims: bool = False ) -> np.ndarray: return (matrix_1 * matrix_2).sum(axis=1, keepdims=keepdims) -def get_occurances(data: np.ndarray) -> np.ndarray: - return np.array(list(set(data.flatten()))) - - -def close_modulo(value: np.ndarray, divider: Optional[int]) -> bool: - if divider is None: - return True - return np.allclose(value % divider, 0.0) or np.allclose(value % divider, divider) - - -def euclidean_norm(vector: np.ndarray, keepdims=False) -> np.ndarray: - data = (vector**2).sum(axis=-1, keepdims=keepdims) - if isinstance(vector, np.ndarray): - return np.sqrt(data) - return data.sqrt() - # return np.linalg.norm(vector, axis=-1) - # return np.sqrt(np.sum(vector ** 2, axis=-1))[..., np.newaxis] - - @numba.njit def euclidean_norm_numba(vector: np.ndarray) -> np.ndarray: data = (vector**2).sum(axis=-1) return np.sqrt(data) -def get_normal(vector: np.ndarray, normal: np.ndarray) -> np.ndarray: - return elementwise_dot(vector, normal, keepdims=True) - - -def get_normal_tangential(vector: np.ndarray, normal: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: - normal_vector = get_normal(vector, normal) - tangential_vector = vector - (normal_vector * normal) - return normal_vector, tangential_vector - - -def get_tangential(vector: np.ndarray, normal: np.ndarray) -> np.ndarray: - _, tangential_vector = get_normal_tangential(vector, normal) - return tangential_vector - - -@numba.njit -def get_tangential_numba(vector: np.ndarray, normal: np.ndarray) -> np.ndarray: - normal_vector = vector @ normal - tangential_vector = vector - (normal_vector * normal) - return tangential_vector - - -def orthogonalize_gram_schmidt(vectors: np.ndarray) -> np.ndarray: - # Gramm-schmidt orthog. - b0 = vectors[0] - if len(vectors) == 1: - return np.array((b0)) - - b1 = vectors[1] - (vectors[1] @ b0) * b0 - if len(vectors) == 2: - return np.array((b0, b1)) - - # MGS for stability - w2 = vectors[2] - (vectors[2] @ b0) * b0 - b2 = w2 - (w2 @ b1) * b1 - # nx = np.cross(ny,nz) - return np.array((b0, b1, b2)) - - -def get_in_base(vectors: np.ndarray, base: np.ndarray) -> np.ndarray: - return vectors @ base.T - - @numba.njit def get_node_index_numba(node, nodes): for i, n in enumerate(nodes): @@ -110,10 +36,6 @@ def get_node_index_numba(node, nodes): raise ArgumentError -def generate_normal(rows: int, columns: int, scale: float) -> np.ndarray: - return np.random.normal(loc=0.0, scale=scale * 0.5, size=[rows, columns]) - - @numba.njit(inline="always") def length(edge, nodes): return np.sqrt( diff --git a/conmech/mesh/boundaries_factory.py b/conmech/mesh/boundaries_factory.py index e1a20a04..811c5b64 100644 --- a/conmech/mesh/boundaries_factory.py +++ b/conmech/mesh/boundaries_factory.py @@ -54,10 +54,25 @@ def apply_predicate_to_surfaces(surfaces, nodes, predicate: Callable): def reorder_boundary_nodes(nodes, elements, is_contact, is_dirichlet): + """ + Node order + + [C C C N N N N N I I I I I I I D D D] + C - contact node + N - neumann node + I - internal node + D - dirichlet node + Dirichlet override other + """ # move boundary nodes to the top nodes, elements, boundary_nodes_count = reorder(nodes, elements, lambda _: True, to_top=True) # then move contact nodes to the top - nodes, elements, contact_nodes_count = reorder(nodes, elements, is_contact, to_top=True) + nodes, elements, contact_nodes_count = reorder( + nodes, + elements, + lambda *args: is_contact(*args) and not is_dirichlet(*args), + to_top=True, + ) # finally move dirichlet nodes to the bottom nodes, elements, dirichlet_nodes_count = reorder(nodes, elements, is_dirichlet, to_top=False) return ( diff --git a/conmech/mesh/interpolators.py b/conmech/mesh/interpolators.py deleted file mode 100644 index db108957..00000000 --- a/conmech/mesh/interpolators.py +++ /dev/null @@ -1,76 +0,0 @@ -# CONMECH @ Jagiellonian University in Kraków -# -# 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 random - -import numpy as np - -from conmech.helpers import nph -from conmech.mesh.utils import interpolate_nodes - - -def decide(scale): - return np.random.uniform(low=0, high=1) < scale - - -def choose(options): - return random.choice(options) - - -def get_mean(dimension, scale): - return np.random.uniform( - low=-scale, - high=scale, - size=(1, dimension), - ) - - -def get_corner_vectors_rotate(dimension, scale): - if dimension != 2: - raise NotImplementedError - # 1 2 - # 0 3 - corner_vector = nph.generate_normal(rows=1, columns=dimension, scale=scale) - corner_vectors = corner_vector * [[1, 1], [-1, 1], [-1, -1], [1, -1]] - return corner_vectors - - -def get_corner_vectors_all(dimension, scale): - corner_vectors = nph.generate_normal(rows=dimension * 2, columns=dimension, scale=scale) - return corner_vectors - - -def scale_nodes_to_square(nodes): - nodes_min = np.min(nodes, axis=0) - nodes_max = np.max(nodes, axis=0) - scaled_nodes = (nodes - nodes_min) / (nodes_max - nodes_min) - return scaled_nodes - - -def get_nodes_interpolation(nodes: np.ndarray, base: np.ndarray, corner_vectors: np.ndarray): - # orthonormal matrix; inverse equals transposition - denormalized_nodes = nph.get_in_base(nodes, base.T) - denormalized_scaled_nodes = denormalized_nodes # scale_nodes_to_square(denormalized_nodes) - - denormalized_nodes_interpolation = interpolate_nodes( - scaled_nodes=denormalized_scaled_nodes, - corner_vectors=corner_vectors, - ) - nodes_interpolation = ( - denormalized_nodes_interpolation # nph.get_in_base(denormalized_nodes_interpolation, base) - ) - return nodes_interpolation diff --git a/conmech/mesh/zoo/rectangle.py b/conmech/mesh/zoo/rectangle.py index 466951cb..e4a62258 100644 --- a/conmech/mesh/zoo/rectangle.py +++ b/conmech/mesh/zoo/rectangle.py @@ -32,8 +32,9 @@ def __init__(self, mesh_descr: RectangleMeshDescription): # pylint: disable=no-member nodes, elements = meshzoo.rectangle_tri( - np.linspace(0.0, scale_x, int(mesh_density[0]) + 1), - np.linspace(0.0, scale_y, int(mesh_density[1]) + 1), + (0.0, 0.0), + (scale_x, scale_y), + n=(int(mesh_density[0]) + 1, int(mesh_density[1]) + 1), variant="zigzag", ) super().__init__(nodes, elements) diff --git a/conmech/plotting/__init__.py b/conmech/plotting/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/conmech/plotting/membrane.py b/conmech/plotting/membrane.py new file mode 100644 index 00000000..dd49f8ec --- /dev/null +++ b/conmech/plotting/membrane.py @@ -0,0 +1,218 @@ +from typing import Callable, List + +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.cm import ScalarMappable, get_cmap +from matplotlib import tri + +from conmech.state.products.intersection_contact_limit_points import ( + VerticalIntersectionContactLimitPoints, +) +from conmech.state.state import State + + +def plot_in_columns(states: List[State], *args, **kwargs): + fig = plt.figure(layout="compressed", figsize=(8.3, 11.7)) + fig.suptitle(kwargs["title"]) + del kwargs["title"] + + in3d = kwargs.get("in3d", False) + if in3d: + extra_kwargs = {"projection": "3d"} + else: + extra_kwargs = {} + kwargs["in3d"] = False + + cols = 2 + rows = len(states) // 2 + i = 1 + for r in range(1, len(states) + 1): + ax = fig.add_subplot(rows, cols, i, **extra_kwargs) + do_plot(fig, states[r - 1], *args, ax=ax, elev=15, azim=30, **kwargs) + i += 1 + # ax = fig.add_subplot(rows, cols, i, projection='3d') + # do_plot(fig, states[r-1], *args, ax=ax, elev=90, azim=90, **kwargs) + # i += 1 + if kwargs["finish"]: + plt.show() + + +def plot_limit_points( + prod: VerticalIntersectionContactLimitPoints, + color="black", + title=None, + label=None, + finish=True, + ylim=(0, 1), +): + buff = np.zeros((2, 1000)) + buff_size = 0 + for time, zeros in prod.data.items(): + for zero in zeros: # significant performance improvement + buff[0, buff_size] = time + buff[1, buff_size] = zero + buff_size += 1 + if buff_size == 1000: + plt.scatter(buff[0, :], buff[1, :], s=1, color=color, label=label) + label = None + buff_size = 0 + if buff_size > 0: + plt.scatter(buff[0, :buff_size], buff[1, :buff_size], s=1, color=color, label=label) + plt.title(title) + plt.ylim(*ylim) + + if finish: + plt.show() + + +def plot(state: State, *args, **kwargs): + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1, projection="3d") + do_plot(state, *args, ax=ax, **kwargs) + plt.show() + + +# pylint: disable=too-many-arguments, too-many-locals +def do_plot( + fig, + state: State, + field="displacement", + vmin=None, + vmax=None, + zmin=None, + zmax=None, + title=None, + ax=None, + elev=0, + azim=-0, + in3d=False, + x=0.0, + **_kwargs, +): + assert state.body.mesh.dimension == 2 # membrane have to be 2D + mesh_node_x = state.body.mesh.nodes[:, 0] + mesh_node_y = state.body.mesh.nodes[:, 1] + + soltri = tri.Triangulation(mesh_node_x, mesh_node_y, triangles=state.body.mesh.elements) + interpol = tri.LinearTriInterpolator # if in3d else tri.CubicTriInterpolator + v = interpol(soltri, getattr(state, field)[:, 0]) + u = interpol(soltri, state.displacement[:, 0]) + + # higher resolution + x_disc = np.linspace(0, 1, 100) + y_disc = np.linspace(0, 1, 100) + node_x, node_y = np.meshgrid(x_disc, y_disc) + node_x = node_x.ravel() + node_y = node_y.ravel() + node_val = np.zeros_like(node_x) + # pylint: disable=consider-using-enumerate + for i in range(len(node_val)): + node_val[i] = u(node_x[i], node_y[i]) + + bound = np.linspace(0, 1, 100) + rev_bound = np.linspace(1, 0, 100) + bound_x = np.concatenate( + (bound, np.ones_like(bound), rev_bound, np.zeros_like(rev_bound), np.zeros(1)) + ) + bound_y = np.concatenate( + (np.zeros_like(bound), bound, np.ones_like(rev_bound), rev_bound, np.zeros(1)) + ) + bound_val = np.empty_like(bound_x) + # pylint: disable=consider-using-enumerate + for i in range(len(bound_x)): + bound_val[i] = u(bound_x[i], bound_y[i]) + + if in3d: + ax.view_init(elev=elev, azim=azim) + plot_surface( + fig, + ax, + node_x, + node_y, + node_val, + bound_x, + bound_y, + bound_val, + lambda x, y, z: v(x, y), + vmin, + vmax, + zmin, + zmax, + title, + ) + else: + plot_intersection(ax, u, x, ymin=zmin, ymax=zmax) + + +def plot_surface( + fig, + ax, + node_x, + node_y, + node_val, + bound_x, + bound_y, + bound_val, + v_func: Callable, + vmin=None, + vmax=None, + zmin=None, + zmax=None, + title=None, +): + p3dc = ax.plot_trisurf(node_x, node_y, node_val, alpha=1) + ax.set_zlim(zmin=zmin, zmax=zmax) + + mappable = map_colors(p3dc, v_func, "coolwarm", vmin, vmax) + plt.title(title) + + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_zlabel("u") + + ax.plot3D(bound_x, bound_y, bound_val, color="blue") + + cbar_ax = fig.add_axes([0.10, 0.957, 0.8, 0.02]) + plt.colorbar(mappable, shrink=0.6, aspect=15, cax=cbar_ax, orientation="horizontal") + + # plt.show() + + +def plot_intersection(ax, valfunc, insec_x, ymin, ymax): + space = np.linspace(0, 1, 100) + ax.set_ylim([ymin, ymax]) + ax.plot(space, valfunc(np.ones_like(space) * insec_x, space)) + + +def map_colors(p3dc, func, cmap="viridis", vmin=None, vmax=None): + """ + Color a tri-mesh according to a function evaluated in each barycentre. + + p3dc: a Poly3DCollection, as returned e.g. by ax.plot_trisurf + func: a single-valued function of 3 arrays: x, y, z + cmap: a colormap NAME, as a string + + Returns a ScalarMappable that can be used to instantiate a colorbar. + """ + + # reconstruct the triangles from internal data + # pylint: disable=protected-access + x, y, z, _ = p3dc._vec + slices = p3dc._segslices + triangles = np.array([np.array((x[s], y[s], z[s])).T for s in slices]) + + # compute the barycentres for each triangle + xb, yb, zb = triangles.mean(axis=1).T + + # compute the function in the barycentres + values = func(xb, yb, zb) + + # usual stuff + norm = plt.Normalize(vmin=vmin, vmax=vmax) + colors = get_cmap(cmap)(norm(values)) + + # set the face colors of the Poly3DCollection + p3dc.set_fc(colors) + + # if the caller wants a colorbar, they need this + return ScalarMappable(cmap=cmap, norm=norm) diff --git a/conmech/properties/body_properties.py b/conmech/properties/body_properties.py index 370208a0..22d62ac9 100644 --- a/conmech/properties/body_properties.py +++ b/conmech/properties/body_properties.py @@ -65,5 +65,5 @@ class ElasticRelaxationProperties(ElasticProperties, RelaxationBodyProperties): @dataclass -class BaseFunctionNormBodyProperties(BodyProperties): - pass +class MembraneProperties(BodyProperties): + propagation: float diff --git a/conmech/properties/mesh_description.py b/conmech/properties/mesh_description.py index f9877820..4dfaa558 100644 --- a/conmech/properties/mesh_description.py +++ b/conmech/properties/mesh_description.py @@ -1,6 +1,6 @@ from abc import ABC from dataclasses import dataclass -from typing import List +from typing import List, Optional import numpy as np import meshio @@ -11,7 +11,7 @@ @dataclass class MeshDescription(ABC): - initial_position: np.ndarray + initial_position: Optional[np.ndarray] def build(self): raise NotImplementedError() diff --git a/conmech/scenarios/problems.py b/conmech/scenarios/problems.py index f3c90729..c9589705 100644 --- a/conmech/scenarios/problems.py +++ b/conmech/scenarios/problems.py @@ -8,29 +8,26 @@ import numpy as np +from conmech.dynamics.contact.contact_law import ContactLaw +from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw from conmech.mesh.boundaries_description import BoundariesDescription from conmech.properties.mesh_description import MeshDescription -# pylint: disable=too-many-ancestors - - -class ContactLaw: - @staticmethod - def potential_normal_direction(u_nu: float) -> float: - raise NotImplementedError() - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - raise NotImplementedError() +from conmech.dynamics.statement import ( + Statement, + StaticDisplacementStatement, + QuasistaticVelocityStatement, + DynamicVelocityWithTemperatureStatement, + TemperatureStatement, + PiezoelectricStatement, + DynamicVelocityStatement, + QuasistaticVelocityWithPiezoelectricStatement, + QuasistaticRelaxationStatement, + StaticPoissonStatement, + WaveStatement, +) - @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - raise NotImplementedError() +# pylint: disable=too-many-ancestors @dataclass @@ -39,6 +36,14 @@ class Problem(ABC): mesh_descr: MeshDescription boundaries: BoundariesDescription + @classmethod + def statement(cls, body) -> Optional[Statement]: + raise NotImplementedError() + + @classmethod + def second_statement(cls, body) -> Optional[Statement]: + return None + @staticmethod def inner_forces( x: np.ndarray, *, v: Optional[np.ndarray] = None, t: Optional[float] = None @@ -69,19 +74,54 @@ class DynamicProblem(TimeDependentProblem, ABC): @dataclass -class PoissonProblem(StaticProblem, ABC): # TODO: rename +class PoissonProblem(StaticProblem, ABC): # pylint: disable=unused-argument + @classmethod + def statement(cls, body) -> Optional[Statement]: + return None + + @classmethod + def second_statement(cls, body) -> Optional[Statement]: + return StaticPoissonStatement(body) + @staticmethod def initial_temperature(x: np.ndarray) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros(len(x)) @staticmethod def internal_temperature(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros(len(x)) @staticmethod def outer_temperature(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros(len(x)) + + +@dataclass +class WaveProblem(DynamicProblem, ABC): + propagation: float + + @classmethod + def statement(cls, body) -> Statement: + return WaveStatement(body) + + @staticmethod + def initial_displacement(x: np.ndarray) -> np.ndarray: + return np.zeros_like(x) + + @staticmethod + def initial_velocity(x: np.ndarray) -> np.ndarray: + return np.zeros_like(x) + + +@dataclass +class ContactWaveProblem(WaveProblem, ABC): + contact_law: ContactLaw + + +@dataclass +class InteriorContactWaveProblem(ContactWaveProblem, ABC): + contact_law: InteriorContactLaw @dataclass @@ -95,15 +135,11 @@ class DisplacementProblem(Problem, ABC): def initial_displacement(x: np.ndarray) -> np.ndarray: return np.zeros_like(x) - @staticmethod - def friction_bound(u_nu: float) -> float: - raise NotImplementedError() - class StaticDisplacementProblem(DisplacementProblem, StaticProblem, ABC): - @staticmethod - def friction_bound(u_nu: float) -> float: - raise NotImplementedError() + @classmethod + def statement(cls, body) -> Statement: + return StaticDisplacementStatement(body) class TimeDependentDisplacementProblem(DisplacementProblem, ABC): @@ -117,47 +153,69 @@ def initial_velocity(x: np.ndarray) -> np.ndarray: class QuasistaticDisplacementProblem(QuasistaticProblem, TimeDependentDisplacementProblem, ABC): - pass + @classmethod + def statement(cls, body) -> Statement: + return QuasistaticVelocityStatement(body) class DynamicDisplacementProblem(DynamicProblem, TimeDependentDisplacementProblem, ABC): - pass + @classmethod + def statement(cls, body) -> Statement: + return DynamicVelocityStatement(body) class TemperatureTimeDependentProblem(TimeDependentDisplacementProblem, ABC): + contact_law_2: ContactLaw thermal_expansion: np.ndarray thermal_conductivity: np.ndarray + @classmethod + def second_statement(cls, body) -> Optional[Statement]: + return TemperatureStatement(body) + @staticmethod def initial_temperature(x: np.ndarray) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros(len(x)) class PiezoelectricTimeDependentProblem(TimeDependentDisplacementProblem, ABC): + contact_law_2: ContactLaw piezoelectricity: np.ndarray permittivity: np.ndarray + @classmethod + def second_statement(cls, body) -> Optional[Statement]: + return PiezoelectricStatement(body) + @staticmethod def initial_electric_potential(x: np.ndarray) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros(len(x)) class RelaxationQuasistaticProblem(QuasistaticDisplacementProblem, ABC): relaxation: Callable[[float], np.ndarray] + @classmethod + def statement(cls, body) -> Statement: + return QuasistaticRelaxationStatement(body) + @staticmethod def initial_absement(x: np.ndarray) -> np.ndarray: - return np.zeros_like(len(x)) + return np.zeros_like(x) class TemperatureDynamicProblem(DynamicProblem, TemperatureTimeDependentProblem, ABC): - pass + @classmethod + def statement(cls, body) -> Statement: + return DynamicVelocityWithTemperatureStatement(body) class PiezoelectricQuasistaticProblem( QuasistaticDisplacementProblem, PiezoelectricTimeDependentProblem, ABC ): - pass + @classmethod + def statement(cls, body) -> Statement: + return QuasistaticVelocityWithPiezoelectricStatement(body) class PiezoelectricDynamicProblem( diff --git a/conmech/simulations/problem_solver.py b/conmech/simulations/problem_solver.py index 077e8aff..13d9037d 100644 --- a/conmech/simulations/problem_solver.py +++ b/conmech/simulations/problem_solver.py @@ -7,18 +7,7 @@ import numpy as np from conmech.dynamics.dynamics import Dynamics -from conmech.dynamics.statement import ( - StaticDisplacementStatement, - QuasistaticVelocityStatement, - DynamicVelocityWithTemperatureStatement, - TemperatureStatement, - PiezoelectricStatement, - DynamicVelocityStatement, - QuasistaticVelocityWithPiezoelectricStatement, - QuasistaticRelaxationStatement, - StaticPoissonStatement, - Variables, -) +from conmech.dynamics.statement import Variables from conmech.mesh.mesh import Mesh from conmech.properties.body_properties import ( BodyProperties, @@ -27,27 +16,23 @@ ElasticRelaxationProperties, ViscoelasticPiezoelectricProperties, ViscoelasticTemperatureProperties, + MembraneProperties, ) from conmech.scenarios.problems import ( TimeDependentProblem, Problem, - StaticProblem, - QuasistaticProblem, PoissonProblem, - DisplacementProblem, StaticDisplacementProblem, TimeDependentDisplacementProblem, TemperatureTimeDependentProblem, - TemperatureDynamicProblem, PiezoelectricTimeDependentProblem, - PiezoelectricQuasistaticProblem, RelaxationQuasistaticProblem, - ContactLaw, + WaveProblem, ) from conmech.solvers import SchurComplementOptimization from conmech.solvers._solvers import SolversRegistry from conmech.solvers.solver import Solver -from conmech.solvers.validator import Validator +from conmech.state.products.product import Product from conmech.state.state import State, TemperatureState, PiezoelectricState @@ -89,12 +74,12 @@ def __init__(self, problem: Problem, body_properties: BodyProperties): self.problem: Problem = problem + # True if only one dimension of displacement is used + self.driving_vector = False + self.coordinates: Optional[str] = None self.step_solver: Optional[Solver] = None self.second_step_solver: Optional[Solver] = None - self.validator: Optional[Validator] = None - - self.penetration = [] self.done = 0 self.to_do = 1 @@ -116,89 +101,61 @@ def __set_step_solver(self, value): solver_class: Type[Solver] = SolversRegistry.get_by_name( solver_name=value, problem=self.problem ) - contact_law: Optional[ContactLaw] - friction_bound: Optional[Callable[[float], float]] - - # TODO: #65 fixed solvers to avoid: th_coef, ze_coef = mu_coef, la_coef - if isinstance(self.problem, DisplacementProblem): - contact_law = self.problem.contact_law - friction_bound = self.problem.friction_bound - if isinstance(self.problem, StaticProblem): - statement = StaticDisplacementStatement(self.body) - elif isinstance(self.problem, TimeDependentProblem): - if isinstance(self.problem, PiezoelectricQuasistaticProblem): - statement = QuasistaticVelocityWithPiezoelectricStatement(self.body) - elif isinstance(self.problem, RelaxationQuasistaticProblem): - statement = QuasistaticRelaxationStatement(self.body) - elif isinstance(self.problem, QuasistaticProblem): - statement = QuasistaticVelocityStatement(self.body) - elif isinstance(self.problem, TemperatureDynamicProblem): - statement = DynamicVelocityWithTemperatureStatement(self.body) - else: - statement = DynamicVelocityStatement(self.body) - else: - raise ValueError(f"Unsupported problem class: {self.problem.__class__}") - elif isinstance(self.problem, PoissonProblem): + statement = self.problem.statement(self.body) + if statement is None: self.step_solver = None - self.validator = None return - else: - raise ValueError(f"Unsupported problem class: {self.problem.__class__}") + self.step_solver = solver_class( statement, self.body, self.time_step, - contact_law, - friction_bound, + getattr(self.problem, "contact_law", None), + driving_vector=self.driving_vector, ) - self.validator = Validator(self.step_solver) def __set_second_step_solver(self, value): second_solver_class: Type[Solver] = SolversRegistry.get_by_name( solver_name=value, problem=self.problem ) - if isinstance(self.problem, TemperatureTimeDependentProblem): - self.second_step_solver = second_solver_class( - TemperatureStatement(self.body), - self.body, - self.time_step, - self.problem.contact_law, - self.problem.friction_bound, - ) - elif isinstance(self.problem, PiezoelectricTimeDependentProblem): - self.second_step_solver = second_solver_class( - PiezoelectricStatement(self.body), - self.body, - self.time_step, - self.problem.contact_law, - self.problem.friction_bound, - ) - elif isinstance(self.problem, PoissonProblem): - self.second_step_solver = second_solver_class( - StaticPoissonStatement(self.body), - self.body, - self.time_step, - self.problem.contact_law if hasattr(self.problem, "contact_law") else None, # TODO - None, - ) - else: + statement = self.problem.second_statement(self.body) + if statement is None: self.second_step_solver = None + return + + self.second_step_solver = second_solver_class( + statement, + self.body, + self.time_step, + getattr(self.problem, "contact_law_2", None), + driving_vector=self.driving_vector, + ) def solve(self, **kwargs): raise NotImplementedError() - def run(self, state: State, n_steps: int, verbose: bool = False, **kwargs): + def run( + self, + state: State, + n_steps: int, + verbose: bool = False, + products: Optional[List[Product]] = None, + **kwargs, + ): """ :param state: :param n_steps: number of steps :param verbose: show prints """ + if products is not None: + for prod in products: + state.inject_product(prod) + for _ in range(n_steps): self.step_solver.current_time += self.step_solver.time_step solution = self.find_solution( state, - self.validator, verbose=verbose, **kwargs, ) @@ -220,17 +177,13 @@ def run(self, state: State, n_steps: int, verbose: bool = False, **kwargs): else: raise ValueError(f"Unknown coordinates: {self.coordinates}") - self.penetration.append((state.time, state.penetration)) self.step_solver.iterate() self.done += 1 print(f"{self.done / self.to_do * 100:.2f}%", end="\r") - def find_solution(self, state, validator, *, verbose=False, **kwargs) -> np.ndarray: - quality = 0 - initial_guess = state[self.coordinates].reshape(state.body.mesh.dimension, -1) + def find_solution(self, state, **kwargs) -> np.ndarray: + initial_guess = state[self.coordinates].T.ravel().reshape(state.body.mesh.dimension, -1) solution = self.step_solver.solve(initial_guess, **kwargs) - # quality = validator.check_quality(state, solution, quality) - self.print_iteration_info(quality, validator.error_tolerance, verbose) return solution def find_solution_uzawa(self, solution, solution_t) -> Tuple[np.ndarray, np.ndarray]: @@ -306,14 +259,6 @@ def find_solution_uzawa(self, solution, solution_t) -> Tuple[np.ndarray, np.ndar return solution, solution_t - @staticmethod - def print_iteration_info(quality: float, error_tolerance: float, verbose: bool) -> None: - qualitative = quality > error_tolerance - sign = ">" if qualitative else "<" - end = "." if qualitative else ", trying again..." - if verbose: - print(f"quality = {quality} {sign} {error_tolerance}{end}") - class PoissonSolver(ProblemSolver): def __init__(self, problem: Problem, solving_method: str): @@ -322,10 +267,9 @@ def __init__(self, problem: Problem, solving_method: str): :param problem: :param solving_method: 'schur', 'optimization', 'direct' """ - body_prop = ElasticProperties( + body_prop = MembraneProperties( mass_density=1.0, - mu=0, - lambda_=0, + propagation=1.0, ) super().__init__(problem, body_prop) @@ -339,6 +283,7 @@ def solve(self, **kwargs) -> TemperatureState: state = TemperatureState(self.body) self.second_step_solver.t_vector[:] = state.temperature.ravel().copy() + self.second_step_solver.iterate() initial_guess = state[self.coordinates] solution = self.second_step_solver.solve(initial_guess=initial_guess, **kwargs) @@ -380,6 +325,7 @@ def solve(self, *, initial_displacement: Callable, verbose: bool = False, **kwar ) self.step_solver.u_vector[:] = state.displacement.T.ravel().copy() + self.step_solver.iterate() self.run(state, n_steps=1, verbose=verbose, **kwargs) return state @@ -444,6 +390,7 @@ def solve( self.step_solver.b_vector[:] = state.absement.T.ravel().copy() self.step_solver.u_vector[:] = state.displacement.T.ravel().copy() + self.step_solver.iterate() output_step = np.diff(output_step) results = [] @@ -510,6 +457,7 @@ def solve( self.step_solver.u_vector[:] = state.displacement.T.ravel().copy() self.step_solver.v_vector[:] = state.velocity.T.ravel().copy() + self.step_solver.iterate() output_step = np.diff(output_step) results = [] @@ -591,6 +539,9 @@ def solve( self.second_step_solver.v_vector[:] = state.velocity.T.ravel().copy() self.second_step_solver.t_vector[:] = state.temperature.T.ravel().copy() + self.step_solver.iterate() + self.second_step_solver.iterate() + output_step = np.diff(output_step) done = 0 for n in output_step: @@ -600,8 +551,6 @@ def solve( self.step_solver.current_time += self.step_solver.time_step self.second_step_solver.current_time += self.second_step_solver.time_step - # solution = self.find_solution(self.step_solver, state, solution, self.validator, - # verbose=verbose) solution, solution_t = self.find_solution_uzawa(solution, solution_t) if self.coordinates == "velocity": @@ -611,7 +560,6 @@ def solve( time=self.step_solver.current_time, ) state.set_temperature(solution_t) - # self.step_solver.iterate(solution) else: raise ValueError(f"Unknown coordinates: {self.coordinates}") yield state.copy() @@ -684,6 +632,9 @@ def solve( self.second_step_solver.v_vector[:] = state.velocity.T.ravel().copy() self.second_step_solver.p_vector[:] = state.electric_potential.T.ravel().copy() + self.step_solver.iterate() + self.second_step_solver.iterate() + output_step = np.diff(output_step) results = [] done = 0 @@ -708,3 +659,73 @@ def solve( results.append(state.copy()) return results + + +class WaveSolver(ProblemSolver): + def __init__( + self, + problem: WaveProblem, + solving_method: str, + ): + """Solves general Contact Mechanics problem. + + :param solving_method: 'schur', 'optimization', 'direct' + """ + body_prop = MembraneProperties(mass_density=1, propagation=problem.propagation) + super().__init__(problem, body_prop) + + _ = State(self.body) # TODO + + self.coordinates = "velocity" + self.driving_vector = True + self.solving_method = solving_method + + # super class method takes **kwargs, so signatures are consistent + # pylint: disable=arguments-differ + def solve( + self, + *, + n_steps: int, + initial_displacement: Callable, + initial_velocity: Callable, + state: Optional[State] = None, + output_step: Optional[iter] = None, + **kwargs, + ) -> List[State]: + """ + :param n_steps: number of time-step in simulation + :param output_step: from which time-step we want to get copy of State, + default (n_steps-1,) + example: for Setup.time-step = 2, n_steps = 10, + output_step = (2, 6, 9) we get 3 shared copy of + State for time-steps 4, 12 and 18 + :param initial_displacement: for the solver + :param initial_velocity: for the solver + :param state: if present, simulation starts from this state, + use this to resume stopped simulation. + :return: state + """ + output_step = (0, *output_step) if output_step else (0, n_steps) # 0 for diff + + if state is None: + state = State(self.body) + state.displacement[:] = initial_displacement( + self.body.mesh.nodes[: self.body.mesh.nodes_count] + ) + state.velocity[:] = initial_velocity(self.body.mesh.nodes[: self.body.mesh.nodes_count]) + + self.step_solver.current_time = state.time + + self.step_solver.u_vector[:] = state.displacement.T.ravel().copy() + self.step_solver.v_vector[:] = state.velocity.T.ravel().copy() + self.step_solver.iterate() + + output_step = np.diff(output_step) + results = [] + self.done = 0 + self.to_do = n_steps + for n in output_step: + self.run(state, n_steps=n, **kwargs) + results.append(state.copy()) + + return results diff --git a/conmech/solvers/direct.py b/conmech/solvers/direct.py index 64aa1976..382775ac 100644 --- a/conmech/solvers/direct.py +++ b/conmech/solvers/direct.py @@ -8,13 +8,14 @@ import scipy.optimize from conmech.dynamics.statement import Statement -from conmech.scenarios.problems import ContactLaw +from conmech.dynamics.contact.contact_law import DirectContactLaw from conmech.scene.body_forces import BodyForces from conmech.solvers._solvers import SolversRegistry from conmech.solvers.solver import Solver from conmech.solvers.solver_methods import make_equation +@SolversRegistry.register("dynamic", "direct") @SolversRegistry.register("static", "direct") class Direct(Solver): def __init__( @@ -22,23 +23,26 @@ def __init__( statement: Statement, body: BodyForces, time_step: float, - contact_law: Optional[ContactLaw] = None, - friction_bound: Optional[Callable[[float], float]] = None, + contact_law: Optional[DirectContactLaw] = None, + driving_vector: bool = False, ): super().__init__( statement, body, time_step, contact_law, - friction_bound, + driving_vector, ) self.equation: Optional[Callable] = None if contact_law is not None: self.equation = make_equation( jn=contact_law.subderivative_normal_direction, - jt=contact_law.regularized_subderivative_tangential_direction, - h_functional=friction_bound, + contact=( + contact_law.general_contact_condition + if hasattr(contact_law, "general_contact_condition") + else None + ), ) def __str__(self) -> str: @@ -46,25 +50,43 @@ def __str__(self) -> str: @property def node_relations(self) -> np.ndarray: - return self.statement.left_hand_side + return self.statement.left_hand_side.data @property def node_forces(self) -> np.ndarray: return self.statement.right_hand_side - def _solve_impl(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: + def _solve_impl( + self, + initial_guess: np.ndarray, + *, + variable_old: np.ndarray, + displacement: np.ndarray, + **kwargs, + ) -> np.ndarray: + displacement = np.squeeze(displacement.copy().reshape(1, -1)) if self.equation is not None: result = scipy.optimize.fsolve( self.equation, initial_guess, args=( + variable_old, self.body.mesh.nodes, self.body.mesh.contact_boundary, self.body.mesh.boundaries.contact_normals, self.node_relations, self.node_forces, + displacement, + self.body.dynamics.acceleration_operator.SM1.data, + self.time_step, ), ) else: result = np.linalg.solve(self.node_relations, self.node_forces) + result_len = len(result) + var_len = len(initial_guess.ravel()) + if result_len < var_len: + result_ = np.zeros(var_len) + result_[:result_len] = result[:] + result = result_ return np.asarray(result) diff --git a/conmech/solvers/optimization/global_optimization.py b/conmech/solvers/optimization/global_optimization.py index 96c71271..d76faab1 100644 --- a/conmech/solvers/optimization/global_optimization.py +++ b/conmech/solvers/optimization/global_optimization.py @@ -15,7 +15,7 @@ def __str__(self): @property def lhs(self) -> np.ndarray: - return self.statement.left_hand_side + return self.statement.left_hand_side.data @property def rhs(self) -> np.ndarray: diff --git a/conmech/solvers/optimization/optimization.py b/conmech/solvers/optimization/optimization.py index e6d3f47a..19f84d24 100644 --- a/conmech/solvers/optimization/optimization.py +++ b/conmech/solvers/optimization/optimization.py @@ -11,13 +11,13 @@ from conmech.dynamics.statement import ( Statement, - TemperatureStatement, - PiezoelectricStatement, StaticPoissonStatement, + WaveStatement, ) -from conmech.scenarios.problems import ContactLaw +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw +from conmech.dynamics.contact.interior_contact_law import InteriorContactLaw from conmech.solvers.solver import Solver -from conmech.solvers.solver_methods import make_cost_functional +from conmech.solvers.solver_methods import make_cost_functional, make_equation class Optimization(Solver): @@ -26,50 +26,39 @@ def __init__( statement: Statement, body: "Body", time_step: float, - contact_law: Optional[ContactLaw], - friction_bound, + contact_law: Optional[PotentialOfContactLaw], + driving_vector, ): super().__init__( statement, body, time_step, contact_law, - friction_bound, + driving_vector, ) - if statement.dimension >= 2: # TODO + self.loss = make_cost_functional( + normal_condition=contact_law.potential_normal_direction, + normal_condition_bound=contact_law.normal_bound, + tangential_condition=contact_law.potential_tangential_direction, + tangential_condition_bound=contact_law.tangential_bound, + variable_dimension=statement.dimension_out, + problem_dimension=statement.dimension_in, + ) + if isinstance(statement, WaveStatement): + if isinstance(contact_law, InteriorContactLaw): + self.loss = make_equation( # TODO! + jn=None, + contact=contact_law.potential_normal_direction, + ) + + if isinstance(statement, StaticPoissonStatement): self.loss = make_cost_functional( - normal_condition=contact_law.potential_normal_direction, - tangential_condition=( - contact_law.potential_tangential_direction - if hasattr(contact_law, "potential_tangential_direction") - else None + normal_condition=( + contact_law.potential_normal_direction if contact_law is not None else None ), - tangential_condition_bound=friction_bound, - variable_dimension=statement.dimension, - problem_dimension=body.mesh.dimension, + variable_dimension=statement.dimension_out, + problem_dimension=statement.dimension_in, ) - elif isinstance(statement, TemperatureStatement): - self.loss = make_cost_functional( - tangential_condition=contact_law.h_temp, - normal_condition=contact_law.temp_exchange, - normal_condition_bound=-1, - ) - elif isinstance(statement, PiezoelectricStatement): - self.loss = make_cost_functional( - tangential_condition=contact_law.electric_charge_tangetial, - tangential_condition_bound=-1, - normal_condition=None, - variable_dimension=statement.dimension, - problem_dimension=body.mesh.dimension, - ) - elif isinstance(statement, StaticPoissonStatement): - self.loss = make_cost_functional( - normal_condition=contact_law.potential_normal_direction, - variable_dimension=statement.dimension, - problem_dimension=body.mesh.dimension, - ) - else: - raise ValueError(f"Unknown statement: {statement}") def __str__(self): raise NotImplementedError() @@ -86,7 +75,7 @@ def _solve_impl( self, initial_guess: np.ndarray, *, - velocity: np.ndarray, + variable_old: np.ndarray, displacement: np.ndarray, method="BFGS", fixed_point_abs_tol: float = math.inf, @@ -100,12 +89,14 @@ def _solve_impl( maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9)) tol = kwargs.get("tol", 1e-12) args = ( + variable_old, self.body.mesh.nodes, self.body.mesh.contact_boundary, self.body.mesh.boundaries.contact_normals, self.lhs, self.rhs, displacement, + np.ascontiguousarray(self.body.dynamics.acceleration_operator.SM1.data), self.time_step, ) @@ -115,7 +106,7 @@ def _solve_impl( loss.append(self.loss(solution, *args)[0]) while norm >= fixed_point_abs_tol: - if method.lower() in ( # TODO + if method.lower() in ( "quasi secant method", "limited memory quasi secant method", "quasi secant method limited memory", @@ -126,7 +117,7 @@ def _solve_impl( from kosopt import qsmlm solution = qsmlm.minimize(self.loss, solution, args=args, maxiter=maxiter) - sols[self.loss(solution, *args)] = solution + sols.append(solution.copy()) elif method.lower() == "constrained": contact_nodes_count = self.body.mesh.boundaries.contact_nodes_count @@ -168,4 +159,5 @@ def constr(x): old_solution = solution.copy() min_index = loss.index(np.min(loss)) solution = sols[min_index] + return solution diff --git a/conmech/solvers/optimization/schur_complement.py b/conmech/solvers/optimization/schur_complement.py index 4942aded..58658815 100644 --- a/conmech/solvers/optimization/schur_complement.py +++ b/conmech/solvers/optimization/schur_complement.py @@ -20,14 +20,14 @@ def __init__( body, time_step, contact_law, - friction_bound, + driving_vector, ): super().__init__( statement, body, time_step, contact_law, - friction_bound, + driving_vector, ) self.contact_ids = slice(0, body.mesh.contact_nodes_count) @@ -44,8 +44,8 @@ def __init__( def recalculate_displacement(self): return calculate_schur_complement_matrices( - matrix=self.statement.left_hand_side, - dimension=self.statement.dimension, + matrix=self.statement.left_hand_side.data, + dimension=self.statement.dimension_in, contact_indices=self.contact_ids, free_indices=self.free_ids, ) @@ -53,13 +53,13 @@ def recalculate_displacement(self): def recalculate_forces(self): node_forces, forces_free = calculate_schur_complement_vector( vector=self.statement.right_hand_side, - dimension=self.statement.dimension, + dimension=self.statement.dimension_in, contact_indices=self.contact_ids, free_indices=self.free_ids, contact_x_free=self.contact_x_free, free_x_free_inverted=self.free_x_free_inverted, ) - if self.statement.dimension == 2: + if self.statement.dimension_in == 2: return node_forces.T, forces_free return node_forces.reshape(-1), forces_free.reshape(-1) @@ -81,16 +81,35 @@ def _solve_impl( fixed_point_abs_tol: float = math.inf, **kwargs, ) -> np.ndarray: - truncated_initial_guess = self.truncate_free_nodes(initial_guess) - solution_contact = super()._solve_impl( - truncated_initial_guess, fixed_point_abs_tol=fixed_point_abs_tol, **kwargs - ) + # initial_guess: [[xc, xf, xd] [yc, yf, yd]] + truncated_ig = self.truncate_free_nodes(initial_guess) + # truncated_ig: [[xc, yc]] + if self.driving_vector: + truncated_ig = truncated_ig[:, : truncated_ig.shape[1] // 2] + # truncated_ig: [[xc]] + solution_contact = super()._solve_impl(truncated_ig, **kwargs) + # solution_contact: [xc, yc] / [xc] solution_free = self.complement_free_nodes(solution_contact) + # solution_free: [[xf], [xd], [yf], [yd]] / [xf, xd] + if self.driving_vector: + length = len(solution_free) + extender = np.zeros(length * 2).reshape(-1, 1) + solution_free = solution_free.reshape(-1, 1) + extender[:length, 0] = solution_free[:, 0] + solution_free = extender + + length = len(solution_contact) + extender = np.zeros(length * 2) + extender[:length] = solution_contact[:] + solution_contact = extender + # solution_free: [[xf,], [xd,], [0,], [0,]] + # solution_contact: [xc, 0] solution = self.merge(solution_contact, solution_free) + # solution: [xc, xf, xd, yc, yf, yd] return solution def truncate_free_nodes(self, initial_guess: np.ndarray) -> np.ndarray: - if self.statement.dimension == 2: + if self.statement.dimension_in == 2 or self.driving_vector: _result = initial_guess.reshape(2, -1) _result = _result[:, self.contact_ids] _result = _result.reshape(1, -1) @@ -99,7 +118,7 @@ def truncate_free_nodes(self, initial_guess: np.ndarray) -> np.ndarray: return initial_guess[self.contact_ids] def complement_free_nodes(self, truncated_solution: np.ndarray) -> np.ndarray: - if self.statement.dimension == 2: + if self.statement.dimension_in == 2: _result = truncated_solution.reshape(-1, 1) else: _result = truncated_solution @@ -110,7 +129,7 @@ def complement_free_nodes(self, truncated_solution: np.ndarray) -> np.ndarray: return result def merge(self, solution_contact: np.ndarray, solution_free: np.ndarray) -> np.ndarray: - if self.statement.dimension == 2: + if self.statement.dimension_in == 2 or self.driving_vector: u_contact = solution_contact.reshape(2, -1) u_free = solution_free.reshape(2, -1) _result = np.concatenate((u_contact, u_free), axis=1) diff --git a/conmech/solvers/solver.py b/conmech/solvers/solver.py index 5d3817d2..9b0fdd06 100644 --- a/conmech/solvers/solver.py +++ b/conmech/solvers/solver.py @@ -2,12 +2,12 @@ Created at 18.02.2021 """ -from typing import Callable, Optional +from typing import Optional import numpy as np from conmech.dynamics.statement import Statement, Variables -from conmech.scenarios.problems import ContactLaw +from conmech.dynamics.contact.contact_law import ContactLaw class Solver: @@ -17,10 +17,9 @@ def __init__( body: "Body", time_step: float, contact_law: Optional[ContactLaw] = None, - friction_bound: Optional[Callable[[float], float]] = None, + driving_vector: bool = False, ): self.contact_law: Optional[ContactLaw] = contact_law - self.friction_bound: Optional[Callable[[float], float]] = friction_bound self.body = body self.statement: Statement = statement @@ -35,6 +34,8 @@ def __init__( self.elasticity: np.ndarray = body.dynamics.elasticity + self.driving_vector = driving_vector + self.statement.update( Variables( absement=self.b_vector, @@ -54,20 +55,28 @@ def iterate(self): pass def _solve_impl( - self, initial_guess, *, velocity: np.ndarray, displacement: np.ndarray, **kwargs + self, + initial_guess, + *, + variable_old: np.ndarray, + displacement: np.ndarray, + **kwargs, ): raise NotImplementedError() def solve(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray: solution = self._solve_impl( - initial_guess, velocity=self.v_vector, displacement=self.u_vector, **kwargs + initial_guess, + variable_old=self.v_vector, + displacement=self.u_vector, + **kwargs, # TODO: FIXME variable_old ) for dirichlet_cond in self.statement.find_dirichlet_conditions(): c = self.body.mesh.boundaries.boundaries[dirichlet_cond].node_condition node_count = self.body.mesh.nodes_count for i, j in self.body.mesh.boundaries.get_all_boundary_indices( - dirichlet_cond, node_count, self.statement.dimension + dirichlet_cond, node_count, self.statement.dimension_in ): solution[i] = c[j] diff --git a/conmech/solvers/solver_methods.py b/conmech/solvers/solver_methods.py index c1eb7366..ee763124 100644 --- a/conmech/solvers/solver_methods.py +++ b/conmech/solvers/solver_methods.py @@ -11,34 +11,59 @@ @numba.njit(inline="always") -def interpolate_node_between(edge, vector, dimension): +def interpolate_node_between(edge, vector, full_vector, dimension): result = np.zeros(dimension) offset = len(vector) // dimension + offset_full = len(full_vector) // dimension for node in edge: for i in range(dimension): if node < offset: # exclude dirichlet nodes (and inner nodes in schur) result[i] += 0.5 * vector[i * offset + node] + else: + # old values + result[i] += 0.5 * full_vector[i * offset_full + node] return result +# pylint: disable=too-many-arguments def make_equation( jn: Optional[callable], - jt: Optional[callable], - h_functional: Optional[callable], + contact: Optional[callable] = None, problem_dimension=2, ) -> callable: - # TODO Make it prettier - if jn is None: + if jn is None or contact is not None: + contact = njit(contact, value=0.0) + # pylint: disable=unused-argument @numba.njit - def equation(u_vector: np.ndarray, _, __, lhs: np.ndarray, rhs: np.ndarray) -> np.ndarray: - result = np.dot(lhs, u_vector) - rhs + def equation( + var: np.ndarray, + var_old: np.ndarray, + _, + __, + ___, + lhs: np.ndarray, + rhs: np.ndarray, + displacement: np.ndarray, + volume_multiplier: np.ndarray, + time_step, + ) -> np.ndarray: + ind = lhs.shape[0] + response = np.zeros(ind) + for i in range(ind): + response[i] = contact(var[i], displacement[i], time_step) + res = ( + 0.5 * np.dot(np.dot(lhs, var[:ind]), var[:ind]) + - np.dot(rhs, var[:ind]) + + 0.5 * np.dot(np.dot(volume_multiplier, response), np.ones_like(var[:ind])) + + np.dot(var[ind:], var[ind:].T) + ) + + result = np.asarray(res).ravel() return result else: jn = numba.njit(jn) - jt = numba.njit(jt) - h_functional = numba.njit(h_functional) @numba.njit() def contact_part(u_vector, nodes, contact_boundary, contact_normals): @@ -51,30 +76,12 @@ def contact_part(u_vector, nodes, contact_boundary, contact_normals): normal_vector = contact_normals[ei] # ASSUMING `u_vector` and `nodes` have the same order! - um = interpolate_node_between(edge, u_vector, dimension=problem_dimension) + um = interpolate_node_between(edge, u_vector, u_vector, dimension=problem_dimension) um_normal = (um * normal_vector).sum() - um_tangential = um - um_normal * normal_vector - - v_tau_0 = np.asarray( - [ - 1 - normal_vector[0] * normal_vector[0], - 0 - normal_vector[0] * normal_vector[1], - ] - ) - v_tau_1 = np.asarray( - [ - 0 - normal_vector[1] * normal_vector[0], - 1 - normal_vector[1] * normal_vector[1], - ] - ) edge_len = nph.length(edge, nodes) - j_x = edge_len * 0.5 * (jn(um_normal, normal_vector[0])) + h_functional( - um_normal - ) * jt(um_tangential, v_tau_0) - j_y = edge_len * 0.5 * (jn(um_normal, normal_vector[1])) + h_functional( - um_normal - ) * jt(um_tangential, v_tau_1) + j_x = edge_len * 0.5 * (jn(um_normal, normal_vector[0], 0.0)) + j_y = edge_len * 0.5 * (jn(um_normal, normal_vector[1], 0.0)) if n_id_0 < offset: contact_vector[n_id_0] += j_x @@ -86,17 +93,22 @@ def contact_part(u_vector, nodes, contact_boundary, contact_normals): return contact_vector + # pylint: disable=unused-argument,too-many-arguments @numba.njit def equation( - u_vector: np.ndarray, - vertices: np.ndarray, + var: np.ndarray, + var_old: np.ndarray, + nodes: np.ndarray, contact_boundary: np.ndarray, contact_normals: np.ndarray, lhs: np.ndarray, rhs: np.ndarray, + displacement: np.ndarray, + base_integrals, + time_step, ) -> np.ndarray: - c_part = contact_part(u_vector, vertices, contact_boundary, contact_normals) - result = np.dot(lhs, u_vector) + c_part - rhs + c_part = contact_part(var, nodes, contact_boundary, contact_normals) + result = np.dot(lhs, var) + c_part - rhs return result return equation @@ -107,7 +119,7 @@ def njit(func: Optional[Callable], value: Optional[Any] = 0) -> Callable: ret_val = func if func is not None else value @numba.njit() - def const(_): + def const(_, __, ___): return ret_val return const @@ -132,7 +144,9 @@ def contact_cost(length, normal, normal_bound, tangential, tangential_bound): return length * (normal_bound * normal + tangential_bound * tangential) @numba.njit() - def contact_cost_functional(var, displacement, nodes, contact_boundary, contact_normals): + def contact_cost_functional( + var, var_old, static_displacement, nodes, contact_boundary, contact_normals, dt + ): offset = len(var) // problem_dimension cost = 0.0 @@ -141,7 +155,7 @@ def contact_cost_functional(var, displacement, nodes, contact_boundary, contact_ edge = contact_boundary[ei] normal_vector = contact_normals[ei] # ASSUMING `u_vector` and `nodes` have the same order! - vm = interpolate_node_between(edge, var, dimension=variable_dimension) + vm = interpolate_node_between(edge, var, var_old, dimension=variable_dimension) if variable_dimension == 1: vm_normal = vm[0] vm_tangential = np.empty(0) @@ -149,8 +163,16 @@ def contact_cost_functional(var, displacement, nodes, contact_boundary, contact_ vm_normal = (vm * normal_vector).sum() vm_tangential = vm - vm_normal * normal_vector - um = interpolate_node_between(edge, displacement, dimension=problem_dimension) - um_normal = (um * normal_vector).sum() + static_displacement_mean = interpolate_node_between( + edge, + static_displacement, + static_displacement, + dimension=problem_dimension, + ) + static_displacement_normal = (static_displacement_mean * normal_vector).sum() + static_displacement_tangential = ( + static_displacement_mean - static_displacement_normal * normal_vector + ) for node_id in edge: if node_id >= offset: @@ -158,18 +180,32 @@ def contact_cost_functional(var, displacement, nodes, contact_boundary, contact_ cost += contact_cost( nph.length(edge, nodes), - normal_condition(vm_normal), - normal_condition_bound(um_normal), - tangential_condition(vm_tangential), - tangential_condition_bound(um_normal), + normal_condition(vm_normal, static_displacement_normal, dt), + normal_condition_bound(vm_normal, static_displacement_normal, dt), + tangential_condition(vm_tangential, static_displacement_tangential, dt), + tangential_condition_bound(vm_normal, static_displacement_normal, dt), ) return cost - # pylint: disable=unused-argument # 'dt' + # pylint: disable=too-many-arguments,unused-argument # 'base_integrals' @numba.njit() - def cost_functional(var, nodes, contact_boundary, contact_normals, lhs, rhs, u_vector, dt): - ju = contact_cost_functional(var, u_vector, nodes, contact_boundary, contact_normals) - result = 0.5 * np.dot(np.dot(lhs, var), var) - np.dot(rhs, var) + ju + def cost_functional( + var, + var_old, + nodes, + contact_boundary, + contact_normals, + lhs, + rhs, + u_vector, + base_integrals, + dt, + ): + ju = contact_cost_functional( + var, var_old, u_vector, nodes, contact_boundary, contact_normals, dt + ) + ind = lhs.shape[0] + result = 0.5 * np.dot(np.dot(lhs, var[:ind]), var[:ind]) - np.dot(rhs, var[:ind]) + ju result = np.asarray(result).ravel() return result diff --git a/conmech/solvers/validator.py b/conmech/solvers/validator.py deleted file mode 100644 index daee11e0..00000000 --- a/conmech/solvers/validator.py +++ /dev/null @@ -1,49 +0,0 @@ -""" -Created at 18.02.2021 -""" - -import numpy as np - -from conmech.solvers.solver import Solver -from conmech.solvers.solver_methods import make_equation -from conmech.state.state import State - - -class Validator: - def __init__(self, solver: Solver, error_tolerance: float = 1): - self.error_tolerance: float = error_tolerance - self.elasticity: np.ndarray = solver.elasticity - self.rhs: callable - if solver.contact_law is None: - self.rhs = make_equation(None, None, None) - else: - self.rhs = make_equation( - jn=solver.contact_law.subderivative_normal_direction, - jt=solver.contact_law.regularized_subderivative_tangential_direction, - h_functional=solver.friction_bound, - ) - - def validate(self, state: State, solution: np.ndarray) -> float: - quality_inv = np.linalg.norm( - self.rhs( - solution, - state.body.mesh.nodes, - state.body.mesh.contact_boundary, - state.body.mesh.boundaries.contact_normals, - self.elasticity, - state.body.dynamics.force.integrate(time=state.time), - ) - ) - if quality_inv == 0: - quality = np.inf - else: - quality = quality_inv**-1 - return quality - - def check_quality( - self, state: State, solution: np.ndarray, previous_quality: float = None - ) -> float: - quality = self.validate(state, solution) - if previous_quality is not None and previous_quality == quality: - raise RuntimeError("Can't find a solution! ") - return quality diff --git a/conmech/state/products/__init__.py b/conmech/state/products/__init__.py new file mode 100644 index 00000000..8961ee91 --- /dev/null +++ b/conmech/state/products/__init__.py @@ -0,0 +1,18 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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. diff --git a/conmech/state/products/intersection_contact_limit_points.py b/conmech/state/products/intersection_contact_limit_points.py new file mode 100644 index 00000000..d9f2cfe2 --- /dev/null +++ b/conmech/state/products/intersection_contact_limit_points.py @@ -0,0 +1,67 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 numpy as np +from matplotlib import tri +import scipy.optimize as opt + +from conmech.state.products.product import Product + + +class VerticalIntersectionContactLimitPoints(Product): + def __init__(self, x, obstacle_level): + super().__init__(f"limit points at {x:.2f}") + self.level = x + self.obstacle_level = obstacle_level + + def update(self, state): + y_min = np.min(state.body.mesh.nodes[:, 1]) + y_max = np.max(state.body.mesh.nodes[:, 1]) + step = len(state.body.mesh.nodes) ** -1 / 2 + u = interpolate(state, "displacement") + + def u_intsec(y_coord): + return u(self.level, y_coord) - self.obstacle_level + + self.data[state.time] = estimate_zeros(u_intsec, y_min, y_max, step) + + +def interpolate(state, field): + assert state.body.mesh.dimension == 2 # membrane have to be 2D + X = state.body.mesh.nodes[:, 0] + Y = state.body.mesh.nodes[:, 1] + + soltri = tri.Triangulation(X, Y, triangles=state.body.mesh.elements) + interpol = tri.LinearTriInterpolator + return interpol(soltri, getattr(state, field)[:, 0]) + + +def estimate_zeros(func, xmin, xmax, step) -> tuple: + zeros = [] + x = xmin + sign = np.sign(func(x)) + while x <= xmax: + f_x = func(x) + if f_x == 0: + zeros.append(x) + elif f_x < 0 < sign or f_x > 0 > sign: + zeros.append(opt.brentq(func, x - step, x)) + sign = np.sign(f_x) + x += step + return tuple(zeros) diff --git a/conmech/state/products/penetration.py b/conmech/state/products/penetration.py new file mode 100644 index 00000000..9929e841 --- /dev/null +++ b/conmech/state/products/penetration.py @@ -0,0 +1,39 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 numpy as np + +from conmech.state.products.product import Product + + +class Penetration(Product): + """ + This class assume foundation equals x=0. + """ + + def __init__(self): + super().__init__("penetration") + + def update(self, state): + if len(state.displaced_nodes[state.body.mesh.contact_indices, 1]) != 0: + all_p = state.displaced_nodes[state.body.mesh.contact_indices, 1] + p = np.min(all_p) + else: + p = 0.0 + self.data[state.time] = p diff --git a/conmech/state/products/product.py b/conmech/state/products/product.py new file mode 100644 index 00000000..cd1525ca --- /dev/null +++ b/conmech/state/products/product.py @@ -0,0 +1,42 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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. + + +class Product: + def __init__(self, name): + self.name = name + self.data = {} + + def update(self, state): + raise NotImplementedError() + + def __copy__(self): + return self.copy() + + def copy(self) -> "Product": + result = Product(self.name) + result.data = self.data.copy() + return result + + def range(self, start, stop) -> "Product": + result = self.copy() + for key in self.data.keys(): + if key <= start or key > stop: + del result.data[key] + return result diff --git a/conmech/state/products/verticalintersection.py b/conmech/state/products/verticalintersection.py new file mode 100644 index 00000000..15f6daea --- /dev/null +++ b/conmech/state/products/verticalintersection.py @@ -0,0 +1,42 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 numpy as np +from matplotlib import tri + +from conmech.state.products.product import Product + + +class VerticalIntersection(Product): + def __init__(self, x): + super().__init__(f"intersection at {x:.2f}") + self.level = x + + def update(self, state): + y_min = np.min(state.body.mesh.nodes[:, 1]) + y_max = np.max(state.body.mesh.nodes[:, 1]) + step = len(state.body.mesh.nodes) ** -1 / 2 + space = np.linspace(y_min, y_max, int((y_max - y_min) // step)) + + X = state.body.mesh.nodes[:, 0] + Y = state.body.mesh.nodes[:, 1] + soltri = tri.Triangulation(X, Y, triangles=state.body.mesh.elements) + u = tri.LinearTriInterpolator(soltri, state.displacement[:, 0]) + + self.data[state.time] = (space, u(np.ones_like(space) * self.level, space)) diff --git a/conmech/state/state.py b/conmech/state/state.py index e0513d1f..a4646c8c 100644 --- a/conmech/state/state.py +++ b/conmech/state/state.py @@ -2,6 +2,8 @@ Created at 18.02.2021 """ +import pickle + import numpy as np @@ -10,6 +12,8 @@ def __init__(self, body): self.body = body self.body.state = self + self.products = {} + self.absement: np.ndarray = np.zeros((self.body.mesh.nodes_count, self.body.mesh.dimension)) self.displacement: np.ndarray = np.zeros( (self.body.mesh.nodes_count, self.body.mesh.dimension) @@ -19,10 +23,15 @@ def __init__(self, body): self.setup = None self.__stress: np.ndarray = None self.constitutive_law = None + self.temperature = None # TODO self.time: float = 0 def set_displacement( - self, displacement_vector: np.ndarray, time: float, *, update_absement: bool = False + self, + displacement_vector: np.ndarray, + time: float, + *, + update_absement: bool = False, ): self.displacement = displacement_vector.reshape((self.body.mesh.dimension, -1)).T self.displaced_nodes[: self.body.mesh.nodes_count, :] = ( @@ -32,6 +41,7 @@ def set_displacement( dt = time - self.time self.absement += dt * self.displacement self.time = time + self.update_products() def set_velocity(self, velocity_vector: np.ndarray, time: float, *, update_displacement: bool): self.velocity = velocity_vector.reshape((self.body.mesh.dimension, -1)).T @@ -42,6 +52,7 @@ def set_velocity(self, velocity_vector: np.ndarray, time: float, *, update_displ self.body.mesh.nodes[: self.body.mesh.nodes_count, :] + self.displacement[:, :] ) self.time = time + self.update_products() @property def stress(self): @@ -70,14 +81,12 @@ def stress_x(self): def stress_y(self): return self.stress[:, 1, 1] - @property - def penetration(self): - """ - This method assume foundation equals x=0. - """ - if len(self.displaced_nodes[self.body.mesh.contact_indices, 1]) != 0: - return np.min(self.displaced_nodes[self.body.mesh.contact_indices, 1]) - return 0 + def inject_product(self, product): + self.products[product.name] = product + + def update_products(self): + for prod in self.products.values(): + prod.update(self) def __getitem__(self, item: [int, str]) -> np.ndarray: if item in (0, "displacement"): @@ -97,8 +106,32 @@ def __copy__(self) -> "State": copy.displaced_nodes[:] = self.displaced_nodes copy.velocity[:] = self.velocity copy.time = self.time + copy.products = {n: p.copy() for n, p in self.products.items()} return copy + def save(self, path): + fos = self.body.dynamics.force.outer.source + self.body.dynamics.force.outer.source = None + fis = self.body.dynamics.force.inner.source + self.body.dynamics.force.inner.source = None + s = self.setup + self.setup = None + cl = self.constitutive_law + self.constitutive_law = None + + with open(path, "wb+") as file: + pickle.dump(self, file) + + self.body.dynamics.force.outer.source = fos + self.body.dynamics.force.inner.source = fis + self.setup = s + self.constitutive_law = cl + + @classmethod + def load(cls, path): + with open(path, "rb") as file: + return pickle.load(file) + class TemperatureState(State): def __init__(self, body): diff --git a/conmech/struct/__init__.py b/conmech/struct/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/conmech/struct/stiffness_matrix.py b/conmech/struct/stiffness_matrix.py new file mode 100644 index 00000000..c36ae8a8 --- /dev/null +++ b/conmech/struct/stiffness_matrix.py @@ -0,0 +1,89 @@ +class StiffnessMatrix: + DIMENSION = None + + def __init__(self, data): + self.data = data + + def __iadd__(self, other): + if isinstance(other, StiffnessMatrix): + assert self.DIMENSION == other.DIMENSION + self.data += other.data + else: + self.data += other + + def __add__(self, other): + if isinstance(other, StiffnessMatrix): + assert self.DIMENSION == other.DIMENSION + return type(self)(self.data + other.data) + return type(self)(self.data + other) + + def __radd__(self, other): + return self + other + + def __imul__(self, other): + self.data *= other + + def __mul__(self, other): + return type(self)(self.data * other) + + def __rmul__(self, other): + return self * other + + def copy(self): + return type(self)(self.data.copy()) + + def __getitem__(self, item): + return self.data[item] + + def __matmul__(self, other): + if isinstance(other, StiffnessMatrix): + return self.data @ other.data + return self.data @ other + + def __imatmul__(self, other): + if isinstance(other, StiffnessMatrix): + self.data @= other.data + else: + self.data @= other + + def __rmatmul__(self, other): + return other @ self.data + + # pylint: disable=invalid-name + @property + def T(self): + return self.data.T + + +class SM1(StiffnessMatrix): + DIMENSION = (1, 1) + + +class SM1to2(StiffnessMatrix): + DIMENSION = (1, 2) + + +class SM1to3(StiffnessMatrix): + DIMENSION = (1, 3) + + +class SM2(StiffnessMatrix): + DIMENSION = (2, 2) + + # pylint: disable=invalid-name + @property + def SM1(self) -> SM1: + x_len = self.data.shape[0] // 2 + y_len = self.data.shape[1] // 2 + return SM1(self.data[:x_len, :y_len]) + + +class SM3(StiffnessMatrix): + DIMENSION = (3, 3) + + # pylint: disable=invalid-name + @property + def SM1(self) -> SM1: + x_len = self.data.shape[0] // 3 + y_len = self.data.shape[1] // 3 + return SM1(self.data[:x_len, :y_len]) diff --git a/examples/BOSK_2024.py b/examples/BOSK_2024.py new file mode 100644 index 00000000..090df6f1 --- /dev/null +++ b/examples/BOSK_2024.py @@ -0,0 +1,280 @@ +from dataclasses import dataclass +from typing import Optional + +import numpy as np + +from conmech.dynamics.contact.damped_normal_compliance import make_damped_norm_compl +from conmech.helpers.config import Config +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.plotting.membrane import plot_in_columns, plot_limit_points +from conmech.scenarios.problems import ContactWaveProblem, InteriorContactWaveProblem +from conmech.simulations.problem_solver import WaveSolver +from conmech.properties.mesh_description import CrossMeshDescription +from conmech.state.products.verticalintersection import VerticalIntersection +from conmech.state.products.intersection_contact_limit_points import ( + VerticalIntersectionContactLimitPoints, +) +from conmech.state.state import State + + +@dataclass() +class BoundaryContactMembrane(ContactWaveProblem): + time_step: ... = 1 / 64 + propagation: ... = np.sqrt(4.0) + contact_law: ... = make_damped_norm_compl(obstacle_level=1.0, kappa=10.0, beta=0.1)() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([100]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0,) or x[1] in (0, 1), + contact=lambda x: x[0] == 1, + ) + + +@dataclass() +class InteriorContactMembrane(InteriorContactWaveProblem): + time_step: ... = 1 / 512 + propagation: ... = np.sqrt(100.00) + contact_law: ... = make_damped_norm_compl( + obstacle_level=1.0, kappa=10.0, beta=100.0, interior=True + )() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([2500]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription(dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1)) + + +def runner(config: Config, setup, name, steps, intersect, continuing, solving_method): + print(name) + output_step = ( + int(steps * 0.25), + int(steps * 0.5), + int(steps * 0.75), + steps, + ) + + to_simulate = False + if config.force: + to_simulate = True + else: + for step in output_step: + try: + State.load(f"{config.path}/{name}_step_{step}") + except IOError: + to_simulate = True + + if to_simulate: + runner = WaveSolver(setup, solving_method) + + if continuing: + path = f"{config.path}/{continuing}" + print("Loading:", path) + state = State.load(path) + state.products = {} + else: + state = None + + states = runner.solve( + n_steps=steps, + output_step=output_step, + products=[ + VerticalIntersectionContactLimitPoints( + obstacle_level=setup.contact_law.obstacle_level, x=intersect + ), + VerticalIntersection(x=intersect), + ], + initial_displacement=setup.initial_displacement, + initial_velocity=setup.initial_velocity, + state=state, + verbose=True, + ) + for step, state in enumerate(states): + state.save(f"{config.path}/{name}_step_{step}") + else: + states = [] + for step in output_step: + states.append(State.load(f"{config.path}/{name}_step_{step}")) + + if continuing: + path = f"{config.path}/{continuing}" + print("Loading:", path) + state = State.load(path) + plot_limit_points( + state.products[f"limit points at {intersect:.2f}"], + title=rf"$\kappa={setup.contact_law.kappa}$ " rf"$\beta={setup.contact_law.beta}$", + finish=False, + ) + plot_limit_points( + states[-1].products[f"limit points at {intersect:.2f}"], + title=rf"$\kappa={setup.contact_law.kappa}$, $\beta={setup.contact_law.beta}$", + finish=config.show, + ) + + states_ids = list(range(len(states))) + to_plot = states_ids + vmin = np.inf + vmax = -np.inf + field = "velocity" + zmin = np.inf + zmax = -np.inf + zfield = "displacement" + for i, state in enumerate(states): + if i not in to_plot: + continue + vmax = max(max(getattr(state, field)[:, 0]), vmax) + vmin = min(min(getattr(state, field)[:, 0]), vmin) + zmax = max(max(getattr(state, zfield)[:, 0]), zmax) + zmin = min(min(getattr(state, zfield)[:, 0]), zmin) + prec = 1 + zmax = round(zmax + 0.05 * prec, prec) + zmin = round(zmin - 0.05 * prec, prec) + states_ = [] + for i, state in enumerate(states): + if i not in to_plot: + continue + states_.append(state) + plot_in_columns( + states_, + field=field, + vmin=vmin, + vmax=vmax, + zmin=zmin, + zmax=zmax, + x=intersect, + title=f"intersection at {intersect:.2f}", + finish=config.show, + ) + plot_in_columns( + states_, + field=field, + vmin=vmin, + vmax=vmax, + zmin=zmin, + zmax=zmax, + in3d=True, + title=f"velocity", + finish=config.show, + ) + + +def boundary_contact(config: Config): + precision = 16 if not config.test else 3 + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=1 / precision, scale=[1, 1] + ) + default = BoundaryContactMembrane(mesh_descr) + T = 4.0 if not config.test else default.time_step * 2 + + setups = dict() + to_simulate = [ + "no contact", + "plain", + "beta=0.00", + ] + + setup = default + setups["plain"] = setup + + setup = BoundaryContactMembrane(mesh_descr) + setup.contact_law = make_damped_norm_compl( + obstacle_level=default.contact_law.obstacle_level, kappa=0.00, beta=0.0, interior=False + )() + setups["no contact"] = setup + + setup = BoundaryContactMembrane(mesh_descr) + setup.contact_law = make_damped_norm_compl( + obstacle_level=default.contact_law.obstacle_level, + kappa=default.contact_law.kappa, + beta=0.0, + interior=False, + )() + setups["beta=0.00"] = setup + + for name in to_simulate: + runner( + config, + setups[name], + name="case_2_bdry_" + name, + steps=int(T / setups[name].time_step) + 1, + intersect=1.00, + continuing=None, + solving_method="schur", + ) + + +def interior_contact(config: Config): + precision = 4 if not config.test else 3 + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=1 / precision, scale=[1, 1] + ) + default = InteriorContactMembrane(mesh_descr) + T = 0.5 if not config.test else default.time_step * 2 + + setups = dict() + to_simulate = [ + "plain", + "beta=0.00", + "beta=1000.0", + ] + + setup = default + setups["plain"] = setup + + setup = InteriorContactMembrane(mesh_descr) + setup.contact_law = make_damped_norm_compl( + obstacle_level=default.contact_law.obstacle_level, + kappa=default.contact_law.kappa, + beta=0.0, + interior=True, + )() + setups["beta=0.00"] = setup + + setup = InteriorContactMembrane(mesh_descr) + setup.contact_law = make_damped_norm_compl( + obstacle_level=default.contact_law.obstacle_level, + kappa=default.contact_law.kappa, + beta=1000.0, + interior=True, + )() + setups["beta=1000.0"] = setup + + for name in to_simulate: + runner( + config, + setups[name], + name="case_1_int_" + name, + steps=int(T / setups[name].time_step) + 1, + intersect=0.50, + continuing=None, + solving_method="global", + ) + + +def main(config: Config): + """ + Entrypoint to example. + + To see result of simulation you need to call from python `main(Config().init())`. + """ + boundary_contact(config) + interior_contact(config) + + +if __name__ == "__main__": + main(Config().init()) diff --git a/examples/BOSK_plot_postprocessing.py b/examples/BOSK_plot_postprocessing.py new file mode 100644 index 00000000..4d6c30a7 --- /dev/null +++ b/examples/BOSK_plot_postprocessing.py @@ -0,0 +1,255 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman-Szwarc +# +# 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 numpy as np +from matplotlib import pyplot as plt + +from conmech.plotting.membrane import plot_limit_points +from conmech.state.state import State +from matplotlib.ticker import StrMethodFormatter + + +def case1(): + states = {} + for kappa, beta in ( + (0.0, 0.0), + (0.0, 0.25), + (0.0, 0.5), + (0.0, 0.75), + (0.0, 1.0), + (1.0, 0.0), + (1.0, 0.5), + (10.0, 0.5), + (100.0, 0.5), + )[:]: + print((kappa, beta)) + states[kappa, beta] = State.load(f"output/BOSK.pub/c1_kappa={kappa:.2f};beta={beta:.2f}") + + c1_reference(states, output_path="output/BOSK.pub") + c1_steady_state(states, output_path="output/BOSK.pub") + c1_influ(states, output_path="output/BOSK.pub") + + +def show(output_path, name): + plt.gca().yaxis.set_major_formatter(StrMethodFormatter("{x:,.2f}")) # 2 decimal places + plt.gca().xaxis.set_major_formatter(StrMethodFormatter("{x:,.2f}")) # 2 decimal places + if output_path is None: + plt.show() + else: + plt.savefig(output_path + f"/{name}.png", format="png", dpi=800) + + +def c1_reference(states, output_path=None): + plt.rc("axes", axisbelow=True) + + kappa = 100.0 + beta = 0.5 + timesteps = (0.875, 2.0) + labels = ("b)", "a)") + eps = 0.001 + + all_zeros = [] + + for subnumber, timestep in zip(labels, timesteps): + intersect = states[kappa, beta].products["intersection at 1.00"] + plt.figure(figsize=(5, 4)) + zeros = states[kappa, beta].products["limit points at 1.00"].data[timestep] + all_zeros.extend(zeros) + plt.xticks((0.0, *zeros, 1.0), rotation=90) + plt.yticks((0.0, 1.0, 1.8)) + plt.grid() + plt.axhspan(1.0, 1.8, alpha=0.1, color="lightskyblue") + + for t, v in intersect.data.items(): + if timestep - eps > t or t > timestep + eps: + continue + print(t) + plt.plot(*v, color=f"black") + states[kappa, beta].products["limit points at 1.00"].range(0.00, 8.00) + + plt.scatter(zeros, np.ones_like(zeros), color=f"black", s=10) + plt.ylim(0.0, 1.8) + plt.xlabel("y") + plt.ylabel("z") + plt.title(rf"$\kappa={kappa:.2f}$, $\beta={beta:.2f}$, $t={timestep}$") + plt.title(subnumber, loc="left") + + show(output_path, name=f"boundary_{timestep:.3f}") + + plt.figure(figsize=(10, 4)) + plot_limit_points( + states[kappa, beta].products["limit points at 1.00"].range(0.00, 8.00), + title=None, + color=f"black", + finish=False, + ) + + plt.xticks((0.0, *timesteps, 8.0)) + plt.yticks((0.0, *all_zeros, 1.0)) + plt.grid() + plt.xlabel("time") + plt.ylabel("y") + plt.title(rf"$\kappa={kappa:.2f}$, $\beta={beta:.2f}$") + plt.title("c)", loc="left") + + show(output_path, name="reference") + + +def c1_steady_state(states, output_path=None): + plt.figure(figsize=(10, 4)) + kappa = 0.0 + beta = 0.0 + plot_limit_points( + states[kappa, beta].products["limit points at 1.00"].range(0.00, 8.00), + title=None, + color=f"gray", + finish=False, + label=rf"$\beta={beta:.2f}$", + ) + beta = 1.0 + plot_limit_points( + states[kappa, beta].products["limit points at 1.00"].range(0.00, 8.00), + title=None, + color=f"salmon", + finish=False, + label=rf"$\beta={beta:.2f}$", + ) + plt.xticks((0.0, 0.9, 1.8, 2.7, 3.6, 4.5, 5.4, 6.3, 7.2, 8.0)) + plt.yticks((0.0, 1.0)) + plt.xlabel("time") + plt.ylabel("y") + plt.title(rf"$\kappa={kappa:.2f}$") + + plt.legend(loc="center right") + show(output_path, name="steady_state") + + +def c1_influ(states, output_path=None): + cases = { + ("a)", "kappa"): ( + ("lightskyblue", 100.0, 0.5), + ("yellowgreen", 10.0, 0.5), + ("gold", 1.0, 0.5), + ("salmon", 0.0, 0.5), + ), + ("b)", "beta"): ( + ("lightskyblue", 0.0, 0.25), + ("yellowgreen", 0.0, 0.5), + ("gold", 0.0, 0.75), + ("salmon", 0.0, 1.0), + ), + } + for (subnumber, var), case in cases.items(): + plt.figure(figsize=(6, 4.5)) + for c, kappa, beta in case: + print((kappa, beta)) + variable = kappa if var == "kappa" else beta + plot_limit_points( + states[kappa, beta].products["limit points at 1.00"].range(0.00, 4.00), + title=None, + label=rf"$\{var}={variable:.2f}$", + finish=False, + color=f"{c}", + ) + + plt.legend(loc="center right") + plt.xticks((0.0, 0.92, 1.8, 2.65, 3.6, 4.0)) + plt.yticks((0.0, 0.5, 1.0)) + plt.grid() + plt.xlabel("time") + plt.ylabel("y") + const_name = "kappa" if var == "beta" else "beta" + const_value = kappa if var == "beta" else beta + plt.title(rf"$\{const_name}={const_value:.2f}$") + plt.title(subnumber, loc="left") + + show(output_path, name="var_" + var) + + +def case2(): + output_path = "output/BOSK.pub" + kappa = 1.0 + # for var, subnumber, beta in (('control', 'a)', 0.), ('beta', 'b)', 100.)): + # plt.figure(figsize=(6, 4.5)) + # state = State.load(f"output/BOSK.pub/i2_kappa={kappa:.2f};beta={beta:.2f}") + # plot_limit_points( + # state.products['limit points at 0.50'], + # title=fr'$\kappa={kappa}$ $\beta={beta}$', finish=False) + # state = State.load(f"output/BOSK.pub/ci2_kappa={kappa:.2f};beta={beta:.2f}") + # plot_limit_points( + # state.products['limit points at 0.50'], + # title=fr'$\kappa={kappa}$ $\beta={beta}$', finish=False) + # plt.title(subnumber, loc='left') + # plt.xticks((0.0, 3.6, 7.2, 10.8, 14.4, 16.0)) + # plt.yticks((0.0, 0.5, 1.0)) + # plt.xlabel("time") + # plt.ylabel("y") + # plt.grid() + # + # show(output_path, name='int_' + var) + + beta = 0.0 # TODO + state = State.load(f"output/BOSK.pub/i2_c2_plain") + plot_limit_points( + state.products["limit points at 0.50"].range(0.0, 4.0), + title=rf"$\kappa={kappa}$ $\beta={beta}$", + finish=False, + ) + plt.show() + + timesteps = ( + 1.0 - 0.0625, + 1.75, + ) + labels = ("b)", "a)") + eps = 0.001 + + all_zeros = [] + + for subnumber, timestep in zip(labels, timesteps): + intersect = state.products["intersection at 0.50"] + plt.figure(figsize=(5, 4)) + zeros = state.products["limit points at 0.50"].data[timestep] + all_zeros.extend(zeros) + plt.xticks((0.0, *zeros, 1.0), rotation=90) + plt.yticks((0.0, 1.0, 1.8)) + plt.grid() + plt.axhspan(1.0, 1.8, alpha=0.1, color="lightskyblue") + + for t, v in intersect.data.items(): + if timestep - eps > t or t > timestep + eps: + continue + print(t) + plt.plot(*v, color=f"black") + state.products["limit points at 0.50"].range(0.00, 8.00) + + plt.scatter(zeros, np.ones_like(zeros), color=f"black", s=10) + plt.ylim(0.0, 1.8) + plt.xlabel("y") + plt.ylabel("z") + plt.title(rf"$\kappa={kappa:.2f}$, $\beta={beta:.2f}$, $t={timestep}$") + plt.title(subnumber, loc="left") + + plt.show() + + # show(output_path, name=f'boundary_{timestep:.3f}') + + +if __name__ == "__main__": + # case1() + case2() diff --git a/examples/BOST_2024.py b/examples/BOST_2024.py index 8cdc5a6a..23f483b5 100644 --- a/examples/BOST_2024.py +++ b/examples/BOST_2024.py @@ -25,7 +25,8 @@ from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.scenarios.problems import StaticDisplacementProblem +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.simulations.problem_solver import StaticSolver as StaticProblemSolver from conmech.properties.mesh_description import BOST2023MeshDescription from examples.error_estimates import error_estimates @@ -35,27 +36,16 @@ kappa = 0.42 -class BOST23(ContactLaw): +class BOST23(PotentialOfContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: - raise NotImplementedError() - - @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return np.sum(u_tau * u_tau) ** 0.5 - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + return np.sum(var_tau * var_tau) ** 0.5 @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - return 0 + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + return 1.0 @dataclass @@ -72,12 +62,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu: float) -> float: - if u_nu < 0: - return 0 - return u_nu - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -86,23 +70,27 @@ def friction_bound(u_nu: float) -> float: def prepare_setup(ig, setup): if ig == "inf": - def potential_normal_direction(u_nu: float) -> float: + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: return 0 kwargs = {"method": "constrained"} else: - def potential_normal_direction(u_nu: float) -> float: - u_nu -= GAP + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + var_nu -= GAP # EXAMPLE 10 a = 0.1 b = 0.1 - if u_nu <= 0: + if var_nu <= 0: result = 0.0 - elif u_nu < b: - result = (a + np.exp(-b)) / (2 * b) * u_nu**2 + elif var_nu < b: + result = (a + np.exp(-b)) / (2 * b) * var_nu**2 else: - result = a * u_nu - np.exp(-u_nu) + ((b + 2) * np.exp(-b) - a * b) / 2 + result = a * var_nu - np.exp(-var_nu) + ((b + 2) * np.exp(-b) - a * b) / 2 return ig * result kwargs = {"method": "POWELL"} @@ -111,14 +99,20 @@ def potential_normal_direction(u_nu: float) -> float: return kwargs -def main(config: Config, igs: List[int], highlighted: Iterable): +def main( + config: Config, + igs: Optional[List[int]] = None, + highlighted: Optional[Iterable] = None, +): """ Entrypoint to example. To see result of simulation you need to call from python `main(Config().init())`. """ - PREFIX = "BOST_POWELL" + PREFIX = "BOST" to_simulate = [] + igs = igs or [1024, "inf"] + highlighted = highlighted or igs for ig in igs: try: with open(f"{config.outputs_path}/{PREFIX}_ig_{ig}", "rb") as output: @@ -129,7 +123,10 @@ def main(config: Config, igs: List[int], highlighted: Iterable): if config.force: to_simulate = igs - mesh_descr = BOST2023MeshDescription(initial_position=None, max_element_perimeter=1 / 32) + mesh_descr = BOST2023MeshDescription( + initial_position=None, + max_element_perimeter=1 / 32 if not config.test else 1 / 6, + ) if to_simulate: print("Simulating...") @@ -208,11 +205,13 @@ def main(config: Config, igs: List[int], highlighted: Iterable): X = [ig for ig in igs if not isinstance(ig, str)] Y = list(errors.values())[:-1] Y = [v[0] for v in Y] - Y = -np.asarray(Y) - plot_errors(X, Y, highlighted_id=None, save=f"{config.outputs_path}/convergence.pdf") + # Y = -np.asarray(Y) + print(f"{X=}") + print(f"{Y=}") + plot_errors(config, X, Y, highlighted_id=None, save=f"{config.outputs_path}/convergence.pdf") -def plot_errors(X, Y, highlighted_id, save: Optional[str] = None): +def plot_errors(config, X, Y, highlighted_id, save: Optional[str] = None): plt.plot(X[:], Y[:], marker="o", color="gray") if highlighted_id is not None: plt.plot(X[highlighted_id], Y[highlighted_id], "ro", color="black") @@ -220,131 +219,14 @@ def plot_errors(X, Y, highlighted_id, save: Optional[str] = None): plt.grid(True, which="major", linestyle="--", linewidth=0.5) plt.xlabel(r"$\lambda^{-1}_{\,\, n}$") plt.ylabel(r"$||\mathbf{u}^h_n - \mathbf{u}||$") - if save is None: + if save is None and config.show: plt.show() - else: + elif config.save: plt.savefig(save, format="pdf") if __name__ == "__main__": - # results from the paper highlighted = (1, 1024, 1216, 2048, 5931641, "inf") - highlighted_id = np.asarray((1, 20, 24, 26, -3)) - - X = np.asarray( - [ - 0, - 1, - 2, - 4, - 5, - 8, - 11, - 16, - 22, - 32, - 45, - 64, - 90, - 128, - 181, - 256, - 362, - 512, - 724, - 896, - 1024, - 1088, - 1120, - 1152, - 1216, - 1448, - 2048, - 2896, - 4096, - 5792, - 8192, - 11585, - 16384, - 23170, - 32768, - 46340, - 65536, - 92681, - 131072, - 185363, - 262144, - 370727, - 524288, - 741455, - 1048576, - 1482910, - 2097152, - 2965820, - 4194304, - 5931641, - 8388608, - 11863283, - ] - ) - Y = np.asarray( - [ - 1.957688239442841, - 1.955823800273701, - 1.9550737995244916, - 1.9535713010795799, - 1.9528202182876828, - 1.9505617554633947, - 1.9483005592237832, - 1.9445216877172742, - 1.939976356940407, - 1.9323603026049123, - 1.9224029295661968, - 1.9077076426739854, - 1.887348094110511, - 1.8570040907627985, - 1.8135169713254757, - 1.7493452687885334, - 1.6528260472436915, - 1.5016916054916745, - 1.247485589366432, - 0.9833463022514393, - 0.7219304912004216, - 0.5510492874873285, - 0.44668589874753256, - 0.32141797154720314, - 0.2235180524765496, - 0.18110766227417885, - 0.13217099205041302, - 0.09590338285207281, - 0.06933759603103229, - 0.050067258317548866, - 0.036147047768223, - 0.02613519925974203, - 0.018945653731469527, - 0.013792699500121302, - 0.01010466081623074, - 0.0074703782466904725, - 0.0055919695751946025, - 0.004256100793696483, - 0.003309243652330325, - 0.002641490469903024, - 0.002173650857273765, - 0.0018483357409417435, - 0.0016235440683907585, - 0.0014687173119601173, - 0.0013620445313966081, - 0.0012883643852057848, - 0.0012373372831507256, - 0.0012019371089261616, - 0.0011773413958895045, - 0.0011602208824840082, - 0.0011019352332271976, - 0.001139913361649654, - ] - ) - - plot_errors(X, Y, highlighted_id, save="convergence.pdf") show = True @@ -357,6 +239,6 @@ def plot_errors(X, Y, highlighted_id, save: Optional[str] = None): 1024 + 128 + 64, } igs.update(eigs) - igs = [0] + sorted(list(igs)) + ["inf"] # + igs = [0] + sorted(list(igs)) + ["inf"] igs = igs[:] - main(Config(save=not show, show=show, force=False).init(), igs, highlighted) + main(Config(save=not show, show=show, force=True).init(), igs, highlighted) diff --git a/examples/Jureczka_Ochal_Bartman_2023.py b/examples/Jureczka_Ochal_Bartman_2023.py index 7c4b6f1f..a0af98d9 100644 --- a/examples/Jureczka_Ochal_Bartman_2023.py +++ b/examples/Jureczka_Ochal_Bartman_2023.py @@ -12,36 +12,45 @@ from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import ContactLaw, QuasistaticDisplacementProblem -from conmech.simulations.problem_solver import TimeDependentSolver as QuasistaticProblemSolver +from conmech.scenarios.problems import QuasistaticDisplacementProblem +from conmech.dynamics.contact.contact_law import ContactLaw, PotentialOfContactLaw +from conmech.simulations.problem_solver import ( + TimeDependentSolver as QuasistaticProblemSolver, +) from conmech.properties.mesh_description import JOB2023MeshDescription from examples.utils import viscoelastic_constitutive_law -def make_contact_law(limit_value, limit): - class JureczkaOchalBartman2023(ContactLaw): +def make_contact_law(limit_value, limit, friction_bound): + class JureczkaOchalBartman2023(PotentialOfContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + n = friction_bound + b = 0.1 + if var_nu <= 0: return 0.0 - if u_nu < limit: - return limit_value * u_nu - return limit_value * limit - - @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return -0.3 * np.exp(-np.sum(u_tau * u_tau) ** 0.5) + 0.7 * np.sum(u_tau * u_tau) ** 0.5 - # return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1)\ + if var_nu < b: + return n * var_nu + return n * b @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + if var_nu <= 0: + return 0.0 + if var_nu < limit: + return limit_value * var_nu + return limit_value * limit @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float ) -> float: - return 0 + return ( + -0.3 * np.exp(-np.sum(var_tau * var_tau) ** 0.5) + + 0.7 * np.sum(var_tau * var_tau) ** 0.5 + ) return JureczkaOchalBartman2023 @@ -55,7 +64,7 @@ def regularized_subderivative_tangential_direction( path = "./output/JOB2023" -def make_setup(mesh_descr_, boundaries_, contact_law_, friction_bound_): +def make_setup(mesh_descr_, boundaries_, contact_law_): @dataclass() class QuasistaticSetup(QuasistaticDisplacementProblem): mu_coef: ... = 40 @@ -77,16 +86,6 @@ def outer_forces(x, t=None): return np.array([-4.0, -4.0]) return np.array([0.0, 0.0]) - @staticmethod - def friction_bound(u_nu: float) -> float: - n = friction_bound_ - b = 0.1 - if u_nu <= 0: - return 0.0 - if u_nu < b: - return n * u_nu - return n * b - boundaries: ... = boundaries_ return QuasistaticSetup(mesh_descr=mesh_descr_) @@ -117,9 +116,9 @@ def main(config: Config): contact=lambda x: x[1] == 0, dirichlet=lambda x: (x[0] - x1) ** 2 + (x[1] - y1) ** 2 <= (r + eps) ** 2, ) - soft_foundation = make_contact_law(300, 0.1) - hard_foundation = make_contact_law(3000, 0.1) - friction_bound = 3 + soft_foundation = make_contact_law(300, 0.1, friction_bound=3) + soft_foundation_friction = make_contact_law(300, 0.1, friction_bound=300) + hard_foundation = make_contact_law(3000, 0.1, friction_bound=3) simulate = config.force if not config.test: @@ -148,12 +147,11 @@ def main(config: Config): if name == "hard": contact_law = hard_foundation() if name == "friction": - friction_bound = 300 + contact_law = soft_foundation_friction() setup = make_setup( mesh_descr_=mesh_descr, boundaries_=boundaries, contact_law_=contact_law, - friction_bound_=friction_bound, ) runner = QuasistaticProblemSolver(setup, "schur") states = runner.solve( @@ -180,7 +178,6 @@ def main(config: Config): mesh_descr_=mesh_descr, boundaries_=boundaries, contact_law_=contact_law, - friction_bound_=friction_bound, ) if name == names[0]: steps = (0, *output_steps[1:]) diff --git a/examples/Jureczka_and_Ochal_2019.py b/examples/Jureczka_and_Ochal_2019.py index 1148a4a4..21fa6a34 100644 --- a/examples/Jureczka_and_Ochal_2019.py +++ b/examples/Jureczka_and_Ochal_2019.py @@ -1,46 +1,57 @@ -""" -Created at 21.08.2019 -""" - +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2024 Piotr Bartman-Szwarc +# +# 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 dataclasses import dataclass import numpy as np from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem +from conmech.scenarios.problems import StaticDisplacementProblem +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.simulations.problem_solver import StaticSolver as StaticProblemSolver from conmech.properties.mesh_description import CrossMeshDescription -class JureczkaOchal2019(ContactLaw): +class JureczkaOchal2019(PotentialOfContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + if var_nu <= 0: return 0.0 - if u_nu < 0.1: - return 10 * u_nu * u_nu + if var_nu < 0.1: + return 10 * var_nu * var_nu return 0.1 @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + return np.log(np.sum(var_tau * var_tau) ** 0.5 + 1) @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - return 0 - # regularization = 1 / np.sqrt(u_tau[0] * u_tau[0] + u_tau[1] * u_tau[1] + rho ** 2) - # result = regularization * (u_tau[0] * v_tau[0] + u_tau[1] * v_tau[1]) - # return result + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + if static_displacement_nu < 0: + return 0 + if static_displacement_nu < 0.1: + return 8 * static_displacement_nu + return 0.8 @dataclass @@ -57,14 +68,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu: float) -> float: - if u_nu < 0: - return 0 - if u_nu < 0.1: - return 8 * u_nu - return 0.8 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) diff --git a/examples/Sofonea_Ochal_Bartman_2023.py b/examples/Sofonea_Ochal_Bartman_2023.py index b6249346..306c139d 100644 --- a/examples/Sofonea_Ochal_Bartman_2023.py +++ b/examples/Sofonea_Ochal_Bartman_2023.py @@ -30,8 +30,9 @@ from conmech.scenarios.problems import RelaxationQuasistaticProblem from conmech.simulations.problem_solver import QuasistaticRelaxation from conmech.properties.mesh_description import SOB2023MeshDescription +from conmech.state.products.penetration import Penetration -from examples.p_slope_contact_law import make_const_contact_law +from conmech.dynamics.contact.constant_constact_law import make_const_contact_law from examples.utils import elastic_relaxation_constitutive_law eps = 1e-18 @@ -50,7 +51,7 @@ class QuasistaticSetup(RelaxationQuasistaticProblem): mu_coef: ... = E / (1 + kappa) la_coef: ... = E * kappa / ((1 + kappa) * (1 - 2 * kappa)) time_step: ... = 1 / 128 - contact_law: ... = make_const_contact_law(slope=10) + contact_law: ... = make_const_contact_law(resistance=10) @staticmethod def relaxation(t: float) -> np.ndarray: @@ -70,10 +71,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0.0, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[0] >= 4 and x[1] < eps, dirichlet=lambda x: x[0] <= 1 and x[1] < eps, @@ -84,7 +81,7 @@ def main(config: Config): """ Entrypoint to example. - To see result of simulation you need to call from python `main(Config().init())`. + To see result you need to call from python `main(Config().init())`. """ if not config.test: elements_number = (20, 20) @@ -175,11 +172,12 @@ def zero_relaxation(t=None): n_steps=examples[name]["n_steps"], output_step=examples[name]["output_steps"], verbose=False, + products=[Penetration()], initial_absement=setup.initial_absement, initial_displacement=setup.initial_displacement, - tol=1e-9 if config.test else 1e-3, - fixed_point_abs_tol=1e-9 if config.test else 1e-3, - method="Powell" if config.test else "BFGS", + tol=1e-9 if not config.test else 1e-3, + fixed_point_abs_tol=1e-9 if not config.test else 1e-3, + method="Powell" if not config.test else "BFGS", ) f_max = -np.inf f_min = np.inf @@ -205,7 +203,8 @@ def zero_relaxation(t=None): f"{config.outputs_path}/{name}_h_{h}_penetration", "wb+", ) as output: - pickle.dump(runner.penetration, output) + penetration = states[-1].products["penetration"] + pickle.dump(list((t, p) for t, p in penetration.data.items()), output) with open( f"{config.outputs_path}/{name}_h_{h}_global", "wb+", diff --git a/examples/error_estimates.py b/examples/error_estimates.py index 9b2e7acb..a074cc27 100644 --- a/examples/error_estimates.py +++ b/examples/error_estimates.py @@ -44,10 +44,20 @@ def compare(ref: TemperatureState, sol: TemperatureState): u2dx, u2dy = calculate_dx_dy(x0, u2_0, x1, u2_1, x2, u2_2) # tdx, tdy = calculate_dx_dy(x0, t0, x1, t1, x2, t2) u1hdx, u1hdy = calculate_dx_dy( - x0, u1hi(*x0).compressed(), x1, u1hi(*x1).compressed(), x2, u1hi(*x2).compressed() + x0, + u1hi(*x0).compressed(), + x1, + u1hi(*x1).compressed(), + x2, + u1hi(*x2).compressed(), ) u2hdx, u2hdy = calculate_dx_dy( - x0, u2hi(*x0).compressed(), x1, u2hi(*x1).compressed(), x2, u2hi(*x2).compressed() + x0, + u2hi(*x0).compressed(), + x1, + u2hi(*x1).compressed(), + x2, + u2hi(*x2).compressed(), ) # thdx, thdy = calculate_dx_dy(x0, thi(*x0), x1, thi(*x1), x2, thi(*x2)) diff --git a/examples/example_tarzia_problem.py b/examples/example_Tarzia_problem.py similarity index 81% rename from examples/example_tarzia_problem.py rename to examples/example_Tarzia_problem.py index 7f3d112c..555b080e 100644 --- a/examples/example_tarzia_problem.py +++ b/examples/example_Tarzia_problem.py @@ -6,17 +6,20 @@ from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import PoissonProblem, ContactLaw +from conmech.scenarios.problems import PoissonProblem +from conmech.dynamics.contact.contact_law import ContactLaw, PotentialOfContactLaw from conmech.simulations.problem_solver import PoissonSolver from conmech.properties.mesh_description import CrossMeshDescription def make_slope_contact_law(slope: float) -> Type[ContactLaw]: - class TarziaContactLaw(ContactLaw): + class TarziaContactLaw(PotentialOfContactLaw): @staticmethod - def potential_normal_direction(u_nu: float) -> float: + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: b = 5 - r = u_nu + r = var_nu # EXAMPLE 11 # if r < b: # result = (r - b) ** 2 @@ -27,25 +30,12 @@ def potential_normal_direction(u_nu: float) -> float: result *= slope return result - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - raise NotImplementedError() - - @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - raise NotImplementedError() - return TarziaContactLaw @dataclass() class StaticPoissonSetup(PoissonProblem): - contact_law: ... = make_slope_contact_law(slope=1000) + contact_law_2: ... = make_slope_contact_law(slope=1000) @staticmethod def internal_temperature(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: @@ -76,8 +66,8 @@ def main(config: Config): """ alphas = [0.01, 0.1, 1, 10, 100, 1000, 10000] ihs = [4, 8, 16, 32, 64, 128, 256] - alphas = alphas - ihs = ihs + alphas = alphas if not config.test else alphas[:1] + ihs = ihs if not config.test else ihs[:1] for alpha in alphas: for ih in ihs: @@ -125,7 +115,11 @@ def draw(config, alpha, ih): drawer.cmap = "plasma" drawer.field_name = "temperature" drawer.draw( - show=config.show, save=config.save, foundation=False, field_max=max_, field_min=min_ + show=config.show, + save=config.save, + foundation=False, + field_max=max_, + field_min=min_, ) diff --git a/examples/example_dynamic.py b/examples/example_dynamic.py index b9f66179..8b7ac757 100644 --- a/examples/example_dynamic.py +++ b/examples/example_dynamic.py @@ -12,7 +12,7 @@ from conmech.scenarios.problems import DynamicDisplacementProblem from conmech.simulations.problem_solver import TimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law @dataclass() @@ -35,10 +35,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - def main(config: Config): """ diff --git a/examples/example_dynamic_membrane.py b/examples/example_dynamic_membrane.py new file mode 100644 index 00000000..ac97d1ba --- /dev/null +++ b/examples/example_dynamic_membrane.py @@ -0,0 +1,81 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman-Szwarc +# +# 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 dataclasses import dataclass +from typing import Optional + +import numpy as np +from conmech.helpers.config import Config +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.plotting.drawer import Drawer +from conmech.scenarios.problems import WaveProblem +from conmech.simulations.problem_solver import WaveSolver +from conmech.properties.mesh_description import CrossMeshDescription + + +@dataclass() +class MembraneSetup(WaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([0.2]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription(dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1)) + + +def main(config: Config): + """ + Entrypoint to example. + + To see result of simulation you need to call from python `main(Config().init())`. + """ + max_element_perimeter = 1 / 8 if not config.test else 1 / 3 + mesh_descr = CrossMeshDescription( + initial_position=None, + max_element_perimeter=max_element_perimeter, + scale=[1, 1], + ) + setup = MembraneSetup(mesh_descr) + runner = WaveSolver(setup, "direct") + n_steps = 32 if not config.test else 3 + + states = runner.solve( + n_steps=n_steps, + output_step=(0, n_steps), + initial_displacement=setup.initial_displacement, + initial_velocity=setup.initial_velocity, + verbose=True, + ) + drawer = Drawer(state=states[-1], config=config) + drawer.draw( + show=config.show, + save=config.save, + foundation=False, + ) + + +if __name__ == "__main__": + main(Config().init()) diff --git a/examples/example_imported_mesh_static.py b/examples/example_imported_mesh_static.py index bd7f30dd..949cb494 100644 --- a/examples/example_imported_mesh_static.py +++ b/examples/example_imported_mesh_static.py @@ -13,7 +13,7 @@ from conmech.properties.mesh_description import ImportedMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law E = 10000 @@ -34,10 +34,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, -1]) if x[0] > 1.9 and x[1] < 0.1 else np.zeros(2) - @staticmethod - def friction_bound(u_nu: float) -> float: - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0 and x[0] < 0.5, dirichlet=lambda x: x[0] == 0 ) diff --git a/examples/example_nonhomogenous_density.py b/examples/example_nonhomogenous_density.py index 8af579d1..5671c15a 100644 --- a/examples/example_nonhomogenous_density.py +++ b/examples/example_nonhomogenous_density.py @@ -12,7 +12,7 @@ from conmech.simulations.problem_solver import NonHomogenousSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law @dataclass @@ -27,14 +27,10 @@ def inner_forces(x, t=None): @staticmethod def outer_forces(x, t=None): - return np.array([0, 0.1]) if x[1] < 0.2 and x[0] > 2.2 else np.array([0, 0]) - - @staticmethod - def friction_bound(u_nu: float) -> float: - return 0 + return np.array([0, 0.1]) if x[1] < 0.2 and x[0] >= 2.0 else np.array([0, 0]) boundaries: ... = BoundariesDescription( - contact=lambda x: x[1] == 0 and x[0] < 0.2, dirichlet=lambda x: x[0] == 0 + contact=lambda x: x[1] == 0 and x[0] <= 1.0, dirichlet=lambda x: x[0] == 0 ) diff --git a/examples/example_piezo_quasistatic.py b/examples/example_piezo_quasistatic.py index 4f4f6da9..d4153548 100644 --- a/examples/example_piezo_quasistatic.py +++ b/examples/example_piezo_quasistatic.py @@ -6,6 +6,7 @@ import numpy as np +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.simulations.problem_solver import PiezoelectricTimeDependentSolver @@ -13,25 +14,31 @@ from conmech.plotting.drawer import Drawer from conmech.properties.mesh_description import Barboteu2008MeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law -class PPSlopeContactLaw(make_slope_contact_law(slope=1e1)): +class PPSlopeContactLaw(PotentialOfContactLaw): @staticmethod - def h_nu(uN, t): - raise NotImplementedError() + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + return -1.0 @staticmethod - def h_tau(uN, t): - raise NotImplementedError() + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + """ + electric charge flux - @staticmethod - def electric_charge_tangetial(u_tau): # potential - return 0 + var_nu == charge + """ + return 0.0 @staticmethod - def electric_charge_flux(charge): - return 0 * charge + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + """electric charge tangential""" + return 0 * np.linalg.norm(var_tau) @dataclass() @@ -41,7 +48,8 @@ class PQuasistaticSetup(PiezoelectricQuasistaticProblem): th_coef: ... = 4.5 ze_coef: ... = 10.5 time_step: ... = 0.01 - contact_law: ... = PPSlopeContactLaw + contact_law: ... = make_slope_contact_law(slope=1e1) + contact_law_2: ... = PPSlopeContactLaw piezoelectricity: ... = field( default_factory=lambda: np.array( [ @@ -66,10 +74,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: 0.0 <= x[0] <= 1.0 and 0.0 <= x[1] <= 1.0, dirichlet=lambda x: 1.0 <= x[0] <= 1.5 and 1.0 <= x[1] <= 1.5, @@ -93,8 +97,8 @@ def main(config: Config): mesh_descr = Barboteu2008MeshDescription(initial_position=None, max_element_perimeter=0.5) setup = PQuasistaticSetup(mesh_descr) runner = PiezoelectricTimeDependentSolver(setup, solving_method="global") - steps = 100 if not config.test else 10 - output = steps // 5 + steps = 100 if not config.test else 2 + output = steps // 5 if not config.test else 1 states = runner.solve( n_steps=steps, output_step=range(0, steps + 1, output), diff --git a/examples/example_piezoelectric_dynamic.py b/examples/example_piezoelectric_dynamic.py index 2605b114..086c81e9 100644 --- a/examples/example_piezoelectric_dynamic.py +++ b/examples/example_piezoelectric_dynamic.py @@ -2,37 +2,43 @@ Created at 21.08.2019 """ -from argparse import ArgumentParser from dataclasses import dataclass, field import numpy as np +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer from conmech.scenarios.problems import PiezoelectricDynamicProblem from conmech.simulations.problem_solver import PiezoelectricTimeDependentSolver from conmech.properties.mesh_description import Barboteu2008MeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law # TODO # 48 -class PPSlopeContactLaw(make_slope_contact_law(slope=1e1)): +class PPSlopeContactLaw(PotentialOfContactLaw): @staticmethod - def h_nu(uN, t): - raise NotImplementedError() + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + return -1.0 @staticmethod - def h_tau(uN, t): - raise NotImplementedError() + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + """ + electric charge flux - @staticmethod - def electric_charge_tangetial(u_tau): # potential - return 0 * 0.1 * 0.5 * np.linalg.norm(u_tau) ** 2 + var_nu == charge + """ + return 0.0 @staticmethod - def electric_charge_flux(charge): - return 0 * charge + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + """electric charge tangential""" + return 0 * 0.1 * 0.5 * np.linalg.norm(var_tau) ** 2 @dataclass() @@ -42,7 +48,8 @@ class PDynamicSetup(PiezoelectricDynamicProblem): th_coef: ... = 4.5 ze_coef: ... = 10.5 time_step: ... = 0.01 - contact_law: ... = PPSlopeContactLaw + contact_law: ... = make_slope_contact_law(slope=1e1) + contact_law_2: ... = PPSlopeContactLaw piezoelectricity: ... = field( default_factory=lambda: np.array( [ @@ -67,10 +74,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: 0.0 <= x[0] <= 1.0 and 0.0 <= x[1] <= 1.0, dirichlet=lambda x: 1.0 <= x[0] <= 1.5 and 1.0 <= x[1] <= 1.5, @@ -91,12 +94,14 @@ def main(config: Config): To see result of simulation you need to call from python `main(Config().init())`. """ - mesh_descr = Barboteu2008MeshDescription(initial_position=None, max_element_perimeter=0.25) + mesh_descr = Barboteu2008MeshDescription( + initial_position=None, max_element_perimeter=0.25 if not config.test else 0.5 + ) setup = PDynamicSetup(mesh_descr) runner = PiezoelectricTimeDependentSolver(setup, solving_method="global") - steps = 100 if not config.test else 10 - output = steps // 5 + steps = 100 if not config.test else 2 + output = steps // 5 if not config.test else 1 states = runner.solve( n_steps=steps, output_step=range(0, steps + 1, output), diff --git a/examples/example_poisson.py b/examples/example_poisson.py index 5f770109..9f71ede6 100644 --- a/examples/example_poisson.py +++ b/examples/example_poisson.py @@ -53,7 +53,11 @@ def main(config: Config): drawer.cmap = "plasma" drawer.field_name = "temperature" drawer.draw( - show=config.show, save=config.save, foundation=False, field_max=max_, field_min=min_ + show=config.show, + save=config.save, + foundation=False, + field_max=max_, + field_min=min_, ) diff --git a/examples/example_quasistatic.py b/examples/example_quasistatic.py index 55078955..449a7418 100644 --- a/examples/example_quasistatic.py +++ b/examples/example_quasistatic.py @@ -12,7 +12,7 @@ from conmech.simulations.problem_solver import TimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law @dataclass() @@ -32,10 +32,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu: float) -> float: - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) diff --git a/examples/example_static.py b/examples/example_static.py index 164ed7ac..a861cb7d 100644 --- a/examples/example_static.py +++ b/examples/example_static.py @@ -16,14 +16,11 @@ CubeMeshDescription, ) -from examples.p_slope_contact_law import make_slope_contact_law - -DIMENSION = 2 +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law @dataclass class StaticSetup(StaticDisplacementProblem): - dimension = DIMENSION mu_coef: ... = 4 la_coef: ... = 4 contact_law: ... = make_slope_contact_law(slope=1) @@ -36,25 +33,20 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return 0 * x - @staticmethod - def friction_bound(u_nu: float) -> float: - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) -def main(config: Config): +def main(config: Config, dimension=2): """ Entrypoint to example. To see result of simulation you need to call from python `main(Config().init())`. """ - # mesh_type = "cross" if DIMENSION == 2 else "meshzoo_cube_3d" - solving_method = "schur" if DIMENSION == 2 else "global" # TODO + solving_method = "schur" if dimension == 2 else "global" - if DIMENSION == 2: + if dimension == 2: mesh_descr = RectangleMeshDescription( initial_position=None, max_element_perimeter=0.5, scale=[2.5, 1] ) @@ -65,7 +57,7 @@ def main(config: Config): runner = StaticSolver(setup, solving_method) state = runner.solve(verbose=True, initial_displacement=setup.initial_displacement) - if DIMENSION == 2: + if dimension == 2: Drawer(state=state, config=config).draw(show=config.show, save=config.save) else: fig = plt.figure() @@ -85,4 +77,4 @@ def main(config: Config): if __name__ == "__main__": - main(Config().init()) + main(Config().init(), 2) diff --git a/examples/example_temperature_dynamic.py b/examples/example_temperature_dynamic.py index d76c0dc1..24c3574d 100644 --- a/examples/example_temperature_dynamic.py +++ b/examples/example_temperature_dynamic.py @@ -2,57 +2,41 @@ Created at 21.08.2019 """ -from argparse import ArgumentParser from dataclasses import dataclass, field import numpy as np +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.helpers.config import Config from conmech.mesh.boundaries_description import BoundariesDescription from conmech.plotting.drawer import Drawer from conmech.scenarios.problems import TemperatureDynamicProblem from conmech.simulations.problem_solver import TemperatureTimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law -class TPSlopeContactLaw(make_slope_contact_law(slope=1e1)): - # @staticmethod # TODO # 48 - # def g(t): - # return 10.7 + t * 0.02 - # return 0.5 + t * 0.01 - +class TPSlopeContactLaw(PotentialOfContactLaw): @staticmethod - def h_nu(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 100.0 * (uN - g_t) - return 0 + def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + """ + Direction of heat flux + """ + return -1.0 @staticmethod - def h_tau(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 10.0 * (uN - g_t) - return 0 - - # def jT(self, vTx, vTy): # TODO # 48 - # # return np.log(np.linalg.norm(vTx, vTy)+1) - # return np.linalg.norm(vTx, vTy) + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + """Temperature exchange""" + return 0.0 @staticmethod - def h_temp(u_tau): # potential # TODO # 48 - return 0.1 * 0.5 * np.linalg.norm(u_tau) ** 2 - - # TODO #96 : ContactLaw abstract class - - @staticmethod - def temp_exchange(temp): # potential # TODO # 48 - return 0 * temp - - @staticmethod - def h_temp(u_tau): # potential # TODO # 48 - return 0 * np.linalg.norm(u_tau) + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + """Friction generated temperature""" + return 0 * np.linalg.norm(var_tau) @dataclass() @@ -62,7 +46,8 @@ class TDynamicSetup(TemperatureDynamicProblem): th_coef: ... = 4 ze_coef: ... = 4 time_step: ... = 0.02 - contact_law: ... = TPSlopeContactLaw + contact_law: ... = make_slope_contact_law(slope=1e1) + contact_law_2: ... = TPSlopeContactLaw thermal_expansion: ... = field( default_factory=lambda: np.array([[0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5]]) ) @@ -86,10 +71,6 @@ def outer_forces(x, t=None): return np.array([-48.0 * (0.25 - (x[1] - 0.5) ** 2), 0]) return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription(contact=lambda x: x[1] == 0, dirichlet=lambda x: False) diff --git a/examples/example_temperature_dynamic_2.py b/examples/example_temperature_dynamic_2.py deleted file mode 100644 index 597001ef..00000000 --- a/examples/example_temperature_dynamic_2.py +++ /dev/null @@ -1,227 +0,0 @@ -""" -Created at 21.08.2019 -""" - -from dataclasses import dataclass - -import numpy as np -import pickle -import matplotlib.pyplot as plt - - -from conmech.helpers.config import Config -from conmech.mesh.boundaries_description import BoundariesDescription -from conmech.plotting.drawer import Drawer -from conmech.scenarios.problems import TemperatureDynamicProblem -from conmech.simulations.problem_solver import ( - TemperatureTimeDependentSolver as TDynamicProblemSolver, -) -from conmech.state.state import TemperatureState -from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law -import matplotlib.tri as tri - -# TODO #99 - - -def compute_error(ref: TemperatureState, sol: TemperatureState): - x = sol.body.mesh.nodes[:, 0] - y = sol.body.mesh.nodes[:, 1] - - soltri = tri.Triangulation(x, y, triangles=sol.body.mesh.elements) - u1hi = tri.LinearTriInterpolator(soltri, sol.velocity[:, 0]) - u2hi = tri.LinearTriInterpolator(soltri, sol.velocity[:, 1]) - thi = tri.LinearTriInterpolator(soltri, sol.temperature) - - total_abs_error = np.full_like(ref.temperature, fill_value=np.nan) - - for element in ref.body.mesh.elements: - x0 = ref.body.mesh.nodes[element[0]] - x1 = ref.body.mesh.nodes[element[1]] - x2 = ref.body.mesh.nodes[element[2]] - if total_abs_error[element[0]] != np.nan: - total_abs_error[element[0]] = ( # abs(ref.velocity[element[0], 0] - u1hi(*x0)) - # + abs(ref.velocity[element[0], 1] - u2hi(*x0)) - +abs(ref.temperature[element[0]] - thi(*x0)) - / ref.temperature[element[0]] - ) - if total_abs_error[element[1]] != np.nan: - total_abs_error[element[1]] = ( # abs(ref.velocity[element[1], 0] - u1hi(*x1)) - # + abs(ref.velocity[element[1], 1] - u2hi(*x1)) - +abs(ref.temperature[element[1]] - thi(*x1)) - / ref.temperature[element[1]] - ) - if total_abs_error[element[2]] != np.nan: - total_abs_error[element[2]] = ( # abs(ref.velocity[element[2], 0] - u1hi(*x2)) - # + abs(ref.velocity[element[2], 1] - u2hi(*x2)) - +abs(ref.temperature[element[2]] - thi(*x2)) - / ref.temperature[element[2]] - ) - - return total_abs_error - - -class TPSlopeContactLaw(make_slope_contact_law(slope=1e1)): - @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: - return 0.0 - return (0.5 * 1e3 * u_nu) * u_nu - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - if u_nu <= 0: - return 0 * v_nu - return (1e3 * u_nu) * v_nu - - @staticmethod - def potential_tangential_direction(u_tau: np.ndarray) -> float: - return np.log(np.sum(u_tau * u_tau) ** 0.5 + 1) - - @staticmethod - def h_nu(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 100.0 * (uN - g_t) - return 0 - - @staticmethod - def h_tau(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 10.0 * (uN - g_t) - return 0 - - @staticmethod - def temp_exchange(temp): # potential # TODO # 48 - return 0.1 * 0.5 * (temp - 0.27) ** 2 - - @staticmethod - def h_temp(u_tau): # potential # TODO # 48 - return 0.1 * 0.5 * np.linalg.norm(u_tau) ** 2 - - -@dataclass() -class TDynamicSetup(TemperatureDynamicProblem): - mu_coef: ... = 45 - la_coef: ... = 105 - th_coef: ... = 4.5 - ze_coef: ... = 10.5 - time_step: ... = 0.01 - contact_law: ... = TPSlopeContactLaw - thermal_expansion: ... = np.eye(3) * 0.5 - - thermal_conductivity: ... = np.eye(3) * 0.033 - - @staticmethod - def initial_temperature(x: np.ndarray) -> np.ndarray: - return np.asarray([0.25]) - - @staticmethod - def inner_forces(x, t=None): - return np.array([0.0, -1]) - - @staticmethod - def outer_forces(x, t=None): - if x[0] == 0: - return np.array([48.0 * (0.25 - (x[1] - 0.5) ** 2), 0]) - if x[0] == 1.5: - return np.array([-44.0 * (0.25 - (x[1] - 0.5) ** 2), 0]) - return np.array([0, 0]) - - @staticmethod - def friction_bound(u_nu): - return 1 - - boundaries: ... = BoundariesDescription(contact=lambda x: x[1] == 0) - - -# TODO: #99 -def main(steps, setup, config: Config): - simulate = True - output_step = (2**i for i in range(int(np.log2(steps)))) - - if setup is None: - mesh_descr = CrossMeshDescription( - initial_position=None, max_element_perimeter=0.25, scale=[1, 1] - ) - setup = TDynamicSetup(mesh_descr) - runner = TDynamicProblemSolver(setup, solving_method="schur") - - # for step, state in zip(output_step, runner.solve( - # n_steps=257, - # output_step=output_step, - # verbose=True, - # initial_displacement=setup.initial_displacement, - # initial_velocity=setup.initial_velocity, - # initial_temperature=setup.initial_temperature, - # )): - # with open(f'output/animation/k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}_t_{step}', - # 'wb') as output: - # pickle.dump(state, output) - - # fig, axes = plt.subplots(3, 2) - # for si, step in enumerate([0, 64, 256]): - # with open(f'output/animation/k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}_t_{step}', - # 'rb') as output: - # state6 = pickle.load(output) - # print(f"k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}") - # # Drawer(state=state6, config=config).draw( - # # temp_max=np.max(state6.temperature), temp_min=np.min(state6.temperature), draw_mesh=False, show=show, save=save - # # ) - # with open(f'output/animation/k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))-1}_t_{step}', - # 'rb') as output: - # state5 = pickle.load(output) - # print(f"k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}") - # drawer = Drawer(state=state5, config=config) - # drawer.node_size = 0 - # drawer.draw(#fig_axes=(fig, axes[si, 0]), - # temp_max=0.275, temp_min=0.25, draw_mesh=False, show=show, save=False, title=f"t={step / 512}" - # ) - # state_err = state6.copy() - # state_err.temperature = compute_error(state6, state5) - # drawer = Drawer(state=state_err, config=config, colormap="PuRd") - # drawer.node_size = 0 - # drawer.draw(#fig_axes=(fig, axes[si, 1]), - # temp_max=0.0005, temp_min=0, draw_mesh=False, show=show, - # save=save, title=f"t={step / 512}" - # ) - # plt.show() - - states = runner.solve( - n_steps=steps // 2 + 1, - output_step=(0, 4, 8, 16, 32, 64, 128, 256), - verbose=True, - initial_displacement=setup.initial_displacement, - initial_velocity=setup.initial_velocity, - initial_temperature=setup.initial_temperature, - ) - for state in states: - # with open(f'output/animation/k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}', - # 'wb') as output: - # pickle.dump(state, output) - - # with open(f'output/animation/k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}', - # 'rb') as output: - # state = pickle.load(output) - # print(f"k_{int(np.log2(steps))}_h_{int(np.log2(setup.elements_number[0]))}") - Drawer(state=state, config=config).draw( - field_max=np.max(state.temperature), - field_min=np.min(state.temperature), - show=config.show, - save=config.save, - ) - - -if __name__ == "__main__": - T = 1 - ks = [2**i for i in [2]] - hs = [2**i for i in [2]] - for h in hs: - for k in ks: - mesh_descr = CrossMeshDescription( - initial_position=None, max_element_perimeter=1 / h, scale=[1.5, 1] - ) - setup = TDynamicSetup(mesh_descr) - setup.time_step = T / k - main(setup=setup, steps=k, config=Config().init()) diff --git a/examples/p_slope_contact_law.py b/examples/p_slope_contact_law.py deleted file mode 100644 index 7b1e6e1b..00000000 --- a/examples/p_slope_contact_law.py +++ /dev/null @@ -1,60 +0,0 @@ -""" -Example contact law -""" - -from typing import Type - -import numpy as np -from conmech.scenarios.problems import ContactLaw - - -def make_slope_contact_law(slope: float) -> Type[ContactLaw]: - class PSlopeContactLaw(ContactLaw): - @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: - return 0.0 - return (0.5 * slope * u_nu) * u_nu - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - if u_nu <= 0: - return 0 * v_nu - return (slope * u_nu) * v_nu - - @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - regularization = 1 / np.sqrt(u_tau[0] * u_tau[0] + u_tau[1] * u_tau[1] + rho**2) - result = regularization * (u_tau[0] * v_tau[0] + u_tau[1] * v_tau[1]) - return result - - return PSlopeContactLaw - - -def make_const_contact_law(slope: float) -> Type[ContactLaw]: - class PSlopeContactLaw(ContactLaw): - @staticmethod - def potential_normal_direction(u_nu: float) -> float: - if u_nu <= 0: - return 0.0 - return slope * u_nu - - @staticmethod - def subderivative_normal_direction(u_nu: float, v_nu: float) -> float: - return 0 - - @staticmethod - def regularized_subderivative_tangential_direction( - u_tau: np.ndarray, v_tau: np.ndarray, rho=1e-7 - ) -> float: - """ - Coulomb regularization - """ - return 0 - - return PSlopeContactLaw diff --git a/examples/utils.py b/examples/utils.py index b1578453..04a09f6b 100644 --- a/examples/utils.py +++ b/examples/utils.py @@ -110,7 +110,13 @@ def viscoelastic_constitutive_law(displacement, velocity, setup, elements, nodes def elastic_relaxation_constitutive_law( - displacement: np.ndarray, absement: np.ndarray, setup, elements, nodes, time, **_kwargs + displacement: np.ndarray, + absement: np.ndarray, + setup, + elements, + nodes, + time, + **_kwargs, ): # TODO! grad_x = gradient(elements, nodes, displacement[:, 0]) grad_y = gradient(elements, nodes, displacement[:, 1]) diff --git a/requirements.txt b/requirements.txt index b46d9852..4f539170 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ jupyter==1.0.0 matplotlib==3.7.2 scipy==1.10.1 pytest==7.4.0 -meshzoo==0.9.11 +meshzoo @ git+https://github.com/KOS-UJ/meshzoo@main pygmsh==7.1.17 gmsh==4.11.1 tqdm==4.65.0 diff --git a/tests/test_conmech/regression/std_boundary.py b/tests/test_conmech/regression/std_boundary.py index 1d6404cc..1cf4d44f 100644 --- a/tests/test_conmech/regression/std_boundary.py +++ b/tests/test_conmech/regression/std_boundary.py @@ -1,6 +1,9 @@ import numpy as np -from conmech.mesh.boundaries_factory import get_boundary_surfaces, extract_unique_indices +from conmech.mesh.boundaries_factory import ( + get_boundary_surfaces, + extract_unique_indices, +) def standard_boundary_nodes(nodes, elements): diff --git a/tests/test_conmech/regression/test_Jureczka_and_Ochal_2019.py b/tests/test_conmech/regression/test_Jureczka_and_Ochal_2019.py index 9653535e..8b9fd96f 100644 --- a/tests/test_conmech/regression/test_Jureczka_and_Ochal_2019.py +++ b/tests/test_conmech/regression/test_Jureczka_and_Ochal_2019.py @@ -14,49 +14,25 @@ def solving_method(request): def test(solving_method): expected_displacement_vector = [ [0.0, 0.0], - [0.03522857, 0.02154723], - [0.06134879, 0.02852184], - [0.08393532, 0.03320384], - [0.10426168, 0.03618088], - [0.12257757, 0.03834405], - [0.13894498, 0.04001477], - [0.15337147, 0.04135569], - [0.1658617, 0.04243747], - [0.17642638, 0.04329012], - [0.18508084, 0.04392449], - [0.19184122, 0.04434147], - [0.1967216, 0.04453428], - [0.19973326, 0.04448642], - [0.20088703, 0.0441648], - [0.20020383, 0.04351448], - [0.19774689, 0.04237685], - [0.19891282, 0.05173695], - [0.19972586, 0.05984292], - [0.20010308, 0.06674484], - [0.19996122, 0.07244883], - [0.19923657, 0.07693874], - [0.19788669, 0.08018577], - [0.19589475, 0.08215244], - [0.19327284, 0.08279996], - [0.19239862, 0.07968588], - [0.18976624, 0.07607903], - [0.18540719, 0.07197493], - [0.17934638, 0.06737894], - [0.1716013, 0.0622999], - [0.16218933, 0.05673946], - [0.15113744, 0.05068978], - [0.13849374, 0.04413655], - [0.12434129, 0.03706901], - [0.10881389, 0.02950038], - [0.09211286, 0.02150421], - [0.0745211, 0.01327832], - [0.05640314, 0.00524877], - [0.03817266, -0.00165521], - [0.01990921, -0.00585183], - [0.0, 0.0], - [0.0, 0.0], - [0.0, 0.0], - [0.0, 0.0], + [0.06078453, 0.03250193], + [0.10404238, 0.03644206], + [0.13873971, 0.04036119], + [0.16570003, 0.04259187], + [0.18495207, 0.04404475], + [0.19660148, 0.0446417], + [0.20075845, 0.0442943], + [0.19765212, 0.04232086], + [0.19967536, 0.05988723], + [0.19991491, 0.07250985], + [0.19784972, 0.08020974], + [0.19323891, 0.08276982], + [0.1897699, 0.07613559], + [0.17937263, 0.06748543], + [0.16222708, 0.05688091], + [0.13856141, 0.04428154], + [0.10898344, 0.0296097], + [0.07489846, 0.01345543], + [0.03838516, -0.00183673], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0], @@ -64,12 +40,13 @@ def test(solving_method): ] mesh_descr = CrossMeshDescription( - initial_position=None, max_element_perimeter=0.125, scale=[2, 1] + initial_position=None, max_element_perimeter=0.25, scale=[2, 1] ) setup = StaticSetup(mesh_descr) runner = StaticProblem(setup, solving_method) result = runner.solve( - fixed_point_abs_tol=0.001, initial_displacement=setup.initial_displacement + initial_displacement=setup.initial_displacement, + method="Powell" if solving_method == "schur" else "BFGS", ) displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] @@ -79,7 +56,7 @@ def test(solving_method): np.set_printoptions(precision=8, suppress=True) print(repr(displacement[std_ids])) - precision = 2 if solving_method == "global optimization" else 3 + precision = 1 if solving_method != "global optimization" else 3 np.testing.assert_array_almost_equal( displacement[std_ids], expected_displacement_vector, decimal=precision ) diff --git a/tests/test_conmech/regression/test_density_update.py b/tests/test_conmech/regression/test_density_update.py index f6dcf4ab..80afb3c1 100644 --- a/tests/test_conmech/regression/test_density_update.py +++ b/tests/test_conmech/regression/test_density_update.py @@ -11,7 +11,7 @@ from conmech.scenarios.problems import StaticDisplacementProblem from conmech.simulations.problem_solver import NonHomogenousSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -51,10 +51,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -172,14 +168,10 @@ def inner_forces(x, t=None): @staticmethod def outer_forces(x, t=None): - return np.array([0, 0.1]) if x[1] < 0.2 and x[0] > 2.2 else np.array([0, 0]) - - @staticmethod - def friction_bound(u_nu: float) -> float: - return 0 + return np.array([0, 0.1]) if x[1] < 0.2 and x[0] >= 2.0 else np.array([0, 0]) boundaries: ... = BoundariesDescription( - contact=lambda x: x[1] == 0 and x[0] < 0.2, dirichlet=lambda x: x[0] == 0 + contact=lambda x: x[1] == 0 and x[0] <= 1.0, dirichlet=lambda x: x[0] == 0 ) mesh_descr = CrossMeshDescription( @@ -190,49 +182,49 @@ def friction_bound(u_nu: float) -> float: expected_displacement_vectors_2 = [ [ [0.0, 0.0], - [0.02599293, 0.021492], - [0.04747832, 0.06462195], - [0.12073424, 0.18288368], - [0.16852813, 0.35124948], - [0.19426259, 0.5945299], - [-0.00279276, 0.54908258], - [-0.17303802, 0.53807316], - [-0.16239225, 0.35466071], - [-0.12059053, 0.18422091], - [-0.04751902, 0.06465119], - [-0.02599838, 0.02148409], + [0.04618958, 0.03954446], + [0.08306192, 0.11665549], + [0.20096414, 0.32413195], + [0.26614292, 0.62621382], + [0.29242491, 0.94577544], + [0.00416508, 0.90460469], + [-0.26927487, 0.89272992], + [-0.25649393, 0.60371902], + [-0.1983749, 0.32574654], + [-0.08307566, 0.11683613], + [-0.0462062, 0.03954483], [0.0, 0.0], [0.0, 0.0], ], [ [0.0, 0.0], - [0.02602551, 0.0215782], - [0.04733487, 0.06410662], - [0.07700301, 0.1388405], - [0.09606589, 0.2350403], - [0.10636675, 0.36098208], - [-0.00111825, 0.34281676], - [-0.09787904, 0.33839926], - [-0.09361354, 0.23640489], - [-0.07694793, 0.13937786], - [-0.04736453, 0.06412507], - [-0.02602924, 0.02157238], + [0.04624452, 0.03968952], + [0.08281797, 0.11579863], + [0.13057317, 0.24686878], + [0.15654855, 0.41813314], + [0.16707307, 0.59601605], + [0.00165445, 0.57957066], + [-0.15783674, 0.57479782], + [-0.1527122, 0.40913347], + [-0.12955284, 0.2475101], + [-0.08281856, 0.1159435], + [-0.04625718, 0.03969205], [0.0, 0.0], [0.0, 0.0], ], [ [0.0, 0.0], - [0.05135096, 0.05502034], - [0.09153136, 0.18665172], - [0.11991624, 0.37669138], - [0.13705379, 0.60465541], - [0.14408795, 0.85793818], - [-0.10038042, 0.85018255], - [-0.31362968, 0.83921941], - [-0.30391213, 0.61065957], - [-0.26572229, 0.38741854], - [-0.19960097, 0.20262106], - [-0.1086065, 0.07793528], + [0.09165037, 0.10057917], + [0.16064146, 0.33569807], + [0.2058385, 0.66787976], + [0.22830296, 1.05861752], + [0.2350899, 1.46257643], + [-0.16208351, 1.45532423], + [-0.52254658, 1.44108317], + [-0.50920152, 1.06017022], + [-0.45275785, 0.68336219], + [-0.34627313, 0.3623118], + [-0.19077607, 0.14100217], [0.0, 0.0], [0.0, 0.0], ], diff --git a/tests/test_conmech/regression/test_dynamic.py b/tests/test_conmech/regression/test_dynamic.py index 15170e14..0f9555e5 100644 --- a/tests/test_conmech/regression/test_dynamic.py +++ b/tests/test_conmech/regression/test_dynamic.py @@ -11,7 +11,7 @@ from conmech.scenarios.problems import DynamicDisplacementProblem from conmech.simulations.problem_solver import TimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -42,10 +42,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -138,10 +134,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0.3, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0.0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) diff --git a/tests/test_conmech/regression/test_dynamic_membrane.py b/tests/test_conmech/regression/test_dynamic_membrane.py new file mode 100644 index 00000000..97f8f7b2 --- /dev/null +++ b/tests/test_conmech/regression/test_dynamic_membrane.py @@ -0,0 +1,259 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman-Szwarc +# +# 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 dataclasses import dataclass +from typing import Optional + +import numpy as np +import pytest + +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.scenarios.problems import WaveProblem +from conmech.simulations.problem_solver import WaveSolver +from conmech.properties.mesh_description import CrossMeshDescription + + +@pytest.fixture(params=["direct"]) +def solving_method(request): + return request.param + + +def generate_test_suits(): + test_suites = [] + + @dataclass() + class MembraneSetup(WaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([2.0]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1) + ) + + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=0.125, scale=[1, 1] + ) + setup_1 = MembraneSetup(mesh_descr) + + expected_displacement_vector_1 = [ + [-0.01001162, 0.0], + [-0.0161334, 0.0], + [-0.01735278, 0.0], + [-0.01761153, 0.0], + [-0.01761153, 0.0], + [-0.01735278, 0.0], + [-0.0161334, 0.0], + [-0.01001162, 0.0], + [-0.0161334, 0.0], + [-0.02990235, 0.0], + [-0.03301072, 0.0], + [-0.03370103, 0.0], + [-0.03370103, 0.0], + [-0.03301072, 0.0], + [-0.02990235, 0.0], + [-0.0161334, 0.0], + [-0.01735278, 0.0], + [-0.03301072, 0.0], + [-0.03699519, 0.0], + [-0.03793492, 0.0], + [-0.03793492, 0.0], + [-0.03699519, 0.0], + [-0.03301072, 0.0], + [-0.01735278, 0.0], + [-0.01761153, 0.0], + [-0.03370103, 0.0], + [-0.03793492, 0.0], + [-0.03896977, 0.0], + [-0.03896977, 0.0], + [-0.03793492, 0.0], + [-0.03370103, 0.0], + [-0.01761153, 0.0], + [-0.01761153, 0.0], + [-0.03370103, 0.0], + [-0.03793492, 0.0], + [-0.03896977, 0.0], + [-0.03896977, 0.0], + [-0.03793492, 0.0], + [-0.03370103, 0.0], + [-0.01761153, 0.0], + [-0.01735278, 0.0], + [-0.03301072, 0.0], + [-0.03699519, 0.0], + [-0.03793492, 0.0], + [-0.03793492, 0.0], + [-0.03699519, 0.0], + [-0.03301072, 0.0], + [-0.01735278, 0.0], + [-0.0161334, 0.0], + [-0.02990235, 0.0], + [-0.03301072, 0.0], + [-0.03370103, 0.0], + [-0.03370103, 0.0], + [-0.03301072, 0.0], + [-0.02990235, 0.0], + [-0.0161334, 0.0], + [-0.01001162, 0.0], + [-0.0161334, 0.0], + [-0.01735278, 0.0], + [-0.01761153, 0.0], + [-0.01761153, 0.0], + [-0.01735278, 0.0], + [-0.0161334, 0.0], + [-0.01001162, 0.0], + [-0.02333987, 0.0], + [-0.02789758, 0.0], + [-0.02889672, 0.0], + [-0.02907672, 0.0], + [-0.02889672, 0.0], + [-0.02789758, 0.0], + [-0.02333987, 0.0], + [-0.02789758, 0.0], + [-0.03484893, 0.0], + [-0.0365059, 0.0], + [-0.03681558, 0.0], + [-0.0365059, 0.0], + [-0.03484893, 0.0], + [-0.02789758, 0.0], + [-0.02889672, 0.0], + [-0.0365059, 0.0], + [-0.03844603, 0.0], + [-0.03882171, 0.0], + [-0.03844603, 0.0], + [-0.0365059, 0.0], + [-0.02889672, 0.0], + [-0.02907672, 0.0], + [-0.03681558, 0.0], + [-0.03882171, 0.0], + [-0.03921568, 0.0], + [-0.03882171, 0.0], + [-0.03681558, 0.0], + [-0.02907672, 0.0], + [-0.02889672, 0.0], + [-0.0365059, 0.0], + [-0.03844603, 0.0], + [-0.03882171, 0.0], + [-0.03844603, 0.0], + [-0.0365059, 0.0], + [-0.02889672, 0.0], + [-0.02789758, 0.0], + [-0.03484893, 0.0], + [-0.0365059, 0.0], + [-0.03681558, 0.0], + [-0.0365059, 0.0], + [-0.03484893, 0.0], + [-0.02789758, 0.0], + [-0.02333987, 0.0], + [-0.02789758, 0.0], + [-0.02889672, 0.0], + [-0.02907672, 0.0], + [-0.02889672, 0.0], + [-0.02789758, 0.0], + [-0.02333987, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + ] + + test_suites.append((setup_1, expected_displacement_vector_1)) + + # various changes + + @dataclass() + class MembraneSetup(WaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([0.0]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([1.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1) + ) + + setup_2 = MembraneSetup(mesh_descr) + + expected_displacement_vector_2 = np.zeros((145, 2)) + + test_suites.append((setup_2, expected_displacement_vector_2)) + + return test_suites + + +@pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) +def test_dynamic_membrane(solving_method, setup, expected_displacement_vector): + runner = WaveSolver(setup, solving_method) + results = runner.solve( + n_steps=2, + initial_displacement=setup.initial_displacement, + initial_velocity=setup.initial_velocity, + ) + + displacement = results[-1].body.mesh.nodes[:] - results[-1].displaced_nodes[:] + + # print result + np.set_printoptions(precision=8, suppress=True) + print(repr(displacement)) + + np.testing.assert_array_almost_equal(displacement, expected_displacement_vector, decimal=3) diff --git a/tests/test_conmech/regression/test_membrane_DNC.py b/tests/test_conmech/regression/test_membrane_DNC.py new file mode 100644 index 00000000..d73a11ef --- /dev/null +++ b/tests/test_conmech/regression/test_membrane_DNC.py @@ -0,0 +1,203 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman-Szwarc +# +# 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 dataclasses import dataclass +from typing import Optional + +import numpy as np +import pytest + +from conmech.dynamics.contact.damped_normal_compliance import make_damped_norm_compl +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.scenarios.problems import InteriorContactWaveProblem +from conmech.simulations.problem_solver import WaveSolver +from conmech.properties.mesh_description import CrossMeshDescription + + +@pytest.fixture(params=["global"]) +def solving_method(request): + return request.param + + +def generate_test_suits(): + test_suites = [] + + @dataclass() + class MembraneSetup(InteriorContactWaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + contact_law: ... = make_damped_norm_compl( + obstacle_level=0.01, kappa=1.0, beta=0.5, interior=True + )() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([100]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1) + ) + + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=1 / 4, scale=[1, 1] + ) + setup_1 = MembraneSetup(mesh_descr) + + expected_displacement_vector_1 = [ + [-1.0900725, -0.0], + [-1.6326846, -0.0], + [-1.6326838, -0.0], + [-1.0900737, -0.0], + [-1.6326846, -0.0], + [-2.7679646, -0.0], + [-2.767965, -0.0], + [-1.6326848, -0.0], + [-1.6326846, -0.0], + [-2.7679654, -0.0], + [-2.7679641, -0.0], + [-1.6326849, -0.0], + [-1.0900731, -0.0], + [-1.6326842, -0.0], + [-1.6326833, -0.0], + [-1.0900724, -0.0], + [-2.3944331, -0.0], + [-2.6196592, -0.0], + [-2.3944336, -0.0], + [-2.6196581, -0.0], + [-2.928575, -0.0], + [-2.6196588, -0.0], + [-2.3944323, -0.0], + [-2.6196583, -0.0], + [-2.3944328, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + ] + + test_suites.append((setup_1, expected_displacement_vector_1)) + + # various changes + + @dataclass() + class MembraneSetup(InteriorContactWaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + contact_law: ... = make_damped_norm_compl( + obstacle_level=0.01, kappa=1.0, beta=0.5, interior=True + )() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([-100]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0, 1) or x[1] in (0, 1) + ) + + setup_2 = MembraneSetup(mesh_descr) + + expected_displacement_vector_2 = [ + [1.0933539, -0.0], + [1.6375976, -0.0], + [1.6375984, -0.0], + [1.0933542, -0.0], + [1.6375979, -0.0], + [2.7762924, -0.0], + [2.7762933, -0.0], + [1.6375972, -0.0], + [1.6375985, -0.0], + [2.7762936, -0.0], + [2.7762937, -0.0], + [1.6375977, -0.0], + [1.0933532, -0.0], + [1.6375972, -0.0], + [1.6375986, -0.0], + [1.0933544, -0.0], + [2.401639, -0.0], + [2.6275417, -0.0], + [2.4016391, -0.0], + [2.6275408, -0.0], + [2.9373864, -0.0], + [2.627541, -0.0], + [2.4016392, -0.0], + [2.6275426, -0.0], + [2.4016395, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + [0.0, -0.0], + ] + + test_suites.append((setup_2, expected_displacement_vector_2)) + + return test_suites + + +@pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) +def test_membrane_DNC(solving_method, setup, expected_displacement_vector): + runner = WaveSolver(setup, solving_method) + results = runner.solve( + n_steps=2, + initial_displacement=setup.initial_displacement, + initial_velocity=setup.initial_velocity, + ) + + displacement = results[-1].body.mesh.nodes[:] - results[-1].displaced_nodes[:] + + # print result + np.set_printoptions(precision=7, suppress=True) + print(repr(displacement)) + + np.testing.assert_array_almost_equal(displacement, expected_displacement_vector, decimal=2) diff --git a/tests/test_conmech/regression/test_membrane_boundary_DNC.py b/tests/test_conmech/regression/test_membrane_boundary_DNC.py new file mode 100644 index 00000000..797717af --- /dev/null +++ b/tests/test_conmech/regression/test_membrane_boundary_DNC.py @@ -0,0 +1,249 @@ +# CONMECH @ Jagiellonian University in Kraków +# +# Copyright (C) 2023 Piotr Bartman-Szwarc +# +# 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 dataclasses import dataclass +from typing import Optional + +import numpy as np +import pytest + +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.scenarios.problems import ContactWaveProblem +from conmech.simulations.problem_solver import WaveSolver +from conmech.properties.mesh_description import CrossMeshDescription +from conmech.dynamics.contact.damped_normal_compliance import make_damped_norm_compl + + +@pytest.fixture(params=["schur", "global optimization"]) +def solving_method(request): + return request.param + + +def generate_test_suits(): + test_suites = [] + + @dataclass() + class MembraneSetup(ContactWaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + contact_law: ... = make_damped_norm_compl(0.01, kappa=1.0, beta=0.5)() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([100]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0,) or x[1] in (0, 1), + contact=lambda x: x[0] == 1, + ) + + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=1 / 6, scale=[1, 1] + ) + setup_1 = MembraneSetup(mesh_descr) + + expected_displacement_vector_1 = [ + [-1.50875561, 0.0], + [-1.80948135, 0.0], + [-1.85151694, 0.0], + [-1.80948058, 0.0], + [-1.50875435, 0.0], + [-1.07643442, 0.0], + [-1.21991845, 0.0], + [-1.24919285, 0.0], + [-1.23634646, 0.0], + [-1.13533086, 0.0], + [-0.65005518, 0.0], + [-2.08407485, 0.0], + [-2.4433707, 0.0], + [-2.5194924, 0.0], + [-2.48593163, 0.0], + [-2.23273285, 0.0], + [-1.13503756, 0.0], + [-2.29408448, 0.0], + [-2.74062776, 0.0], + [-2.83936081, 0.0], + [-2.79564395, 0.0], + [-2.48024186, 0.0], + [-1.2338268, 0.0], + [-2.29408387, 0.0], + [-2.74062736, 0.0], + [-2.83936113, 0.0], + [-2.79564248, 0.0], + [-2.48024111, 0.0], + [-1.2338265, 0.0], + [-2.08407363, 0.0], + [-2.44337049, 0.0], + [-2.51949141, 0.0], + [-2.48593195, 0.0], + [-2.23273261, 0.0], + [-1.13503913, 0.0], + [-1.07643343, 0.0], + [-1.21991743, 0.0], + [-1.24919236, 0.0], + [-1.23634806, 0.0], + [-1.13533196, 0.0], + [-0.65005674, 0.0], + [-1.99885568, 0.0], + [-2.12194468, 0.0], + [-2.13325662, 0.0], + [-2.06107777, 0.0], + [-1.66772572, 0.0], + [-2.54556707, 0.0], + [-2.74205763, 0.0], + [-2.76078498, 0.0], + [-2.64390479, 0.0], + [-2.05966075, 0.0], + [-2.63717937, 0.0], + [-2.85454336, 0.0], + [-2.87565865, 0.0], + [-2.74557381, 0.0], + [-2.12255092, 0.0], + [-2.54556523, 0.0], + [-2.74205761, 0.0], + [-2.76078459, 0.0], + [-2.64390367, 0.0], + [-2.05966124, 0.0], + [-1.99885497, 0.0], + [-2.12194321, 0.0], + [-2.13325616, 0.0], + [-2.0610787, 0.0], + [-1.66772549, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + ] + + test_suites.append((setup_1, expected_displacement_vector_1)) + + # various changes + + @dataclass() + class MembraneSetup(ContactWaveProblem): + time_step: ... = 0.1 + propagation: ... = 1.0 + contact_law: ... = make_damped_norm_compl(0.01, kappa=10.0, beta=0.5)() + + @staticmethod + def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray: + return np.array([-100]) + + @staticmethod + def outer_forces( + x: np.ndarray, v: Optional[np.ndarray] = None, t: Optional[float] = None + ) -> np.ndarray: + return np.array([0.0]) + + boundaries: ... = BoundariesDescription( + dirichlet=lambda x: x[0] in (0,) or x[1] in (0, 1), + contact=lambda x: x[0] == 1, + ) + + mesh_descr = CrossMeshDescription( + initial_position=None, max_element_perimeter=1 / 4, scale=[1, 1] + ) + setup_2 = MembraneSetup(mesh_descr) + + expected_displacement_vector_2 = [ + [2.64317675, 0.0], + [2.96235934, 0.0], + [2.64317508, 0.0], + [1.67732622, 0.0], + [1.67610052, 0.0], + [1.63886416, 0.0], + [1.09339295, 0.0], + [2.87365321, 0.0], + [2.87056248, 0.0], + [2.77947225, 0.0], + [1.63768291, 0.0], + [2.87365323, 0.0], + [2.87056375, 0.0], + [2.77947193, 0.0], + [1.63768443, 0.0], + [1.67732504, 0.0], + [1.67610041, 0.0], + [1.63886298, 0.0], + [1.09339442, 0.0], + [2.64291874, 0.0], + [2.63536374, 0.0], + [2.40190544, 0.0], + [2.96209659, 0.0], + [2.94988312, 0.0], + [2.62782706, 0.0], + [2.64291756, 0.0], + [2.63536304, 0.0], + [2.40190459, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + [0.0, 0.0], + ] + + test_suites.append((setup_2, expected_displacement_vector_2)) + + return test_suites + + +@pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) +def test_boundary_DNC(solving_method, setup, expected_displacement_vector): + runner = WaveSolver(setup, solving_method) + results = runner.solve( + n_steps=2, + initial_displacement=setup.initial_displacement, + initial_velocity=setup.initial_velocity, + ) + + displacement = results[-1].body.mesh.nodes[:] - results[-1].displaced_nodes[:] + + # print result + np.set_printoptions(precision=8, suppress=True) + print(repr(displacement)) + + precision = 2 if solving_method != "global optimization" else 3 + np.testing.assert_array_almost_equal(displacement, expected_displacement_vector, decimal=2) diff --git a/tests/test_conmech/regression/test_piezoelectric_quasistatic.py b/tests/test_conmech/regression/test_piezoelectric_quasistatic.py index b1376de2..17417b96 100644 --- a/tests/test_conmech/regression/test_piezoelectric_quasistatic.py +++ b/tests/test_conmech/regression/test_piezoelectric_quasistatic.py @@ -7,14 +7,14 @@ import numpy as np import pytest +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.mesh.boundaries_description import BoundariesDescription from conmech.scenarios.problems import ( PiezoelectricQuasistaticProblem, - PiezoelectricDynamicProblem, ) from conmech.simulations.problem_solver import PiezoelectricTimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -23,25 +23,28 @@ def solving_method(request): return request.param -def make_slope_contact_law_piezo(slope): - class PPSlopeContactLaw(make_slope_contact_law(slope=slope)): - @staticmethod - def h_nu(uN, t): - raise NotImplementedError() +class PPSlopeContactLaw(PotentialOfContactLaw): + @staticmethod + def tangential_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + return -1.0 - @staticmethod - def h_tau(uN, t): - raise NotImplementedError() + @staticmethod + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + """ + electric charge flux - @staticmethod - def electric_charge_tangetial(u_tau): # potential - return 0.1 * 0.5 * np.linalg.norm(u_tau) ** 2 - - @staticmethod - def electric_charge_flux(charge): - return 0 * charge + var_nu == charge + """ + return 0.0 - return PPSlopeContactLaw + @staticmethod + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + """electric charge tangential""" + return 0.1 * 0.5 * np.linalg.norm(var_tau) ** 2 def generate_test_suits(): @@ -56,7 +59,8 @@ class QuasistaticSetup_1(PiezoelectricQuasistaticProblem): th_coef: ... = 4 ze_coef: ... = 4 time_step: ... = 0.1 - contact_law: ... = make_slope_contact_law_piezo(1e1) + contact_law: ... = make_slope_contact_law(1e1) + contact_law_2: ... = PPSlopeContactLaw() piezoelectricity: ... = field( default_factory=lambda: np.array( [ @@ -77,10 +81,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -93,37 +93,37 @@ def friction_bound(u_nu): expected_displacement_vector_m02_m02 = np.asarray( [ [0.0, 0.0], - [0.02739129, 0.01551078], - [0.04193194, 0.0344827], - [0.04966183, 0.05129362], - [0.05372157, 0.06523827], - [0.05497221, 0.07868415], - [0.03990124, 0.08298064], - [0.02610119, 0.08494335], - [0.02444053, 0.07036843], - [0.01887287, 0.05538343], - [0.0105034, 0.03823251], - [0.00236344, 0.01854232], + [0.01284572, 0.00727844], + [0.01998403, 0.01155948], + [0.02406325, 0.01469925], + [0.0261574, 0.01657194], + [0.02644623, 0.01809549], + [0.02433448, 0.02099327], + [0.02235926, 0.02215269], + [0.02133084, 0.01972781], + [0.01788636, 0.01699742], + [0.01220893, 0.01293959], + [0.00540516, 0.00662938], [0.0, 0.0], [0.0, 0.0], ] ) expected_temperature_vector_m02_m02 = np.asarray( [ - -0.17213276, - -0.08602799, - -0.0221841, - 0.01236338, - 0.03289128, - 0.04208222, - 0.01226481, - 0.00902652, - 0.02242403, - 0.04431714, - 0.05231734, - 0.03319498, - -0.0243407, - -0.06355553, + -0.09317176, + -0.05082419, + -0.01792647, + 0.00284861, + 0.01673696, + 0.02415654, + -0.00054601, + -0.00339339, + 0.00680334, + 0.02489092, + 0.04042651, + 0.04127354, + 0.01658142, + -0.02771255, ] ) @@ -138,7 +138,7 @@ def friction_bound(u_nu): # p = 0 and opposite forces setup_0_02_p_0 = QuasistaticSetup_1(mesh_descr) - setup_0_02_p_0.contact_law = make_slope_contact_law_piezo(0) + setup_0_02_p_0.contact_law = make_slope_contact_law(0) def inner_forces(x, time=None): return np.array([0, 0.2]) @@ -148,37 +148,37 @@ def inner_forces(x, time=None): expected_displacement_vector_0_02_p_0 = np.asarray( [ [0.0, 0.0], - [-0.07690522, -0.06648437], - [-0.12262936, -0.19051381], - [-0.1462746, -0.34206085], - [-0.15508485, -0.50302121], - [-0.15648751, -0.66223258], - [0.00000041, -0.66244246], - [0.15648834, -0.66223257], - [0.15508565, -0.5030212], - [0.14627531, -0.34206083], - [0.12262989, -0.19051379], - [0.07690553, -0.06648434], + [-0.05119819, -0.04903926], + [-0.08203236, -0.13236709], + [-0.09784363, -0.23465962], + [-0.10375805, -0.34281109], + [-0.10470465, -0.44949937], + [0.00000025, -0.44956951], + [0.10470515, -0.44949937], + [0.10375853, -0.34281111], + [0.09784406, -0.23465962], + [0.08203271, -0.13236709], + [0.05119838, -0.04903929], [0.0, 0.0], [0.0, 0.0], ] ) expected_temperature_vector_0_02_p_0 = np.asarray( [ - 0.43120988, - 0.2234753, - 0.05332727, - -0.0674368, - -0.1399693, - -0.1641408, - -0.16740288, - -0.16414068, - -0.13996898, - -0.06743615, - 0.05332823, - 0.22347662, - 0.43121148, - 0.21898649, + 0.32217749, + 0.17962827, + 0.04606892, + -0.04874313, + -0.10557165, + -0.12452466, + -0.12761607, + -0.12452451, + -0.10557137, + -0.04874258, + 0.04606974, + 0.17962935, + 0.32217874, + 0.16140012, ] ) @@ -193,7 +193,7 @@ def inner_forces(x, time=None): # p = 0 setup_0_m02_p_0 = QuasistaticSetup_1(mesh_descr) - setup_0_m02_p_0.contact_law = make_slope_contact_law_piezo(0) + setup_0_m02_p_0.contact_law = make_slope_contact_law(0) def inner_forces(x, time=None): return np.array([0, -0.2]) @@ -218,7 +218,8 @@ class QuasistaticSetup_2(PiezoelectricQuasistaticProblem): th_coef: ... = 2.11 ze_coef: ... = 4.99 time_step: ... = 0.1 - contact_law: ... = make_slope_contact_law_piezo(2.71) + contact_law: ... = make_slope_contact_law(2.71) + contact_law_2: ... = PPSlopeContactLaw piezoelectricity: ... = field( default_factory=lambda: np.array( [ @@ -239,10 +240,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0.3, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0.0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -254,17 +251,17 @@ def friction_bound(u_nu): expected_displacement_vector_var = np.asarray( [ [0.0, 0.0], - [0.02454485, 0.05521995], - [0.03263034, 0.15296714], - [0.0326597, 0.28043961], - [0.02912746, 0.42943956], - [0.02194101, 0.58807768], - [-0.13470771, 0.59642012], - [-0.31508619, 0.59818791], - [-0.29559302, 0.43025764], - [-0.25711724, 0.27975335], - [-0.20030387, 0.14778518], - [-0.11928462, 0.04488547], + [-0.00767607, 0.01454776], + [-0.0191455, 0.04235841], + [-0.0297773, 0.07752592], + [-0.03753689, 0.12139559], + [-0.0456317, 0.17298394], + [-0.09733228, 0.18388321], + [-0.16801203, 0.18878844], + [-0.15417002, 0.1300507], + [-0.13068551, 0.08537821], + [-0.10023372, 0.04869301], + [-0.05979862, 0.02043574], [0.0, 0.0], [0.0, 0.0], ] @@ -272,20 +269,20 @@ def friction_bound(u_nu): expected_temperature_vector_var = np.asarray( [ - -0.00912808, - 0.09603883, - 0.14079504, - 0.14273884, - 0.13046765, - 0.15567758, - -0.02148364, - -0.17076294, - -0.1409075, - -0.14751809, - -0.17838383, - -0.2639388, - -0.39721732, - -0.00154481, + 0.18017535, + 0.20951466, + 0.19322749, + 0.13329264, + 0.08185755, + 0.09077227, + -0.10696354, + -0.26170111, + -0.21741867, + -0.17859079, + -0.14847378, + -0.1462184, + -0.17605883, + 0.08383385, ] ) @@ -308,7 +305,7 @@ def test_piezoelectric_time_dependent_solver( ): runner = PiezoelectricTimeDependentSolver(setup, solving_method) results = runner.solve( - n_steps=32, + n_steps=8, initial_displacement=setup.initial_displacement, initial_velocity=setup.initial_velocity, initial_electric_potential=setup.initial_electric_potential, diff --git a/tests/test_conmech/regression/test_poisson.py b/tests/test_conmech/regression/test_poisson.py index 61154f9c..911a5e9e 100644 --- a/tests/test_conmech/regression/test_poisson.py +++ b/tests/test_conmech/regression/test_poisson.py @@ -101,7 +101,9 @@ def outer_temperature(x, t=None): ) mesh_descr = CrossMeshDescription( - initial_position=None, max_element_perimeter=1.37 / 3, scale=[1.37 * 5 / 3, 1.37] + initial_position=None, + max_element_perimeter=1.37 / 3, + scale=[1.37 * 5 / 3, 1.37], ) setup_2 = StaticSetup(mesh_descr) expected_temperature_vector_2 = [ diff --git a/tests/test_conmech/regression/test_quasistatic.py b/tests/test_conmech/regression/test_quasistatic.py index d371e641..6908fc5f 100644 --- a/tests/test_conmech/regression/test_quasistatic.py +++ b/tests/test_conmech/regression/test_quasistatic.py @@ -11,7 +11,7 @@ from conmech.scenarios.problems import QuasistaticDisplacementProblem from conmech.simulations.problem_solver import TimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -42,10 +42,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -138,10 +134,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0.3, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0.0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) diff --git a/tests/test_conmech/regression/test_static.py b/tests/test_conmech/regression/test_static.py index 538c7a4e..2efe6f60 100644 --- a/tests/test_conmech/regression/test_static.py +++ b/tests/test_conmech/regression/test_static.py @@ -11,7 +11,7 @@ from conmech.scenarios.problems import StaticDisplacementProblem from conmech.simulations.problem_solver import StaticSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -39,10 +39,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -132,10 +128,6 @@ def inner_forces(x, t=None): def outer_forces(x, t=None): return np.array([0.3, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0.0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) diff --git a/tests/test_conmech/regression/test_static_3d.py b/tests/test_conmech/regression/test_static_3d.py new file mode 100644 index 00000000..69b82909 --- /dev/null +++ b/tests/test_conmech/regression/test_static_3d.py @@ -0,0 +1,181 @@ +""" +Created at 21.08.2019 +""" + +from dataclasses import dataclass + +import numpy as np +import pytest + +from conmech.mesh.boundaries_description import BoundariesDescription +from conmech.scenarios.problems import StaticDisplacementProblem +from conmech.simulations.problem_solver import StaticSolver +from conmech.properties.mesh_description import ( + CrossMeshDescription, + CubeMeshDescription, +) +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law +from tests.test_conmech.regression.std_boundary import standard_boundary_nodes + + +@pytest.fixture(params=["global optimization"]) +def solving_method(request): + return request.param + + +def generate_test_suits(): + test_suites = [] + + # Simple example + + @dataclass() + class StaticSetup(StaticDisplacementProblem): + mu_coef: ... = 4 + la_coef: ... = 4 + contact_law: ... = make_slope_contact_law(slope=1) + + @staticmethod + def inner_forces(x, t=None): + return -0.2 * x + + @staticmethod + def outer_forces(x, t=None): + return 0 * x + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + mesh_descr = CubeMeshDescription(initial_position=None) + setup_m02_m02 = StaticSetup(mesh_descr) + + expected_displacement_vector_m02_m02 = [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.00886805, 0.00759973, 0.01309122], + [0.01251517, 0.01304646, 0.02180133], + [0.0047388, 0.01377955, 0.02134457], + [-0.00179084, 0.01481824, 0.02136016], + [0.00211174, 0.01387662, 0.02174084], + [0.00165348, 0.00748375, 0.01153221], + [-0.00196903, 0.00853421, 0.01211016], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [-0.00548426, 0.00962111, 0.01202977], + [-0.0060029, 0.01529378, 0.02094435], + [0.00087541, 0.01483005, 0.02075665], + [-0.00010425, 0.00850936, 0.01158434], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.00577241, 0.00703007, 0.01156364], + [0.00711386, 0.0132022, 0.02045272], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.01268809, 0.0081164, 0.0137591], + [0.01424769, 0.01004692, 0.02280081], + [0.00928219, 0.01244232, 0.02200091], + [0.00655036, 0.00806369, 0.01241398], + ] + + test_suites.append((setup_m02_m02, expected_displacement_vector_m02_m02)) + + # p = 0 and opposite forces + + setup_0_02_p_0 = StaticSetup(mesh_descr) + setup_0_02_p_0.contact_law = make_slope_contact_law(slope=0) + + def inner_forces(x, t=None): + result = 0.2 * x + result[0] = 0.0 + return result + + setup_0_02_p_0.inner_forces = inner_forces + + expected_displacement_vector_0_02_p_0 = [ + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [-0.00521917, -0.01197369, -0.01232288], + [-0.00671187, -0.02140411, -0.02124774], + [0.00035836, -0.02143418, -0.02143418], + [0.0075367, -0.02140411, -0.02218005], + [0.00039165, -0.02055906, -0.02130339], + [0.00016529, -0.01087314, -0.01172428], + [0.00563059, -0.01197369, -0.01291526], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.01132733, -0.01349373, -0.01349373], + [0.01346888, -0.02202116, -0.02202116], + [0.0075367, -0.02218005, -0.02140411], + [0.00563059, -0.01291526, -0.01197369], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [0.00016529, -0.01172428, -0.01087314], + [0.00039165, -0.02130339, -0.02055906], + [0.0, 0.0, 0.0], + [0.0, 0.0, 0.0], + [-0.01099676, -0.01264259, -0.01264259], + [-0.01268559, -0.02127684, -0.02127684], + [-0.00671187, -0.02124774, -0.02140411], + [-0.00521917, -0.01232288, -0.01197369], + ] + + test_suites.append((setup_0_02_p_0, expected_displacement_vector_0_02_p_0)) + + # p = 0 + + setup_0_m02_p_0 = StaticSetup(mesh_descr) + setup_0_m02_p_0.contact_law = make_slope_contact_law(slope=0) + + def inner_forces(x, t=None): + result = -0.2 * x + result[0] = 0.0 + return result + + setup_0_m02_p_0.inner_forces = inner_forces + + expected_displacement_vector_0_m02_p_0 = [ + [-v[0], -v[1], -v[2]] for v in expected_displacement_vector_0_02_p_0 + ] + + test_suites.append((setup_0_m02_p_0, expected_displacement_vector_0_m02_p_0)) + + # various changes + + @dataclass() + class StaticSetup(StaticDisplacementProblem): + mu_coef: ... = 4.58 + la_coef: ... = 3.33 + contact_law: ... = make_slope_contact_law(slope=2.71) + + @staticmethod + def inner_forces(x, t=None): + return np.array([0, -0.2]) + + @staticmethod + def outer_forces(x, t=None): + return np.array([0.3, 0.0]) + + boundaries: ... = BoundariesDescription( + contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 + ) + + return test_suites + + +@pytest.mark.parametrize("setup, expected_displacement_vector", generate_test_suits()) +def test_static_solver(solving_method, setup, expected_displacement_vector): + runner = StaticSolver(setup, solving_method) + result = runner.solve(initial_displacement=setup.initial_displacement) + + displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:] + std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements) + + # print result + np.set_printoptions(precision=8, suppress=True) + print(repr(displacement[std_ids])) + + np.testing.assert_array_almost_equal( + displacement[std_ids], expected_displacement_vector, decimal=3 + ) diff --git a/tests/test_conmech/regression/test_temperature_dynamic.py b/tests/test_conmech/regression/test_temperature_dynamic.py index 85e6ee37..a7a4f014 100644 --- a/tests/test_conmech/regression/test_temperature_dynamic.py +++ b/tests/test_conmech/regression/test_temperature_dynamic.py @@ -7,11 +7,12 @@ import numpy as np import pytest +from conmech.dynamics.contact.contact_law import PotentialOfContactLaw from conmech.mesh.boundaries_description import BoundariesDescription from conmech.scenarios.problems import TemperatureDynamicProblem from conmech.simulations.problem_solver import TemperatureTimeDependentSolver from conmech.properties.mesh_description import CrossMeshDescription -from examples.p_slope_contact_law import make_slope_contact_law +from conmech.dynamics.contact.relu_slope_contact_law import make_slope_contact_law from tests.test_conmech.regression.std_boundary import standard_boundary_nodes @@ -20,37 +21,27 @@ def solving_method(request): return request.param -def make_slope_contact_law_temp(slope): - class TPSlopeContactLaw(make_slope_contact_law(slope=slope)): - @staticmethod - def h_nu(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 100.0 * (uN - g_t) - return 0 - - @staticmethod - def h_tau(uN, t): - g_t = 10.7 + t * 0.02 - if uN > g_t: - return 10.0 * (uN - g_t) - return 0 +class TPSlopeContactLaw(PotentialOfContactLaw): + @staticmethod + def normal_bound(var_nu: float, static_displacement_nu: float, dt: float) -> float: + """ + Direction of heat flux + """ + return -1.0 - @staticmethod - def h_temp(vT): - return 0.1 * np.linalg.norm(vT) + @staticmethod + def potential_normal_direction( + var_nu: float, static_displacement_nu: float, dt: float + ) -> float: + """Temperature exchange""" + return 0.0 - # TODO #96: ContactLaw abstract class - - @staticmethod - def temp_exchange(temp): # potential # TODO # 48 - return 0 * temp - - @staticmethod - def h_temp(u_tau): # potential # TODO # 48 - return 0 * np.linalg.norm(u_tau) - - return TPSlopeContactLaw + @staticmethod + def potential_tangential_direction( + var_tau: float, static_displacement_tau: float, dt: float + ) -> float: + """Friction generated temperature""" + return 0.0 def generate_test_suits(): @@ -65,7 +56,8 @@ class DynamicSetup(TemperatureDynamicProblem): th_coef: ... = 4 ze_coef: ... = 4 time_step: ... = 0.1 - contact_law: ... = make_slope_contact_law_temp(1e1) + contact_law: ... = make_slope_contact_law(1e1) + contact_law_2: ... = TPSlopeContactLaw() thermal_expansion: ... = field( default_factory=lambda: np.array([[0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5]]) ) @@ -81,10 +73,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0, 0]) - @staticmethod - def friction_bound(u_nu): - return 0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -97,17 +85,17 @@ def friction_bound(u_nu): expected_displacement_vector_m02_m02 = np.asarray( [ [0.0, 0.0], - [0.03142428, 0.02691658], - [0.04839312, 0.04238928], - [0.05698615, 0.05458973], - [0.06022096, 0.06307318], - [0.0591116, 0.07018657], - [0.05128066, 0.07756148], - [0.04143363, 0.07994437], - [0.03895115, 0.06888961], - [0.03117291, 0.05560756], - [0.01913418, 0.03800219], - [0.00700475, 0.01647643], + [0.00570683, 0.00431404], + [0.00857057, 0.00526714], + [0.01005346, 0.00586446], + [0.01058508, 0.00597765], + [0.01028078, 0.00577634], + [0.01051063, 0.0074202], + [0.01045592, 0.00795605], + [0.0099967, 0.00760991], + [0.00853351, 0.00684343], + [0.00600395, 0.00531373], + [0.00282262, 0.00263098], [0.0, 0.0], [0.0, 0.0], ] @@ -115,17 +103,17 @@ def friction_bound(u_nu): expected_temperature_vector_m02_m02 = np.asarray( [ 0.0, - 0.00735363, - 0.00913178, - 0.00828201, - 0.00694812, - 0.00636948, - 0.00607767, - 0.00578943, - 0.00648897, - 0.00789741, - 0.00825697, - 0.00542543, + 0.00407167, + 0.00284784, + 0.00232436, + 0.00175827, + 0.0014074, + 0.00094456, + 0.00049082, + 0.00088726, + 0.00154845, + 0.00210648, + 0.00257013, 0.0, 0.0, ] @@ -142,7 +130,7 @@ def friction_bound(u_nu): # p = 0 and opposite forces setup_0_02_p_0 = DynamicSetup(mesh_descr) - setup_0_02_p_0.contact_law = make_slope_contact_law_temp(0) + setup_0_02_p_0.contact_law = make_slope_contact_law(0) def inner_forces(x, time=None): return np.array([0, 0.2]) @@ -152,17 +140,17 @@ def inner_forces(x, time=None): expected_displacement_vector_0_02_p_0 = np.asarray( [ [0.0, 0.0], - [-0.09962716, -0.10791196], - [-0.16199068, -0.27160136], - [-0.19306197, -0.47359855], - [-0.20451057, -0.68641614], - [-0.20646559, -0.89647558], - [0.00000059, -0.89632143], - [0.20646675, -0.89647557], - [0.20451168, -0.68641617], - [0.19306293, -0.47359861], - [0.16199141, -0.27160145], - [0.09962757, -0.10791209], + [-0.00357335, -0.00721697], + [-0.00400228, -0.01306104], + [-0.00328983, -0.01728046], + [-0.00266257, -0.02009906], + [-0.00246619, -0.02248587], + [0.00000012, -0.02249418], + [0.00246643, -0.02248586], + [0.00266281, -0.02009906], + [0.00329004, -0.01728048], + [0.00400244, -0.01306106], + [0.00357344, -0.007217], [0.0, 0.0], [0.0, 0.0], ] @@ -170,17 +158,17 @@ def inner_forces(x, time=None): expected_temperature_vector_0_02_p_0 = np.asarray( [ 0.0, - -0.01368016, - -0.01005095, - -0.00546105, - -0.00217319, - -0.00098241, - -0.00000016, - 0.00098209, - 0.00217289, - 0.0054608, - 0.01005075, - 0.01368001, + -0.00125683, + 0.00007371, + 0.00030325, + 0.00019259, + 0.00008153, + -0.00000002, + -0.00008157, + -0.00019262, + -0.00030328, + -0.00007372, + 0.00125684, 0.0, 0.0, ] @@ -197,7 +185,7 @@ def inner_forces(x, time=None): # p = 0 setup_0_m02_p_0 = DynamicSetup(mesh_descr) - setup_0_m02_p_0.contact_law = make_slope_contact_law_temp(0) + setup_0_m02_p_0.contact_law = make_slope_contact_law(0) def inner_forces(x, time=None): return np.array([0, -0.2]) @@ -222,7 +210,8 @@ class DynamicSetup(TemperatureDynamicProblem): th_coef: ... = 2.11 ze_coef: ... = 4.99 time_step: ... = 0.1 - contact_law: ... = make_slope_contact_law_temp(2.71) + contact_law: ... = make_slope_contact_law(2.71) + contact_law_2: ... = TPSlopeContactLaw thermal_expansion: ... = field( default_factory=lambda: np.array([[0.5, 0.0, 0.0], [0.0, 0.5, 0.0], [0.0, 0.0, 0.5]]) ) @@ -238,10 +227,6 @@ def inner_forces(x, time=None): def outer_forces(x, time=None): return np.array([0.3, 0.0]) - @staticmethod - def friction_bound(u_nu): - return 0.0 - boundaries: ... = BoundariesDescription( contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0 ) @@ -253,17 +238,17 @@ def friction_bound(u_nu): expected_displacement_vector_var = np.asarray( [ [0.0, 0.0], - [0.00826122, 0.0421944], - [0.00410242, 0.13421482], - [-0.00990836, 0.25783836], - [-0.02692933, 0.40249172], - [-0.04900262, 0.56156278], - [-0.20968181, 0.59213072], - [-0.39831803, 0.6111903], - [-0.36578831, 0.43996181], - [-0.31149975, 0.29375818], - [-0.23806523, 0.16985506], - [-0.13971534, 0.07516066], + [-0.00224574, 0.00303146], + [-0.00635254, 0.00724873], + [-0.01133355, 0.01098753], + [-0.01611561, 0.01621815], + [-0.02404372, 0.02594285], + [-0.03629391, 0.03524607], + [-0.06458808, 0.04365407], + [-0.05360272, 0.02541212], + [-0.04325996, 0.01882516], + [-0.0340562, 0.01503337], + [-0.0222147, 0.01126095], [0.0, 0.0], [0.0, 0.0], ] @@ -272,17 +257,17 @@ def friction_bound(u_nu): expected_temperature_vector_var = np.asarray( [ 0.0, - -0.00870148, - -0.01319496, - -0.00947956, - -0.00501162, - -0.00299865, - -0.00606658, - -0.00913445, - -0.01288869, - -0.02153761, - -0.03067854, - -0.03005361, + -0.00160508, + -0.00133705, + -0.00092485, + -0.00083582, + -0.00073482, + -0.00167299, + -0.00278935, + -0.00325494, + -0.00339006, + -0.00344784, + -0.00747268, 0.0, 0.0, ] @@ -304,7 +289,7 @@ def test_temperature_time_dependent_solver( ): runner = TemperatureTimeDependentSolver(setup, solving_method) result_generator = runner.solve( - n_steps=32, + n_steps=4, initial_displacement=setup.initial_displacement, initial_velocity=setup.initial_velocity, initial_temperature=setup.initial_temperature, diff --git a/tests/test_conmech/smoke/test_examples.py b/tests/test_conmech/smoke/test_examples.py index 852675ed..e88c3a49 100644 --- a/tests/test_conmech/smoke/test_examples.py +++ b/tests/test_conmech/smoke/test_examples.py @@ -13,9 +13,13 @@ from examples.example_dynamic import main as dynamic from examples.example_piezoelectric_dynamic import main as dynamic_piezo from examples.example_temperature_dynamic import main as dynamic_temp +from examples.example_dynamic_membrane import main as dynamic_membrane from examples.Jureczka_and_Ochal_2019 import main as Jureczka_and_Ochal_2019 from examples.Jureczka_Ochal_Bartman_2023 import main as Jureczka_Ochal_Bartman_2023 from examples.Sofonea_Ochal_Bartman_2023 import main as Sofonea_Ochal_Bartman_2023 +from examples.BOST_2024 import main as BOST_2024 +from examples.BOSK_2024 import main as BOSK_2024 +from examples.example_Tarzia_problem import main as Tarzia_problem from examples.example_poisson import main as poisson from examples.example_nonhomogenous_density import main as nonhomogenous_density from examples.example_imported_mesh_static import main as imported_mesh @@ -26,6 +30,7 @@ test_suits = { "poisson": lambda: poisson(Config(**default_args).init()), + "dynamic_membrane": lambda: dynamic_membrane(Config(**default_args).init()), "static": lambda: static(Config(**default_args).init()), "quasistatic": lambda: quasistatic(Config(**default_args).init()), "quasi_piezo": lambda: quasi_piezo(Config(**default_args).init()), @@ -39,6 +44,11 @@ "Sofonea_Ochal_Bartman_2023": lambda: Sofonea_Ochal_Bartman_2023( Config(outputs_path="./output/SOB2023", **default_args).init() ), + "BOST_2024": lambda: BOST_2024(Config(outputs_path="./output/BOST2024", **default_args).init()), + "BOSK_2024": lambda: BOSK_2024(Config(outputs_path="./output/BOSK2024", **default_args).init()), + "Tarzia_problem": lambda: Tarzia_problem( + Config(outputs_path="./output/BOST2024", **default_args).init() + ), "nonhomogenous_density": lambda: nonhomogenous_density(Config(**default_args).init()), "imported_mesh": lambda: imported_mesh( Config(**default_args).init(), diff --git a/tests/test_conmech/unit_test/test_matrices.py b/tests/test_conmech/unit_test/test_matrices.py index 59644414..e455a440 100644 --- a/tests/test_conmech/unit_test/test_matrices.py +++ b/tests/test_conmech/unit_test/test_matrices.py @@ -10,7 +10,10 @@ from conmech.simulations.problem_solver import Body from conmech.mesh.mesh import Mesh from conmech.mesh import mesh_builders -from conmech.properties.mesh_description import RectangleMeshDescription, CubeMeshDescription +from conmech.properties.mesh_description import ( + RectangleMeshDescription, + CubeMeshDescription, +) def test_matrices_2d_integrals(): @@ -20,7 +23,9 @@ def test_matrices_2d_integrals(): area = scale_x * scale_y initial_nodes, elements = mesh_builders.build_mesh( mesh_descr=RectangleMeshDescription( - initial_position=None, max_element_perimeter=(scale_x / 3), scale=[scale_x, scale_y] + initial_position=None, + max_element_perimeter=(scale_x / 3), + scale=[scale_x, scale_y], ) ) @@ -81,7 +86,9 @@ def test_local_stiff_mats_assembly(): scale_y = 3 initial_nodes, elements = mesh_builders.build_mesh( mesh_descr=RectangleMeshDescription( - initial_position=None, max_element_perimeter=(scale_x / 3), scale=[scale_x, scale_y] + initial_position=None, + max_element_perimeter=(scale_x / 3), + scale=[scale_x, scale_y], ) ) edges_features_matrix, _, local_stiff_mats = sut_2d(elements=elements, nodes=initial_nodes) diff --git a/tests/test_conmech/unit_test/test_mesh.py b/tests/test_conmech/unit_test/test_mesh.py index 8aa4ecb9..820d1dfa 100644 --- a/tests/test_conmech/unit_test/test_mesh.py +++ b/tests/test_conmech/unit_test/test_mesh.py @@ -7,7 +7,9 @@ from conmech.mesh.boundaries_description import BoundariesDescription from conmech.mesh.boundaries_factory import BoundariesFactory -from tests.test_conmech.regression.std_boundary import extract_boundary_paths_from_elements +from tests.test_conmech.regression.std_boundary import ( + extract_boundary_paths_from_elements, +) def test_identify_surfaces(): @@ -215,7 +217,9 @@ def test_condition_boundaries(_test_name_, params): elements, boundaries_data, ) = BoundariesFactory.identify_boundaries_and_reorder_nodes( - unordered_nodes, unordered_elements, boundaries_description=boundaries_description + unordered_nodes, + unordered_elements, + boundaries_description=boundaries_description, ) # Assert diff --git a/tests/test_conmech/unit_test/test_setting_mesh.py b/tests/test_conmech/unit_test/test_setting_mesh.py index b4564557..b7c85b59 100644 --- a/tests/test_conmech/unit_test/test_setting_mesh.py +++ b/tests/test_conmech/unit_test/test_setting_mesh.py @@ -5,7 +5,10 @@ from conmech.mesh import mesh from conmech.mesh.boundaries_description import BoundariesDescription from conmech.mesh.mesh import Mesh -from conmech.properties.mesh_description import RectangleMeshDescription, CubeMeshDescription +from conmech.properties.mesh_description import ( + RectangleMeshDescription, + CubeMeshDescription, +) from conmech.simulations.problem_solver import Body