Mechanics#

Solid Bodies

`SolidBody`(umat, field[, statevars])	A SolidBody with methods for the assembly of sparse vectors/matrices.
`SolidBodyNearlyIncompressible`(umat, field, bulk)	A (nearly) incompressible solid body with methods for the assembly of sparse vectors/matrices.
`SolidBodyPressure`(field[, pressure])	A hydrostatic pressure boundary on a solid body.
`SolidBodyGravity`(field[, gravity, density])	A gravity (body) force on a solid body.

State Variables

StateNearlyIncompressible(field)

A State with internal fields for (nearly) incompressible solid bodies.

Steps and Jobs

`Step`(items[, ramp, boundaries])	A Step with multiple substeps, subsequently depending on the solution of the previous substep.
`Job`(steps[, callback])	A job with a list of steps and a method to evaluate them.
`CharacteristicCurve`(steps, boundary[, ...])	A job with a list of steps and a method to evaluate them.

Point Load and Multi-Point Constraints

`PointLoad`(field, points[, values, apply_on, ...])	A point load with methods for the assembly of sparse vectors/matrices.
`MultiPointConstraint`(field, points, centerpoint)
`MultiPointContact`(field, points, centerpoint)

Detailed API Reference

class felupe.SolidBody(umat, field, statevars=None)[source]#

Bases: Solid

A SolidBody with methods for the assembly of sparse vectors/matrices.

Parameters:

umat (class) – A class which provides methods for evaluating the gradient and the hessian of the strain energy density function per unit undeformed volume. The function signatures must be dψdF, ζ_new = umat.gradient([F, ζ]) for the gradient and d2ψdFdF = umat.hessian([F, ζ]) for the hessian of the strain energy density function \(\psi(\boldsymbol{F})\), where \(\boldsymbol{F}\) is the deformation gradient tensor and \(\zeta\) holds the array of internal state variables.
field (FieldContainer) – A field container with one or more fields.
statevars (ndarray or None, optional) – Array of initial internal state variables (default is None).

Notes

\[ \begin{align}\begin{aligned}W_{int} &= -\int_V \psi(\boldsymbol{F}) \ dV\\\delta W_{int} &= -\int_V \delta \boldsymbol{F} : \frac{\partial \psi}{\partial \boldsymbol{F}}\ dV\\\Delta\delta W_{int} &= -\int_V \delta \boldsymbol{F} : \frac{\partial^2 \psi}{\partial \boldsymbol{F}\ \partial \boldsymbol{F}} : \Delta \boldsymbol{F} \ dV\end{aligned}\end{align} \]

Examples

>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> region = fem.RegionHexahedron(mesh)
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
>>>
>>> umat = fem.NeoHooke(mu=1, bulk=2)
>>> solid = fem.SolidBody(umat, field)
>>>
>>> table = fem.math.linsteps([0, 1], num=5)
>>> step = fem.Step(
>>>     items=[solid],
>>>     ramp={boundaries["move"]: table},
>>>     boundaries=boundaries,
>>> )
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress")

class felupe.SolidBodyNearlyIncompressible(umat, field, bulk, state=None, statevars=None)[source]#

Bases: Solid

A (nearly) incompressible solid body with methods for the assembly of sparse vectors/matrices. The constitutive material definition must provide the distortional part of the strain energy density function per unit undeformed volume only. The volumetric part of the strain energy density function is automatically added.

Parameters:

umat (class) – A class which provides methods for evaluating the gradient and the hessian of the isochoric part of the strain energy density function per unit undeformed volume \(\hat\psi(\boldsymbol{F})\). The function signatures must be dψdF, ζ_new = umat.gradient([F, ζ]) for the gradient and d2ψdFdF = umat.hessian([F, ζ]) for the hessian of the strain energy density function \(\psi(\boldsymbol{F})\), where \(\boldsymbol{F}\) is the deformation gradient tensor and \(\zeta\) holds the array of internal state variables.
field (FieldContainer) – A field container with one or more fields.
bulk (float) – The bulk modulus of the volumetric material behaviour (\(U(\bar{J})=K(\bar{J}-1)^2/2\)).
state (StateNearlyIncompressible or None, optional) – 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).

Notes

\[ \begin{align}\begin{aligned}W_{int} &= -\int_V \hat\psi(\boldsymbol{F}) \ dV + \int_V U(\bar{J}) \ dV + \int_V p (J - \bar{J}) \ dV\\\hat{W}_{int} &= -\int_V \hat\psi(\boldsymbol{F}) \ dV\\\delta \hat{W}_{int} &= -\int_V \delta \boldsymbol{F} : \frac{\partial \hat\psi}{\partial \boldsymbol{F}} \ dV\\\Delta\delta \hat{W}_{int} &= -\int_V \delta \boldsymbol{F} : \frac{\partial^2 \hat\psi}{\partial\boldsymbol{F}\ \partial\boldsymbol{F}} : \Delta \boldsymbol{F} \ dV\end{aligned}\end{align} \]

The volumetric material behaviour is hard-coded and is defined by the strain energy function.

\[U(\bar{J}) = \frac{K}{2} (\bar{J} - 1)^2\]

Hu-Washizu Three-Field-Variational Principle

The Three-Field-Variation \((\boldsymbol{u},p,\bar{J})\) leads to a linearized equation system with nine sub block-matrices. Due to the fact that the equation system is derived by a potential, the matrix is symmetric and hence, only six independent sub-matrices have to be evaluated. Furthermore, by the application of the mean dilatation technique, two of the remaining six sub-matrices are identified to be zero. That means four sub-matrices are left to be evaluated, where two non-zero sub-matrices are scalar-valued entries.

\[\begin{split}\begin{bmatrix} \boldsymbol{A} & \boldsymbol{b} & \boldsymbol{0} \\ \boldsymbol{b}^T & 0 & -c \\ \boldsymbol{0}^T & -c & d \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{x} \\ y \\ z \end{bmatrix} = \begin{bmatrix} \boldsymbol{u} \\ v \\ w \end{bmatrix}\end{split}\]

An alternative representation of the equation system, only dependent on the primary unknowns \(\boldsymbol{u}\) is carried out. To do so, the second line is multiplied by \(\frac{d}{c}\).

\[\begin{split}\begin{bmatrix} \boldsymbol{A} & \boldsymbol{b} & \boldsymbol{0} \\ \frac{d}{c}~\boldsymbol{b}^T & 0 & -c \\ \boldsymbol{0}^T & -c & d \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{x} \\ y \\ z \end{bmatrix} = \begin{bmatrix} \boldsymbol{u} \\ \frac{d}{c}~v \\ -w \end{bmatrix}\end{split}\]

Now, equations two and three are summed up. This eliminates one of the three unknowns.

\[\begin{split}\begin{bmatrix} \boldsymbol{A} & \boldsymbol{b} \\ \frac{d}{c}~\boldsymbol{b}^T & -c \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{x} \\ y \end{bmatrix} = \begin{bmatrix} \boldsymbol{u} \\ \frac{d}{c}~v + w \end{bmatrix}\end{split}\]

Next, the second equation is left-multiplied by \(\frac{1}{c}~\boldsymbol{b}\) and both equations are summed up again.

\[\begin{bmatrix} \boldsymbol{A} + \frac{d}{c^2}~\boldsymbol{b} \otimes \boldsymbol{b} \end{bmatrix} \cdot \begin{bmatrix} \boldsymbol{x} \end{bmatrix} = \begin{bmatrix} \boldsymbol{u} + \frac{d}{c^2}~\boldsymbol{b}~v + \frac{1}{c}~\boldsymbol{b}~w \end{bmatrix}\]

The secondary unknowns are evaluated after solving the primary unknowns.

\[ \begin{align}\begin{aligned}z &= \frac{1}{c}~\boldsymbol{b}^T \boldsymbol{x} - \frac{1}{c}~v\\y &= \frac{d}{c}~z - \frac{1}{c}~w\end{aligned}\end{align} \]

For the mean-dilatation technique, the variables, equations as well as sub-matrices are evaluated. Note that the pairs of indices \((ai)\) and \((bk)\) have to be treated as 1d-vectors.

\[ \begin{align}\begin{aligned}A_{aibk} &= \int_V \frac{\partial h_a}{\partial X_J} \left( \frac{\partial^2 \overset{\wedge}{\psi}}{\partial F_{iJ} \partial F_{kL}} + p \frac{\partial^2 J}{\partial F_{iJ} \partial F_{kL}} \right) \frac{\partial h_b}{\partial X_L} \ dV\\b_{ai} &= \int_V \frac{\partial h_a}{\partial X_J} \frac{\partial J}{\partial F_{iJ}} \ dV\\c &= \int_V \ dV = V\\d &= \int_V \frac{\partial^2 U(\bar{J})}{\partial \bar{J} \partial \bar{J}} \ dV = \bar{U}'' V\end{aligned}\end{align} \]

and

\[ \begin{align}\begin{aligned}x_{ai} &= \delta {u}_{ai}\\y &= \delta p\\z &= \delta \bar{J}\end{aligned}\end{align} \]

as well as

\[ \begin{align}\begin{aligned}u_{ai} (= -r_{ai}) &= -\int_V \frac{\partial h_a}{\partial X_J} \left( \frac{\partial \overset{\wedge}{\psi}}{\partial F_{iJ}} + p \frac{\partial J}{\partial F_{iJ}} \right) \ dV\\v &= -\int_V (J - \bar{J}) \ dV = \bar{J} V - v\\w &= -\int_V (\bar{U}' - p) \ dV = p V - \bar{U}' V\end{aligned}\end{align} \]

Examples

>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> region = fem.RegionHexahedron(mesh)
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>> boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
>>>
>>> umat = fem.NeoHooke(mu=1)
>>> solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
>>>
>>> table = fem.math.linsteps([0, 1], num=5)
>>> step = fem.Step(
>>>     items=[solid],
>>>     ramp={boundaries["move"]: table},
>>>     boundaries=boundaries,
>>> )
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress")

class felupe.StateNearlyIncompressible(field)[source]#

Bases: object

A State with internal fields for (nearly) incompressible solid bodies.

h(parallel=False)[source]#: Integrated shape-function gradient w.r.t. the deformed coordinates.

\[\int_V \frac{\partial J}{\partial \boldsymbol{F}} : \delta \boldsymbol{F} ~ dV\]

v()[source]#: Cell volumes of the deformed configuration.

\[v = \int_V J ~ dV\]

class felupe.SolidBodyGravity(field, gravity=None, density=1.0)[source]#

Bases: object

A gravity (body) force on a solid body.

Parameters:

field (FieldContainer) – A field container with fields created on a boundary region.
gravity (ndarray or None, optional) – The prescribed values of gravity \(\boldsymbol{g}\) (default is None). If None, the gravity vector is set to zero (the dimension of the gravity vector is derived from the first field of the field container).
density (float, optional) – The density \(\rho\) of the solid body (default is 1.0).

Notes

\[\delta W_{ext} = \int_V \delta \boldsymbol{u} \cdot \rho \boldsymbol{g} \ dV\]

Examples

>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> region = fem.RegionHexahedron(mesh)
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>> boundaries = fem.dof.symmetry(field[0])
>>>
>>> umat = fem.NeoHooke(mu=1, bulk=2)
>>> solid = fem.SolidBody(umat, field)
>>> gravity = fem.SolidBodyGravity(field, density=1.0)
>>>
>>> table = fem.math.linsteps([0, 1], num=5, axis=0, axes=3)
>>> step = fem.Step(
>>>     items=[solid, gravity],
>>>     ramp={gravity: 2 * table},
>>>     boundaries=boundaries,
>>> )
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress")

update(gravity)[source]#

class felupe.SolidBodyPressure(field, pressure=None)[source]#

Bases: object

A hydrostatic pressure boundary on a solid body.

Parameters:

field (FieldContainer) – A field container with fields created on a boundary region.
pressure (float or ndarray or None, optional) – A scaling factor for the prescribed pressure \(p\) (default is None). If None, the pressure is set to zero.

Notes

\[\delta W_{ext} = \int_{\partial V} \delta \boldsymbol{u} \cdot p J \boldsymbol{F}^{-T} \ d\boldsymbol{A}\]

Examples

>>> 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_pressure = 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_pressure, dim=2)])
>>> pressure = fem.SolidBodyPressure(field=field_boundary)
>>>
>>> table = fem.math.linsteps([0, 1], num=5)
>>> step = fem.Step(
>>>     items=[solid, pressure], ramp={pressure: 1 * table}, boundaries=boundaries
>>> )
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress", component=2)

update(pressure)[source]#

class felupe.PointLoad(field, points, values=None, apply_on=0, axisymmetric=False)[source]#

Bases: object

A point load with methods for the assembly of sparse vectors/matrices.

update(values)[source]#

class felupe.Step(items, ramp=None, boundaries=None)[source]#

Bases: object

A Step with multiple substeps, subsequently depending on the solution of the previous substep.

Parameters:

items (list of SolidBody, SolidBodyNearlyIncompressible, SolidBodyPressure, SolidBodyGravity, PointLoad, MultiPointConstraint or MultiPointContact) – A list of items with methods for the assembly of sparse vectors/matrices.
ramp (dict, optional) – A dict with Boundary or item-keys which holds the array of values to ramp (default is None). If None, only one substep is evaluated.
boundaries (dict of Boundary, optional) – A dict with Boundary conditions (default is None).

Examples

>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> region = fem.RegionHexahedron(mesh)
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>>
>>> boundaries = fem.dof.symmetry(field[0])
>>> boundaries["clamped"] = fem.Boundary(field[0], fx=1, skip=(True, False, False))
>>> boundaries["move"] = fem.Boundary(field[0], fx=1, skip=(False, True, True))
>>>
>>> umat = fem.NeoHooke(mu=1, bulk=2)
>>> solid = fem.SolidBody(umat, field)
>>>
>>> move = fem.math.linsteps([0, 1], num=5)
>>> step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress")

generate(**kwargs)[source]#: Yield all generated substeps.

class felupe.Job(steps, callback=<function Job.<lambda>>)[source]#

Bases: object

A job with a list of steps and a method to evaluate them.

Parameters:

steps (list of Step) – A list with steps, where each step subsequently depends on the solution of the previous step.
callback (callable, optional) – A callable which is called after each completed substep. Function signature must be lambda stepnumber, substepnumber, substep: None, where substep is an instance of NewtonResult. THe field container of the completed substep is available as substep.x. Default callback is lambda stepnumber, substepnumber, substep: None.

steps#

A list with steps, where each step subsequently depends on the solution of the previous step.

Type:: list of Step

nsteps#

The number of steps.

Type:: int

callback#

A callable which is called after each completed substep. Function signature must be lambda stepnumber, substepnumber, substep: None, where substep is an instance of NewtonResult. THe field container of the completed substep is available as substep.x. Default callback is lambda stepnumber, substepnumber, substep: None.

Type:: callable, optional

timetrack#

A list with times at which the results are written to the XDMF result file.

Type:: list of int

fnorms#

List with norms of the objective function for each completed substep of each step. See also class:~felupe.tools.NewtonResult.

Type:: list of list of float

Examples

>>> import felupe as fem
>>>
>>> mesh = fem.Cube(n=6)
>>> region = fem.RegionHexahedron(mesh)
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>>
>>> boundaries = fem.dof.symmetry(field[0])
>>> boundaries["clamped"] = fem.Boundary(field[0], fx=1, skip=(True, False, False))
>>> boundaries["move"] = fem.Boundary(field[0], fx=1, skip=(False, True, True))
>>>
>>> umat = fem.NeoHooke(mu=1, bulk=2)
>>> solid = fem.SolidBody(umat, field)
>>>
>>> move = fem.math.linsteps([0, 1], num=5)
>>> step = fem.Step(items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries)
>>>
>>> job = fem.Job(steps=[step]).evaluate()
>>> ax = solid.imshow("Principal Values of Cauchy Stress")