From 3806bc457fbf2da65dc4ad6c5b7f3178a679ec13 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Wed, 27 Nov 2024 17:32:57 +0100 Subject: [PATCH 01/10] Add `FreeVibration(items, boundaries)` with methods to `evaluate()` and `extract()` the n-th eigenvector to a field --- src/felupe/__init__.py | 2 + src/felupe/mechanics/__init__.py | 2 + src/felupe/mechanics/_free_vibration.py | 118 ++++++++++++++++++ src/felupe/mechanics/_solidbody.py | 13 +- .../mechanics/_solidbody_incompressible.py | 12 +- 5 files changed, 143 insertions(+), 4 deletions(-) create mode 100644 src/felupe/mechanics/_free_vibration.py diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 91427396..dd1c2fcd 100644 --- a/src/felupe/__init__.py +++ b/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", diff --git a/src/felupe/mechanics/__init__.py b/src/felupe/mechanics/__init__.py index 9b14d26b..69266ada 100644 --- a/src/felupe/mechanics/__init__.py +++ b/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", diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py new file mode 100644 index 00000000..04dc13d5 --- /dev/null +++ b/src/felupe/mechanics/_free_vibration.py @@ -0,0 +1,118 @@ +# -*- 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 . +""" + +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. + + 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). + + Examples + -------- + .. pyvista-plot:: + + >>> import felupe as fem + >>> import numpy as np + >>> from scipy.sparse.linalg import eigsh + >>> + >>> 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)) + >>> dof0, dof1 = fem.dof.partition(field, boundaries) + >>> + >>> solid = fem.SolidBody(umat=fem.LinearElastic(E=2.5, nu=0.25), field=field) + """ + + 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, k=6, solver=eigsh, tol=0, **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(**kwargs) + 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] + + self.eigenvalues, self.eigenvectors = solver(A=K, k=k, M=M, sigma=0, tol=tol) + + return self + + def extract(self, n=0, x0=None, inplace=True): + if x0 is not None: + x = x0 + else: + x = self.items[0].field + + if not inplace: + x = x.copy() + + frequency = np.sqrt(self.eigenvalues[n]) / (2 * np.pi) + x[0].values.ravel()[self.dof1] = self.eigenvectors[:, n] + + return x, frequency diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index e7b9cd7b..66abd50c 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -132,6 +132,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 ----- @@ -209,9 +211,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) @@ -364,9 +367,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, @@ -375,4 +383,5 @@ def _mass(self, density=1.0): grad_v=[False], grad_u=[False], ) + return form.assemble() diff --git a/src/felupe/mechanics/_solidbody_incompressible.py b/src/felupe/mechanics/_solidbody_incompressible.py index dfa5264b..e52aa69c 100644 --- a/src/felupe/mechanics/_solidbody_incompressible.py +++ b/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() From 94d0b8a050619de3482f7704a272376f72c0bb39 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 16:17:58 +0100 Subject: [PATCH 02/10] Update _free_vibration.py --- src/felupe/mechanics/_free_vibration.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py index 04dc13d5..42e101d0 100644 --- a/src/felupe/mechanics/_free_vibration.py +++ b/src/felupe/mechanics/_free_vibration.py @@ -61,7 +61,7 @@ def __init__(self, items, boundaries=None): self.eigenvalues = None self.eigenvectors = None - def evaluate(self, x0=None, k=6, solver=eigsh, tol=0, **kwargs): + def evaluate(self, x0=None, solver=eigsh, parallel=False, **kwargs): if x0 is not None: x = x0 else: @@ -77,7 +77,7 @@ def evaluate(self, x0=None, k=6, solver=eigsh, tol=0, **kwargs): # assemble matrices for item in self.items: - K = item.assemble.matrix(**kwargs) + K = item.assemble.matrix(parallel=parallel) M = item.assemble.mass() if item.assemble.multiplier is not None: @@ -99,20 +99,26 @@ def evaluate(self, x0=None, k=6, solver=eigsh, tol=0, **kwargs): K = stiffness[self.dof1][:, self.dof1] M = mass[self.dof1][:, self.dof1] - self.eigenvalues, self.eigenvectors = solver(A=K, k=k, M=M, sigma=0, tol=tol) + 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: - x = x0 + field = x0 else: - x = self.items[0].field + field = self.items[0].field if not inplace: - x = x.copy() + field = field.copy() + + values = np.zeros(sum(field.fieldsizes)) + values[self.dof1] = self.eigenvectors[:, n] frequency = np.sqrt(self.eigenvalues[n]) / (2 * np.pi) - x[0].values.ravel()[self.dof1] = self.eigenvectors[:, n] - return x, frequency + [f.fill(0) for f in field] + field += values + + return field, frequency From 0939857bb8e02592dfd0f28db5d446545f457711 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 17:13:04 +0100 Subject: [PATCH 03/10] Create test_free_vibration.py --- tests/test_free_vibration.py | 40 ++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 tests/test_free_vibration.py diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py new file mode 100644 index 00000000..8c02a18f --- /dev/null +++ b/tests/test_free_vibration.py @@ -0,0 +1,40 @@ +import numpy as np +import felupe as fem + + +def test_free_vibration(): + 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) + + assert np.isclose(new_field[0].values.max(), 358.08157) + + +if __name__ == "__main__": + test_free_vibration() From 4548a8f465075135e10f5f4fb6563b670f9c2297 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 17:13:27 +0100 Subject: [PATCH 04/10] add test --- src/felupe/mechanics/_free_vibration.py | 1 + tests/test_free_vibration.py | 1 + 2 files changed, 2 insertions(+) diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py index 42e101d0..dc28ff6f 100644 --- a/src/felupe/mechanics/_free_vibration.py +++ b/src/felupe/mechanics/_free_vibration.py @@ -19,6 +19,7 @@ import numpy as np from scipy.sparse import csr_matrix from scipy.sparse.linalg import eigsh + from ..dof import partition diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index 8c02a18f..c4aec5d1 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -1,4 +1,5 @@ import numpy as np + import felupe as fem From bae3834d2c6c8eaee97b3a54d9eed09fe0d24e19 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 19:20:57 +0100 Subject: [PATCH 05/10] Update test_mechanics.py --- tests/test_mechanics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 4039afc5..60baa44c 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -235,7 +235,7 @@ 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() From 5692cacf60927d4cd6fbba35c7e9eb5a53739960 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 19:22:01 +0100 Subject: [PATCH 06/10] Update test_mechanics.py --- tests/test_mechanics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index 60baa44c..f27ed6a0 100644 --- a/tests/test_mechanics.py +++ b/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() From d5b28a202a0b7f1e72ab6fef57eac227a92193c4 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 20:54:19 +0100 Subject: [PATCH 07/10] add test, format black --- src/felupe/mechanics/_solidbody.py | 8 ++++---- tests/test_free_vibration.py | 17 +++++++++++++++-- tests/test_mechanics.py | 6 +++++- 3 files changed, 24 insertions(+), 7 deletions(-) diff --git a/src/felupe/mechanics/_solidbody.py b/src/felupe/mechanics/_solidbody.py index 7bdaa4af..fa806639 100644 --- a/src/felupe/mechanics/_solidbody.py +++ b/src/felupe/mechanics/_solidbody.py @@ -398,13 +398,13 @@ def _cauchy_stress(self, field=None): return dot(P, transpose(F)) / J 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, @@ -413,5 +413,5 @@ def _mass(self, density=None): grad_v=[False], grad_u=[False], ) - + return form.assemble() diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index c4aec5d1..e74129a3 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -4,6 +4,20 @@ 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(x0=field) + new_field, frequency = job.extract(x0=field, 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)), @@ -34,8 +48,7 @@ def test_free_vibration(): job = fem.FreeVibration(solids, boundaries).evaluate(x0=field) new_field, frequency = job.extract(x0=field, n=-1, inplace=False) - assert np.isclose(new_field[0].values.max(), 358.08157) - if __name__ == "__main__": test_free_vibration() + test_free_vibration_mixed() diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py index f27ed6a0..fdf3daa8 100644 --- a/tests/test_mechanics.py +++ b/tests/test_mechanics.py @@ -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), density=1.0 + umat=umat, + field=u, + bulk=5000, + state=fem.StateNearlyIncompressible(u), + density=1.0, ) M = b.assemble.mass() From 3ce9893ce3f5ca5f6e18988900714b12bf51fbe6 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 20:54:40 +0100 Subject: [PATCH 08/10] Update test_free_vibration.py --- tests/test_free_vibration.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index e74129a3..43c9f271 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -1,3 +1,30 @@ +# -*- 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 . + +""" + import numpy as np import felupe as fem From ad36bc111b455511e1c6606a46fd3063c329412c Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 21:08:22 +0100 Subject: [PATCH 09/10] update the free-vibration test --- src/felupe/mechanics/_free_vibration.py | 17 ++++++++++++----- tests/test_free_vibration.py | 4 ++-- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/felupe/mechanics/_free_vibration.py b/src/felupe/mechanics/_free_vibration.py index dc28ff6f..fa149b88 100644 --- a/src/felupe/mechanics/_free_vibration.py +++ b/src/felupe/mechanics/_free_vibration.py @@ -24,7 +24,7 @@ class FreeVibration: - """A Free-Vibration Step. + """A Free-Vibration Step/Job. Parameters ---------- @@ -33,6 +33,12 @@ class FreeVibration: 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 -------- @@ -40,16 +46,17 @@ class FreeVibration: >>> import felupe as fem >>> import numpy as np - >>> from scipy.sparse.linalg import eigsh >>> >>> 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)) - >>> dof0, dof1 = fem.dof.partition(field, boundaries) - >>> - >>> solid = fem.SolidBody(umat=fem.LinearElastic(E=2.5, nu=0.25), field=field) + >>> 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): diff --git a/tests/test_free_vibration.py b/tests/test_free_vibration.py index 43c9f271..83708194 100644 --- a/tests/test_free_vibration.py +++ b/tests/test_free_vibration.py @@ -39,8 +39,8 @@ def test_free_vibration(): umat = fem.NeoHooke(mu=1, bulk=2) solid = fem.SolidBody(umat=umat, field=field, density=1.5e-9) - job = fem.FreeVibration([solid]).evaluate(x0=field) - new_field, frequency = job.extract(x0=field, n=-1, inplace=False) + job = fem.FreeVibration([solid]).evaluate() + new_field, frequency = job.extract(n=-1, inplace=False) def test_free_vibration_mixed(): From 95c965d29da5da9424fe0a10fa2f4e55e46ee000 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 28 Nov 2024 21:09:36 +0100 Subject: [PATCH 10/10] Update mechanics.rst --- docs/felupe/mechanics.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/felupe/mechanics.rst b/docs/felupe/mechanics.rst index 7afe0792..2f9d4d8e 100644 --- a/docs/felupe/mechanics.rst +++ b/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: