Refactor time-integrator storage resolution

fedoo/constitutivelaw/beam.py

@@ -61,6 +61,40 @@ def __init__(self, material, A, Jx, Iyy, Izz, k=0, name=""):
 
         ConstitutiveLaw.__init__(self, name)  # heritage
 
+    @property
+    def linear_density(self):
+        """Mass per unit beam length."""
+        if not hasattr(self, "_linear_density"):
+            self._linear_density = self.compute_linear_density()
+        return self._linear_density
+
+    @property
+    def rotary_density(self):
+        """Mass moments of inertia per unit beam length."""
+        if not hasattr(self, "_rotary_density"):
+            self._rotary_density = self.compute_rotary_density()
+        return self._rotary_density
+
+    def compute_linear_density(self):
+        density = self._material_density()
+        return density * self.A
+
+    def compute_rotary_density(self):
+        density = self._material_density()
+        return [density * self.Jx, density * self.Iyy, density * self.Izz]
+
+    def _material_density(self):
+        density = getattr(self.material, "density", None)
+        if density is None:
+            material_name = getattr(self.material, "name", type(self.material).__name__)
+            raise ValueError(
+                f"Beam properties {self.name!r} need density on material "
+                f"{material_name!r} for dynamic analysis. Set it with "
+                "material.set_density(rho), or attach storage explicitly "
+                "with weakform.set_inertia(...)."
+            )
+        return density
+
     def get_beam_rigidity(self):
         E = self.material.E
         G = self.material.G

fedoo/constitutivelaw/shell.py

@@ -22,6 +22,43 @@ def get_shell_stiffness_matrix(self):
             '"get_shell_stiffness_matrix" not implemented, contact developer.'
         )
 
+    @property
+    def area_density(self):
+        """Mass per unit midsurface area."""
+        if not hasattr(self, "_area_density"):
+            self._area_density = self.compute_area_density()
+        return self._area_density
+
+    @property
+    def rotary_density(self):
+        """Mass moment of inertia per unit midsurface area."""
+        if not hasattr(self, "_rotary_density"):
+            self._rotary_density = self.compute_rotary_density()
+        return self._rotary_density
+
+    def compute_area_density(self):
+        raise NotImplementedError(
+            f"{type(self).__name__} does not define shell area density."
+        )
+
+    def compute_rotary_density(self):
+        raise NotImplementedError(
+            f"{type(self).__name__} does not define shell rotary density."
+        )
+
+    @staticmethod
+    def _material_density(material, owner_name):
+        density = getattr(material, "density", None)
+        if density is None:
+            material_name = getattr(material, "name", type(material).__name__)
+            raise ValueError(
+                f"Shell law {owner_name!r} needs density on material "
+                f"{material_name!r} for dynamic analysis. Set it with "
+                "material.set_density(rho), or attach storage explicitly "
+                "with weakform.set_inertia(...)."
+            )
+        return density
+
     def update(self, assembly, pb):
         # disp = pb.get_disp()
         # rot = pb.get_rot()
@@ -100,6 +137,14 @@ def __init__(self, material, thickness, k=1, name=""):
 
         self.material = material
 
+    def compute_area_density(self):
+        density = self._material_density(self.material, self.name)
+        return density * self.thickness
+
+    def compute_rotary_density(self):
+        density = self._material_density(self.material, self.name)
+        return density * self.thickness**3 / 12.0
+
     def get_shell_stiffness_matrix(self):
         Hplane = self.material.get_elastic_matrix(
             "2Dstress"
@@ -221,6 +266,20 @@ def __init__(self, listMat, list_thickness, k=1, name=""):
 
         ShellBase.__init__(self, thickness, k, name)  # heritage
 
+    def compute_area_density(self):
+        return sum(
+            self._material_density(material, self.name) * thickness
+            for material, thickness in zip(self.__listMat, self.list_thickness)
+        )
+
+    def compute_rotary_density(self):
+        return sum(
+            self._material_density(material, self.name)
+            * (self.__layer[i + 1] ** 3 - self.__layer[i] ** 3)
+            / 3.0
+            for i, material in enumerate(self.__listMat)
+        )
+
     def get_shell_stiffness_matrix(self):
         H = np.zeros((8, 8), dtype="object")
         for i in range(len(self.list_thickness)):

fedoo/core/weakform.py

@@ -141,6 +141,15 @@ def set_storage(self, storage, evolution=None):
             self.time_evolution = normalize_time_evolution(evolution)
         return self
 
+    def get_storage(self):
+        """Return the optional storage weakform for a time integrator.
+
+        Subclasses may derive storage from their constitutive data (for
+        example inertia from material density). Return ``None`` when the
+        weakform has no transient storage contribution.
+        """
+        return self.storage
+
     def set_dissipation(self, dissipation=None, evolution=None, **kargs):
         """Attach a dissipative contribution.
 
@@ -174,6 +183,10 @@ def set_dissipation(self, dissipation=None, evolution=None, **kargs):
             self.time_evolution = normalize_time_evolution(evolution)
         return self
 
+    def get_dissipation(self):
+        """Return the optional dissipative provider for a time integrator."""
+        return self.dissipation
+
     def set_inertia(self, density_or_storage):
         """Mechanical alias for ``set_storage``.
 

fedoo/time/backward_euler.py

@@ -14,14 +14,14 @@ class BackwardEuler(TimeIntegratorBase):
 
     def _integrate_leaf(self, weakform):
         """Return a backward-Euler weakform sum when storage is declared."""
-        if getattr(weakform, "dissipation", None) is not None:
+        if weakform.get_dissipation() is not None:
             raise NotImplementedError(
                 "BackwardEuler does not integrate dissipative terms yet. "
                 "Remove set_dissipation() from this first-order weakform, or "
                 "fold the dissipation into the storage term."
             )
 
-        storage = getattr(weakform, "storage", None)
+        storage = weakform.get_storage()
         if storage is None:
             # A first-order weakform without a storage term has no transient
             # contribution: keep it as a plain static term.

fedoo/time/generalized_alpha.py

@@ -76,7 +76,7 @@ def _integrate_leaf(self, weakform):
                 "weakform.set_inertia(density_or_weakform)."
             )
 
-        dissipation = getattr(weakform, "dissipation", None)
+        dissipation = self._resolve_dissipation(weakform)
         if dissipation is not None and not isinstance(
             dissipation, (RayleighDamping, WeakFormBase)
         ):
@@ -89,19 +89,15 @@ def _integrate_leaf(self, weakform):
         return self._wrap_static_weakform(weakform, storage, dissipation)
 
     def _resolve_storage(self, weakform):
-        storage = getattr(weakform, "storage", None)
-        if storage is not None:
-            if isinstance(storage, WeakFormBase):
-                return storage
-            return Inertia(storage, space=weakform.space)
-
-        # Fall back to the material density, read at compile time so that a
-        # set_density() call made after the weakform was built is still honored.
-        constitutivelaw = getattr(weakform, "constitutivelaw", None)
-        density = getattr(constitutivelaw, "density", None)
-        if density is not None:
-            return Inertia(density, space=weakform.space)
-        return None
+        storage = weakform.get_storage()
+        if isinstance(storage, WeakFormBase):
+            return storage
+        return Inertia(storage, space=weakform.space)
+
+    def _resolve_dissipation(self, weakform):
+        if getattr(weakform, "dissipation", None) is not None:
+            return weakform.dissipation
+        return weakform.get_dissipation()
 
     def _wrap_static_weakform(self, weakform, storage, dissipation):
         # Decorate a *copy* of the user's weakform with the generalized-alpha
@@ -283,6 +279,15 @@ def get_weak_equation(self, assembly, pb):
 class GeneralizedAlphaWeakFormSum(WeakFormSum):
     """WeakFormSum with Rayleigh damping accessors for second-order terms."""
 
+    def __getattr__(self, name):
+        if name.startswith("__"):
+            raise AttributeError(name)
+        try:
+            stiffness = self.list_weakform[0]
+        except (AttributeError, IndexError):
+            raise AttributeError(name) from None
+        return getattr(stiffness, name)
+
     def _rayleigh_terms(self):
         """Return the [stiffness, inertia] terms if they carry damping_coef.
 

fedoo/weakform/beam.py

@@ -75,26 +75,6 @@ def __init__(
                 material, A, Jx, Iyy, Izz, k, name + "_properties"
             )
         self.time_evolution = SECOND_ORDER
-        density = getattr(self.properties.material, "density", None)
-        if density is not None:
-            translational_inertia = Inertia(
-                density * self.properties.A, space=self.space
-            )
-            if self.space.ndim == 3:
-                rotary_inertia = RotaryInertia(
-                    [
-                        density * self.properties.Jx,
-                        density * self.properties.Iyy,
-                        density * self.properties.Izz,
-                    ],
-                    space=self.space,
-                )
-            else:
-                rotary_inertia = RotaryInertia(
-                    density * self.properties.Izz,
-                    space=self.space,
-                )
-            self.storage = WeakFormSum([translational_inertia, rotary_inertia])
         self.assembly_options.set("elm_type", "beam", elm_type="lin2")
 
         self.nlgeom = nlgeom  # geometric non linearities -> False, True, 'UL' or 'TL' (True or 'UL': updated lagrangian - 'TL': total lagrangian)
@@ -106,6 +86,18 @@ def __init__(
               initial mesh with initial displacement effet)
         """
 
+    def get_storage(self):
+        if self.storage is not None:
+            return self.storage
+        translational_inertia = Inertia(
+            self.properties.linear_density, space=self.space
+        )
+        rotary_density = self.properties.rotary_density
+        if self.space.ndim == 2:
+            rotary_density = rotary_density[2]
+        rotary_inertia = RotaryInertia(rotary_density, space=self.space)
+        return WeakFormSum([translational_inertia, rotary_inertia])
+
     def initialize(self, assembly, pb):
         assembly.sv["BeamStrain"] = 0
         assembly.sv["BeamStress"] = 0

fedoo/weakform/plate.py

@@ -1,5 +1,7 @@
-from fedoo.core.weakform import WeakFormBase
+from fedoo.core.weakform import WeakFormBase, WeakFormSum
 from fedoo.core.base import ConstitutiveLaw
+from fedoo.core.time_evolution import SECOND_ORDER
+from fedoo.weakform.inertia import Inertia, RotaryInertia
 from scipy.spatial.transform import Rotation
 import numpy as np
 
@@ -89,6 +91,16 @@ def __init__(
         self.true_drilling_rotation = true_drilling_rotation
         self.drill_stiffness_coefficient = drill_stiffness_coefficient
         self._store_local_pos = False
+        if type(self).__name__ in ("PlateEquilibriumFI", "PlateEquilibrium"):
+            self.time_evolution = SECOND_ORDER
+
+    def get_storage(self):
+        if self.storage is not None:
+            return self.storage
+        shell = self.constitutivelaw
+        translational_inertia = Inertia(shell.area_density, space=self.space)
+        rotary_inertia = RotaryInertia(shell.rotary_density, space=self.space)
+        return WeakFormSum([translational_inertia, rotary_inertia])
 
     def initialize(self, assembly, pb):
         assembly.sv["ShellStrain"] = 0

fedoo/weakform/stress_equilibrium.py

@@ -3,6 +3,7 @@
 from fedoo.core.weakform import WeakFormBase
 from fedoo.core.base import ConstitutiveLaw
 from fedoo.core.time_evolution import SECOND_ORDER
+from fedoo.weakform.inertia import Inertia
 from fedoo.util.voigt_tensors import StressTensorList, StrainTensorList
 import numpy as np
 import simcoon as sim
@@ -89,6 +90,22 @@ def __init__(self, constitutivelaw, name="", nlgeom=None, space=None):
         # (if TangentMatrix is symmetric)
         # -> need to be checked for general case
 
+    def get_storage(self):
+        if self.storage is not None:
+            return self.storage
+        density = getattr(self.constitutivelaw, "density", None)
+        if density is None:
+            material_name = getattr(
+                self.constitutivelaw, "name", type(self.constitutivelaw).__name__
+            )
+            raise ValueError(
+                "StressEquilibrium requires a material density for dynamic "
+                f"analysis, but material {material_name!r} has no density. "
+                "Set it with material.set_density(rho), or attach inertia "
+                "explicitly with weakform.set_inertia(density_or_weakform)."
+            )
+        return Inertia(density, space=self.space)
+
     def get_weak_equation(self, assembly, pb):
         """Get the weak equation related to the current problem state."""
         if assembly._nlgeom == "TL":  # add initial displacement effect
@@ -666,52 +683,3 @@ def _comp_gn_strain(wf, assembly, pb):
     )
     assembly.sv["DR"] = DR
     assembly.sv["DStrain"] = StrainTensorList(DStrain)
-
-
-# def _comp_linear_strain_pgd(wf, assembly, pb):
-#     # may be compatible with other methods like PGD
-#     # but not compatible with simcoon
-#     assert not (
-#         wf.nlgeom
-#     ), "the current strain measure isn't adapted for finite strain"
-#     grad_values = assembly.sv["DispGradient"]
-
-#     strain = [grad_values[i][i] for i in range(3)]
-#     strain += [
-#         grad_values[0][1] + grad_values[1][0],
-#         grad_values[0][2] + grad_values[2][0],
-#         grad_values[1][2] + grad_values[2][1],
-#     ]
-#     assembly.sv["Strain"] = StrainTensorList(strain)
-
-
-# def _comp_gl_strain(wf, assembly, pb):
-#     # not compatible with simcoon
-#     if not (wf.nlgeom):
-#         return _comp_linear_strain_pgd(wf, assembly, pb)
-#     else:
-#         grad_values = assembly.sv["DispGradient"]
-#         # GL strain tensor
-#         # possibility to be improve from simcoon functions
-#         # to get the logarithmic strain tensor...
-#         strain = [
-#             grad_values[i][i]
-#             + 0.5 * sum([grad_values[k][i] ** 2 for k in range(3)])
-#             for i in range(3)
-#         ]
-#         strain += [
-#             grad_values[0][1]
-#             + grad_values[1][0]
-#             + sum([grad_values[k][0] * grad_values[k][1] for k in range(3)])
-#         ]
-#         strain += [
-#             grad_values[0][2]
-#             + grad_values[2][0]
-#             + sum([grad_values[k][0] * grad_values[k][2] for k in range(3)])
-#         ]
-#         strain += [
-#             grad_values[1][2]
-#             + grad_values[2][1]
-#             + sum([grad_values[k][1] * grad_values[k][2] for k in range(3)])
-#         ]
-#         return StrainTensorList(strain)