Density

conmech/dynamics/dynamics.py

@@ -89,7 +89,8 @@ def __init__(
         self.acceleration_operator: np.ndarray
         self.elasticity: np.ndarray
         self.viscosity: np.ndarray
-        self._w_matrix = None
+        self._w_matrix: np.ndarray = None
+        self._local_stifness_matrices: np.ndarray = None
         self.__relaxation: Optional[np.ndarray] = None
         self.__relaxation_tensor: Optional[float] = None
         self.thermal_expansion: np.ndarray
@@ -107,14 +108,27 @@ def remesh(self, boundaries_description, create_in_subprocess):
         super().mesh.remesh(boundaries_description, create_in_subprocess)
         self.reinitialize_matrices()
 
-    def reinitialize_matrices(self):
+    def reinitialize_matrices(self, elements_density: Optional[np.ndarray] = None):
         (
             self.element_initial_volume,
             self.volume_at_nodes,
             U,
             V,
             self._w_matrix,
+            self._local_stifness_matrices,
         ) = get_basic_matrices(elements=self.mesh.elements, nodes=self.moved_nodes)
+
+        if elements_density is not None:
+            self._w_matrix.fill(0)
+            for element_index, element in enumerate(self.mesh.elements):
+                for i, global_i in enumerate(element):
+                    for j, global_j in enumerate(element):
+                        for w_idx in range(4):
+                            self._w_matrix[w_idx // 2][w_idx % 2][global_i, global_j] += (
+                                elements_density[element_index]
+                                * self._local_stifness_matrices[w_idx, element_index, i, j]
+                            )
+
         (
             self.acceleration_operator,
             self.elasticity,

conmech/dynamics/factory/_dynamics_factory_2d.py

@@ -23,6 +23,7 @@ def get_edges_features_matrix_numba(elements, nodes):
         (FEATURE_MATRIX_COUNT, nodes_count, nodes_count), dtype=np.double
     )
     element_initial_volume = np.zeros(elements_count)
+    local_stifness_matrices = np.empty((4, elements_count, element_size, element_size))
 
     for element_index in range(elements_count):  # TODO: #65 prange?
         element = elements[element_index]
@@ -85,6 +86,10 @@ def get_edges_features_matrix_numba(elements, nodes):
 
                 w = [[i_d_phi * j_d_phi for j_d_phi in j_d_phi_vec] for i_d_phi in i_d_phi_vec]
 
+                local_stifness_matrices[:, element_index, i, j] = (
+                    element_volume * np.asarray(w).flatten()
+                )
+
                 edges_features_matrix[:, element[i], element[j]] += element_volume * np.array(
                     [
                         volume_at_nodes,
@@ -100,7 +105,7 @@ def get_edges_features_matrix_numba(elements, nodes):
                 )
 
     # Performance TIP: we need only sparse, triangular matrix (?)
-    return edges_features_matrix, element_initial_volume
+    return edges_features_matrix, element_initial_volume, local_stifness_matrices
 
 
 @numba.njit

conmech/dynamics/factory/_dynamics_factory_3d.py

@@ -68,7 +68,7 @@ def get_edges_features_matrix_numba(elements, nodes):
                     ]
                 )
 
-    return edges_features_matrix, element_initial_volume
+    return edges_features_matrix, element_initial_volume, None
 
 
 @numba.njit

conmech/dynamics/factory/dynamics_factory_method.py

@@ -25,9 +25,11 @@ def get_basic_matrices(elements: np.ndarray, nodes: np.ndarray):
     dimension = len(elements[0]) - 1
     factory = get_factory(dimension)
 
-    edges_features_matrix, element_initial_volume = factory.get_edges_features_matrix(
-        elements, nodes
-    )
+    (
+        edges_features_matrix,
+        element_initial_volume,
+        local_stifness_matrices,
+    ) = factory.get_edges_features_matrix(elements, nodes)
 
     volume_at_nodes = edges_features_matrix[0]
     U = edges_features_matrix[1]
@@ -42,7 +44,7 @@ def get_basic_matrices(elements: np.ndarray, nodes: np.ndarray):
             for k in range(factory.dimension)
         ]
     )
-    return element_initial_volume, volume_at_nodes, U, V, W
+    return element_initial_volume, volume_at_nodes, U, V, W, local_stifness_matrices
 
 
 def get_dynamics(elements: np.ndarray, body_prop: BodyProperties, U, V, W):

conmech/simulations/problem_solver.py

@@ -109,6 +109,10 @@ def solving_method(self, value: str) -> None:
         self.__set_step_solver(value)
         self.__set_second_step_solver(value)
 
+    def refresh_solvers(self):
+        self.__set_step_solver(self.solving_method)
+        self.__set_second_step_solver(self.solving_method)
+
     def __set_step_solver(self, value):
         solver_class: Type[Solver] = SolversRegistry.get_by_name(
             solver_name=value, problem=self.problem
@@ -376,6 +380,12 @@ def solve(self, *, initial_displacement: Callable, verbose: bool = False, **kwar
         return state
 
 
+class NonHomogenousSolver(StaticSolver):
+    def update_density(self, density: np.ndarray):
+        self.body.reinitialize_matrices(density)
+        self.refresh_solvers()
+
+
 class QuasistaticRelaxation(ProblemSolver):
     def __init__(self, setup: RelaxationQuasistaticProblem, solving_method: str):
         """Solves general Contact Mechanics problem.

examples/example_density.py

@@ -0,0 +1,63 @@
+"""
+Created at 02.09.2023
+"""
+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 StaticDisplacementProblem
+from conmech.simulations.problem_solver import NonHomogenousSolver
+
+from examples.p_slope_contact_law import make_slope_contact_law
+
+
+@dataclass
+class StaticSetup(StaticDisplacementProblem):
+    grid_height: ... = 1.0
+    elements_number: ... = (2, 5)
+    mu_coef: ... = 4
+    la_coef: ... = 4
+    contact_law: ... = make_slope_contact_law(slope=1)
+
+    @staticmethod
+    def inner_forces(x, t=None):
+        return np.array([0, 0])
+
+    @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
+
+    boundaries: ... = BoundariesDescription(
+        contact=lambda x: x[1] == 0 and x[0] < 0.2, dirichlet=lambda x: x[0] == 0
+    )
+
+
+def main(config: Config):
+    """
+    Entrypoint to example.
+
+    To see result of simulation you need to call from python `main(Config().init())`.
+    """
+    setup = StaticSetup(mesh_type="cross")
+    runner = NonHomogenousSolver(setup, "schur")
+
+    elem_centers = np.empty(shape=(len(runner.body.mesh.elements), 2))
+    for idx, elem in enumerate(runner.body.mesh.elements):
+        verts = runner.body.mesh.initial_nodes[elem]
+        elem_centers[idx] = np.sum(verts, axis=0) / len(elem)
+    density = np.asarray([1 if x[0] < 1 else 0.2 for x in elem_centers])
+
+    runner.update_density(density)
+
+    state = runner.solve(verbose=True, initial_displacement=setup.initial_displacement)
+    Drawer(state=state, config=config).draw(show=config.show, save=config.save)
+
+
+if __name__ == "__main__":
+    main(Config().init())

tests/test_conmech/unit_test/test_matrices.py

@@ -22,7 +22,9 @@ def test_matrices_2d_integrals():
     )
 
     # Act
-    edges_features_matrix, element_initial_volume = sut_2d(elements=elements, nodes=initial_nodes)
+    edges_features_matrix, element_initial_volume, _ = sut_2d(
+        elements=elements, nodes=initial_nodes
+    )
 
     # Assert
     np.testing.assert_allclose(element_initial_volume.sum(), area)
@@ -50,7 +52,9 @@ def test_matrices_3d_integrals():
     )
 
     # Act
-    edges_features_matrix, element_initial_volume = sut_3d(elements=elements, nodes=initial_nodes)
+    edges_features_matrix, element_initial_volume, _ = sut_3d(
+        elements=elements, nodes=initial_nodes
+    )
 
     # Assert
     np.testing.assert_allclose(element_initial_volume.sum(), 1)