Add FreeVibration(items, boundaries) (#909)

* Add `FreeVibration(items, boundaries)`

with methods to `evaluate()` and `extract()` the n-th eigenvector to a field

* Update _free_vibration.py

* Create test_free_vibration.py

* add test

* Update test_mechanics.py

* Update test_mechanics.py

* add test, format black

* Update test_free_vibration.py

* update the free-vibration test

* Update mechanics.rst

docs/felupe/mechanics.rst

@@ -22,6 +22,7 @@ Mechanics
    Step
    Job
    CharacteristicCurve
+   FreeVibration
 
 **Point Load and Multi-Point Constraints**
 
@@ -87,6 +88,11 @@ Mechanics
    :undoc-members:
    :show-inheritance:
 
+.. autoclass:: felupe.FreeVibration
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 .. autoclass:: felupe.MultiPointConstraint
    :members:
    :undoc-members:

src/felupe/__init__.py

@@ -70,6 +70,7 @@
 from .mechanics import (
     CharacteristicCurve,
     FormItem,
+    FreeVibration,
     Job,
     MultiPointConstraint,
     MultiPointContact,
@@ -133,6 +134,7 @@
     "tools",
     "Form",
     "FormItem",
+    "FreeVibration",
     "IntegralForm",
     "Basis",
     "Field",

src/felupe/mechanics/__init__.py

@@ -1,4 +1,5 @@
 from ._curve import CharacteristicCurve
+from ._free_vibration import FreeVibration
 from ._helpers import Assemble, Evaluate, Results, StateNearlyIncompressible
 from ._item import FormItem
 from ._job import Job
@@ -16,6 +17,7 @@
     "Assemble",
     "CharacteristicCurve",
     "Evaluate",
+    "FreeVibration",
     "FormItem",
     "StateNearlyIncompressible",
     "Job",

src/felupe/mechanics/_free_vibration.py

@@ -0,0 +1,132 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe 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.
+
+FElupe 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 FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import numpy as np
+from scipy.sparse import csr_matrix
+from scipy.sparse.linalg import eigsh
+
+from ..dof import partition
+
+
+class FreeVibration:
+    """A Free-Vibration Step/Job.
+
+    Parameters
+    ----------
+    items : list of SolidBody or SolidBodyNearlyIncompressible
+        A list of items with methods for the assembly of sparse stiffness and mass
+        matrices.
+    boundaries : dict of Boundary, optional
+        A dict with :class:`~felupe.Boundary` conditions (default is None).
+    
+    Notes
+    -----
+    ..  note::
+        
+        Boundary conditions with non-zero values are not supported.
+
+    Examples
+    --------
+    ..  pyvista-plot::
+
+        >>> import felupe as fem
+        >>> import numpy as np
+        >>>
+        >>> mesh = fem.Rectangle(b=(5, 1), n=(50, 10))
+        >>> region = fem.RegionQuad(mesh)
+        >>> field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])
+        >>>
+        >>> boundaries = dict(left=fem.Boundary(field[0], fx=0))
+        >>> solid = fem.SolidBody(fem.LinearElastic(E=2.5, nu=0.25), field, density=1.0)
+        >>> modal = fem.FreeVibration(items=[solid], boundaries=boundaries).evaluate()
+        >>> 
+        >>> eigenvector, frequency = modal.extract(n=4, inplace=True)
+        >>> solid.plot("Stress", component=0).show()
+    """
+
+    def __init__(self, items, boundaries=None):
+        self.items = items
+
+        if boundaries is None:
+            boundaries = {}
+
+        self.boundaries = boundaries
+        self.eigenvalues = None
+        self.eigenvectors = None
+
+    def evaluate(self, x0=None, solver=eigsh, parallel=False, **kwargs):
+        if x0 is not None:
+            x = x0
+        else:
+            x = self.items[0].field
+
+        # link field of items with global field
+        [item.field.link(x) for item in self.items]
+
+        # init matrix with shape from global field
+        shape = (np.sum(x.fieldsizes), np.sum(x.fieldsizes))
+        stiffness = csr_matrix(shape)
+        mass = csr_matrix(shape)
+
+        # assemble matrices
+        for item in self.items:
+            K = item.assemble.matrix(parallel=parallel)
+            M = item.assemble.mass()
+
+            if item.assemble.multiplier is not None:
+                K *= item.assemble.multiplier
+
+            # check and reshape matrices
+            if K.shape != stiffness.shape:
+                K.resize(*shape)
+
+            if M.shape != mass.shape:
+                M.resize(*shape)
+
+            # add matrices
+            stiffness += K
+            mass += M
+
+        dof0, self.dof1 = partition(x, self.boundaries)
+
+        K = stiffness[self.dof1][:, self.dof1]
+        M = mass[self.dof1][:, self.dof1]
+
+        sigma = kwargs.get("sigma", 0)
+        self.eigenvalues, self.eigenvectors = solver(A=K, M=M, sigma=sigma, **kwargs)
+
+        return self
+
+    def extract(self, n=0, x0=None, inplace=True):
+        if x0 is not None:
+            field = x0
+        else:
+            field = self.items[0].field
+
+        if not inplace:
+            field = field.copy()
+
+        values = np.zeros(sum(field.fieldsizes))
+        values[self.dof1] = self.eigenvectors[:, n]
+
+        frequency = np.sqrt(self.eigenvalues[n]) / (2 * np.pi)
+
+        [f.fill(0) for f in field]
+        field += values
+
+        return field, frequency

src/felupe/mechanics/_solidbody.py

@@ -162,6 +162,8 @@ class SolidBody(Solid):
         A field container with one or more fields.
     statevars : ndarray or None, optional
         Array of initial internal state variables (default is None).
+    density : float or None, optional
+        The density of the solid body.
 
     Notes
     -----
@@ -239,9 +241,10 @@ class SolidBody(Solid):
         methods for the assembly of sparse vectors/matrices.
     """
 
-    def __init__(self, umat, field, statevars=None):
+    def __init__(self, umat, field, statevars=None, density=None):
         self.umat = umat
         self.field = field
+        self.density = density
 
         self.results = Results(stress=True, elasticity=True)
         self.results.kinematics = self._extract(self.field)
@@ -394,9 +397,14 @@ def _cauchy_stress(self, field=None):
 
         return dot(P, transpose(F)) / J
 
-    def _mass(self, density=1.0):
+    def _mass(self, density=None):
+
+        if density is None:
+            density = self.density
+
         field = self.field[0].as_container()
         dim = field[0].dim
+
         form = self._form(
             fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)],
             v=field,
@@ -405,4 +413,5 @@ def _mass(self, density=1.0):
             grad_v=[False],
             grad_u=[False],
         )
+
         return form.assemble()

src/felupe/mechanics/_solidbody_incompressible.py

@@ -56,6 +56,8 @@ class SolidBodyNearlyIncompressible(Solid):
         A valid initial state for a (nearly) incompressible solid (default is None).
     statevars : ndarray or None, optional
         Array of initial internal state variables (default is None).
+    density : float or None, optional
+        The density of the solid body.
 
     Notes
     -----
@@ -315,10 +317,11 @@ class SolidBodyNearlyIncompressible(Solid):
         for (nearly) incompressible solid bodies.
     """
 
-    def __init__(self, umat, field, bulk, state=None, statevars=None):
+    def __init__(self, umat, field, bulk, state=None, statevars=None, density=None):
         self.umat = umat
         self.field = field
         self.bulk = bulk
+        self.density = density
 
         self._area_change = AreaChange()
         self._form = IntegralForm
@@ -540,9 +543,13 @@ def _cauchy_stress(self, field=None):
 
         return dot(P, transpose(F)) / J
 
-    def _mass(self, density=1.0):
+    def _mass(self, density=None):
+        if density is None:
+            density = self.density
+
         field = self.field[0].as_container()
         dim = field[0].dim
+
         form = self._form(
             fun=[density * np.eye(dim).reshape(dim, dim, 1, 1)],
             v=field,
@@ -551,4 +558,5 @@ def _mass(self, density=1.0):
             grad_v=[False],
             grad_u=[False],
         )
+
         return form.assemble()

tests/test_free_vibration.py

@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+"""
+ _______  _______  ___      __   __  _______  _______ 
+|       ||       ||   |    |  | |  ||       ||       |
+|    ___||    ___||   |    |  | |  ||    _  ||    ___|
+|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___ 
+|    ___||    ___||   |___ |       ||    ___||    ___|
+|   |    |   |___ |       ||       ||   |    |   |___ 
+|___|    |_______||_______||_______||___|    |_______|
+
+This file is part of felupe.
+
+Felupe 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.
+
+Felupe 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 Felupe.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+
+import numpy as np
+
+import felupe as fem
+
+
+def test_free_vibration():
+
+    mesh = fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4))
+    region = fem.RegionHexahedron(mesh)
+    field = fem.FieldContainer([fem.Field(region, dim=3)])
+
+    umat = fem.NeoHooke(mu=1, bulk=2)
+    solid = fem.SolidBody(umat=umat, field=field, density=1.5e-9)
+
+    job = fem.FreeVibration([solid]).evaluate()
+    new_field, frequency = job.extract(n=-1, inplace=False)
+
+
+def test_free_vibration_mixed():
+
+    meshes = [
+        fem.Cube(a=(0, 0, 30), b=(50, 100, 35), n=(3, 6, 2)),
+        fem.Cube(a=(0, 0, 5), b=(50, 100, 30), n=(3, 6, 4)),
+        fem.Cube(a=(0, 0, 0), b=(50, 100, 5), n=(3, 6, 2)),
+    ]
+    container = fem.MeshContainer(meshes, merge=True)
+    mesh = container.stack()
+
+    regions = [fem.RegionHexahedron(m) for m in container.meshes]
+    fields = [
+        fem.FieldsMixed(regions[0], n=1),
+        fem.FieldsMixed(regions[1], n=3),
+        fem.FieldsMixed(regions[2], n=1),
+    ]
+
+    region = fem.RegionHexahedron(mesh)
+    field = fem.FieldContainer([fem.Field(region, dim=3), *fields[1][1:]])
+
+    boundaries = dict(left=fem.Boundary(field[0], fx=0))
+    rubber = fem.ThreeFieldVariation(fem.NeoHooke(mu=1, bulk=5000))
+    steel = fem.LinearElasticLargeStrain(2.1e5, 0.3)
+    solids = [
+        fem.SolidBody(umat=steel, field=fields[0], density=7.85e-9),
+        fem.SolidBody(umat=rubber, field=fields[1], density=1.5e-9),
+        fem.SolidBody(umat=steel, field=fields[2], density=7.85e-9),
+    ]
+
+    job = fem.FreeVibration(solids, boundaries).evaluate(x0=field)
+    new_field, frequency = job.extract(x0=field, n=-1, inplace=False)
+
+
+if __name__ == "__main__":
+    test_free_vibration()
+    test_free_vibration_mixed()

tests/test_mechanics.py

@@ -39,7 +39,7 @@ def test_simple():
     r = fem.RegionHexahedron(m, uniform=True)
     u = fem.FieldContainer([fem.Field(r, dim=3)])
 
-    b = fem.SolidBody(umat, u)
+    b = fem.SolidBody(umat, u, density=1.0)
     r = b.assemble.vector()
 
     K = b.assemble.matrix()
@@ -235,7 +235,11 @@ def test_solidbody_incompressible():
 
     umat = fem.OgdenRoxburgh(fem.NeoHooke(mu=1), r=3, m=1, beta=0)
     b = fem.SolidBodyNearlyIncompressible(
-        umat=umat, field=u, bulk=5000, state=fem.StateNearlyIncompressible(u)
+        umat=umat,
+        field=u,
+        bulk=5000,
+        state=fem.StateNearlyIncompressible(u),
+        density=1.0,
     )
 
     M = b.assemble.mass()