Skip to content

Commit

Permalink
Add SolidBodyCauchyStress (#842)
Browse files Browse the repository at this point in the history
* Add `SolidBodyCauchyStress`

* Update test_mechanics.py

* Update mechanics.rst

* remove `cauchy_stress`-arg from `SolidBodyCauchyStress._gradient`

* Update test_mechanics.py

* Update test_mechanics.py
  • Loading branch information
adtzlr authored Aug 31, 2024
1 parent 5ea2566 commit 6fe8d85
Show file tree
Hide file tree
Showing 6 changed files with 214 additions and 1 deletion.
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ All notable changes to this project will be documented in this file. The format
- Add a new submodule `view` which contains the `View...` classes like `ViewSolid` or `ViewField`, previously located at `tools._plot`.
- Add the `grad`- and `hess`-arguments to the `reload()`- and `copy()`-methods of a `Region`, i.e. `Region.reload(grad=None, hess=None)`.
- Add `LinearElasticOrthotropic`.
- Add `SolidBodyCauchyStress` in addition to `SolidBodyPressure`.

### Changed
- Change the internal initialization of `field = Field(region, values=1, dtype=None)` values from `field.values = np.ones(shape) * values` to `field = np.full(shape, fill_value=values, dtype=dtype)`. This enforces `field = Field(region, values=1)` to return the gradient array with data-type `int` which was of type `float` before.
Expand Down
8 changes: 7 additions & 1 deletion docs/felupe/mechanics.rst
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,9 @@ Mechanics

SolidBody
SolidBodyNearlyIncompressible
SolidBodyPressure
SolidBodyForce
SolidBodyPressure
SolidBodyCauchyStress

**Steps and Jobs**

Expand Down Expand Up @@ -61,6 +62,11 @@ Mechanics
:undoc-members:
:show-inheritance:

.. autoclass:: felupe.SolidBodyCauchyStress
:members:
:undoc-members:
:show-inheritance:

.. autoclass:: felupe.PointLoad
:members:
:undoc-members:
Expand Down
2 changes: 2 additions & 0 deletions src/felupe/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -99,6 +99,7 @@
MultiPointContact,
PointLoad,
SolidBody,
SolidBodyCauchyStress,
SolidBodyForce,
SolidBodyGravity,
SolidBodyNearlyIncompressible,
Expand Down Expand Up @@ -243,6 +244,7 @@
"Job",
"PointLoad",
"SolidBody",
"SolidBodyCauchyStress",
"SolidBodyGravity",
"SolidBodyForce",
"SolidBodyNearlyIncompressible",
Expand Down
2 changes: 2 additions & 0 deletions src/felupe/mechanics/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
from ._multipoint import MultiPointConstraint, MultiPointContact
from ._pointload import PointLoad
from ._solidbody import SolidBody
from ._solidbody_cauchy_stress import SolidBodyCauchyStress
from ._solidbody_force import SolidBodyForce
from ._solidbody_gravity import SolidBodyGravity
from ._solidbody_incompressible import SolidBodyNearlyIncompressible
Expand All @@ -21,6 +22,7 @@
"PointLoad",
"Results",
"SolidBody",
"SolidBodyCauchyStress",
"SolidBodyGravity",
"SolidBodyForce",
"SolidBodyNearlyIncompressible",
Expand Down
162 changes: 162 additions & 0 deletions src/felupe/mechanics/_solidbody_cauchy_stress.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
# -*- 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 ..assembly import IntegralForm
from ..constitution import AreaChange
from ..math import dot
from ._helpers import Assemble, Results


class SolidBodyCauchyStress:
r"""A Cauchy stress boundary on a solid body.
Parameters
----------
field : FieldContainer
A field container with fields created on a boundary region.
stress : ndarray of shape (3, 3, ...) or None, optional
The prescribed Cauchy stress components :math:`\sigma_{ij}` (default is None).
If None, all Cauchy stress components are set to zero.
Notes
-----
.. math::
\delta W_{ext} = \int_{\partial V}
\delta \boldsymbol{u} \cdot \boldsymbol{\sigma}\ d\boldsymbol{a} =
\int_{\partial V}
\delta \boldsymbol{u} \cdot \boldsymbol{\sigma} J \boldsymbol{F}^{-T}\
d\boldsymbol{A}
Examples
--------
.. pyvista-plot::
>>> import felupe as fem
>>>
>>> mesh = fem.Rectangle(n=6)
>>> region = fem.RegionQuad(mesh)
>>> field = fem.FieldContainer([fem.FieldAxisymmetric(region, dim=2)])
>>> boundaries = fem.dof.symmetry(field[0])
>>> umat = fem.NeoHooke(mu=1, bulk=2)
>>> solid = fem.SolidBody(umat, field)
>>>
>>> region_stress = fem.RegionQuadBoundary(
... mesh=mesh,
... only_surface=True, # select only faces on the outline
... mask=mesh.points[:, 0] == 1, # select a subset of faces on the surface
... ensure_3d=True, # requires True for axisymmetric/plane strain, otherwise False
... )
>>> field_boundary = fem.FieldContainer([fem.FieldAxisymmetric(region_stress, dim=2)])
>>> stress = fem.SolidBodyCauchyStress(field=field_boundary)
>>>
>>> table = fem.math.linsteps([0, 1], num=5, axis=0, axes=9).reshape(-1, 3, 3)
>>> step = fem.Step(
... items=[solid, stress], ramp={stress: 1 * table}, boundaries=boundaries
... )
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> solid.plot(
... "Principal Values of Cauchy Stress", component=2, clim=[-1.01, -0.99]
... ).show()
"""

def __init__(self, field, cauchy_stress=None):
self.field = field
self._normals = self.field.region.normals

self.results = Results()
self.results.kinematics = self._extract(self.field)

self.results.cauchy_stress = np.zeros((3, 3))
if cauchy_stress is not None:
self.results.cauchy_stress = cauchy_stress

self.assemble = Assemble(
vector=self._vector, matrix=self._matrix, multiplier=-1.0
)
self._area_change = AreaChange()

def update(self, cauchy_stress):
self.__init__(self.field, cauchy_stress)

def _extract(self, field):
self.field = field
self.results.kinematics = self.field.extract()

return self.results.kinematics

def _vector(self, field=None, parallel=False, resize=None):
if field is not None:
self._update(field)
self.results.kinematics = self._extract(self.field)

fun = self._area_change.function(
self.results.kinematics,
self._normals,
parallel=parallel,
)

fun[0] = dot(self.results.cauchy_stress, fun[0], mode=(2, 2), out=fun[0])

self.results.force = IntegralForm(
fun=fun, v=self.field, dV=self.field.region.dV, grad_v=[False]
).assemble(parallel=parallel)

if resize is not None:
self.results.force.resize(*resize.shape)

return self.results.force

def _matrix(self, field=None, parallel=False, resize=None):
if field is not None:
self._update(field)
self.results.kinematics = self._extract(self.field)

fun = self._area_change.gradient(
self.results.kinematics,
self._normals,
parallel=parallel,
)

fun[0] = dot(self.results.cauchy_stress, fun[0], mode=(2, 4), out=fun[0])

self.results.stiffness = IntegralForm(
fun=fun,
v=self.field,
u=self.field,
dV=self.field.region.dV,
grad_v=[False],
grad_u=[True],
).assemble(parallel=parallel)

if resize is not None:
self.results.stiffness.resize(*resize.shape)

return self.results.stiffness

def _update(self, other_field, field=None):
if field is not None:
self.field = field

self.field[0].values = other_field[0].values
self.results.kinematics = self._extract(self.field)

return self.field
40 changes: 40 additions & 0 deletions tests/test_mechanics.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@

import numpy as np
import pytest
from scipy.sparse import csr_matrix

import felupe as fem

Expand Down Expand Up @@ -540,9 +541,48 @@ def test_threefield():
assert np.isclose(job.fnorms[0][-1], 0)


def test_solidbody_cauchy_stress():
field = fem.FieldsMixed(fem.RegionHexahedron(fem.Cube(n=3)), n=3)
region_stress = fem.RegionHexahedronBoundary(
mesh=field.region.mesh,
only_surface=True,
mask=field.region.mesh.x == 1,
)
field_boundary = fem.FieldContainer([fem.Field(region_stress, dim=3)])
boundaries = dict(left=fem.Boundary(field[0], fx=0))
umat = fem.NearlyIncompressible(fem.NeoHooke(mu=1), bulk=5000)
solid = fem.SolidBody(umat, field)

resize_matrix = csr_matrix(([0.0], ([0], [0])), shape=(100, 100))
resize_vector = csr_matrix(([0.0], ([0], [0])), shape=(100, 1))

for cauchy_stress in [None, np.zeros((3, 3))]:
stress = fem.SolidBodyCauchyStress(
field=field_boundary, cauchy_stress=cauchy_stress
)
matrix = stress.assemble.matrix(field, resize=resize_matrix)
vector = stress.assemble.vector(field, resize=resize_vector)
stress._update(other_field=field, field=field_boundary)

assert matrix.shape == (100, 100)
assert vector.shape == (100, 1)

table = (
fem.math.linsteps([1], num=0, axis=2, axes=9)
+ fem.math.linsteps([1], num=0, axis=6, axes=9)
).reshape(-1, 3, 3)
step = fem.Step(
items=[solid, stress], ramp={stress: 1 * table}, boundaries=boundaries
)
fem.Job(steps=[step]).evaluate()

assert np.isclose(field[0].values.max(), 0.971866)


if __name__ == "__main__":
test_simple()
test_solidbody()
test_solidbody_cauchy_stress()
test_solidbody_incompressible()
test_solidbody_axi()
test_solidbody_axi_incompressible()
Expand Down

0 comments on commit 6fe8d85

Please sign in to comment.