Tools#

Optimization

`newtonraphson`([x0, fun, jac, solve, ...])	Find a root of a real function using the Newton-Raphson method.
`tools.IterationState`([iteration, fnorm, ...])	A class to keep track of the state of an iteration during evaluation.
`tools.NewtonResult`(x[, fun, jac, success, ...])	A data class which represents the result found by Newton's method.
`NewtonConvergenceError`	An error which is raised if the Newton-Raphson method fails to converge.

Event Dispatcher

`tools.EventDispatcher`([plugins])	A class to dispatch events to plugins during evaluation.
`tools.Context`([job, step, substep, items, ...])	A class to keep track of the context of a Job during evaluation.

Export of Results

save(region, field[, forces, gradient, ...])

Write field-data to a VTU file.

Convert Cell-Data to Point-Data

`topoints`(values, region[, average, mean])	Shift array of values located at quadrature points of cells to mesh-points.
`project`(values, region[, average, mean, dV, ...])	Project given values at quadrature-points to mesh-points.
`tools.extrapolate`(values, region[, average, ...])	Extrapolate (and optionally average) an array with values at quadrature points to mesh-points.

Reaction-Force and -Moment

`tools.force`(field, forces, boundary)	Evaluate the force vector sum on points of a boundary.
`tools.moment`(field, forces, boundary[, ...])	Evaluate the moment vector sum on points of a boundary at a given center point.

Detailed API Reference

felupe.newtonraphson(x0=None, fun=<function fun>, jac=<function jac>, solve=<function solve>, maxiter=16, update=<function update>, check=<function check>, args=(), kwargs=None, tol=np.float64(1.4901161193847656e-08), items=None, dof1=None, dof0=None, ext0=None, solver=<function spsolve>, verbose=None, callback=None, callback_kwargs=None, progress_bar=None, tqdm='tqdm', dispatcher=None, plugins=None)[source]#

Find a root of a real function using the Newton-Raphson method.

Parameters:

x0 (felupe.FieldContainer, ndarray or None (optional)) – Array or Field container with values of unknowns at a valid starting point (default is None).
fun (callable, optional) – Callable which assembles the vector-valued objective function. Additional args and kwargs are passed. Function signature of a user-defined function has to be fun = lambda x, *args, **kwargs: f.
jac (callable, optional) – Callable which assembles the matrix-valued Jacobian. Additional args and kwargs are passed. Function signature of a user-defined Jacobian has to be jac = lambda x, *args, **kwargs: K.
solve (callable, optional) – Callable which prepares the linear equation system and solves it. If a keyword- argument from the list ["x", "dof1", "dof0", "ext0", "solver"] is found in the function-signature, then these arguments are passed to solve.
maxiter (int, optional) – Maximum number of function iterations (default is 16).
update (callable, optional) – Callable to update the unknowns. Function signature must be update = lambda x0, dx: x.
check (callable, optional) – Callable to check the result with signature check = lambda xnorm, fnorm, success: dx, x, **kwargs.
tol (float, optional) – Tolerance value to check if the function has converged (default is 1.490e-8).
items (list or None, optional) – List with items which provide methods for assembly, e.g. like felupe.SolidBody (default is None).
dof1 (ndarray or None, optional) – 1d-array of int with all active degrees of freedom (default is None).
dof0 (ndarray or None, optional) – 1d-array of int with all prescribed degrees of freedom (default is None).
ext0 (ndarray or None, optional) – Field values at mesh-points for the prescribed components of the unknowns based on dof0 (default is None).
solver (callable, optional) – A sparse or dense solver (default is scipy.sparse.linalg.spsolve()). For a more performant alternative install PyPardiso and use pypardiso.spsolve().
verbose (bool or int or None, optional) – Verbosity level to control how messages are printed during evaluation. If 1 or True and tqdm is installed, a progress bar is shown. If tqdm is missing or verbose is 2, more detailed text-based messages are printed. Default is None. If None, verbosity is set to True. If None and the environmental variable FELUPE_VERBOSE is set and its value is not true, then logging is turned off.
callback (callable or None, optional) – An optional callback function with function signature callback = lambda dx, x, iteration, xnorm, fnorm, success: None, which is called after each completed iteration. Default is None.
progress_bar (tqdm or None, optional) – Use an existing instance of a progress bar if verbose is True. If None and verbose is True, a new bar is created. Default is None.
tqdm (str, optional) – If verbose is True, choose a backend for tqdm ("tqdm", "auto" or "notebook"). Default is "tqdm".
dispatcher (EventDispatcher or None, optional) – An optional EventDispatcher to trigger events during evaluation. Default is None.
plugins (list or None, optional) – A list of plugins with hooks to be used during evaluation, e.g. a LinesearchPlugin. If a dispatcher is given, the plugins are added to a copy of the dispatcher for this evaluation. Default is None.

Returns:

The result object.

Return type:

felupe.tools.NewtonResult

Raises:

NewtonConvergenceError – If the solution of the linearized equation system contains NaN values or if the maximum number of iterations is reached (a subclass of ValueError).

Notes

Nonlinear equilibrium equations \(f(x)\) as a function of the unknowns \(x\) are solved by linearization of \(f\) at a valid starting point of given unknowns \(x_0\), see Eq. (1).

(1)#\[f(x_0) = 0\]

The linearization is given in Eq. (2)

(2)#\[f(x_0 + dx) \approx f(x_0) + K(x_0) \ dx \ (= 0)\]

with the Jacobian as in Eq. (3), evaluated at given unknowns \(x_0\),

(3)#\[K(x_0) = \frac{\partial f}{\partial x}(x_0)\]

and is rearranged to an equation system with left- and right-hand sides, see Eq. (4).

(4)#\[K(x_0) \ dx = -f(x_0)\]

After a solution is found, the unknowns are updated, see Eq. (5).

(5)#\[ \begin{align}\begin{aligned}dx &= \text{solve} \left( K(x_0), -f(x_0) \right) \nonumber\\x &= x_0 + dx\end{aligned}\end{align} \]

Repeated evaluations lead to an incrementally updated solution of \(x\). Herein, \(x_n\) refer to the initial unknowns whereas \(x\) are the updated unknowns, see Eq. (6).

Note

The subscript \((\bullet)_{n+1}\) is dropped for easier readability.

(6)#\[ \begin{align}\begin{aligned}dx &= \text{solve} \left( K(x_n), -f(x_n) \right) \nonumber\\ x &= x_n + dx\end{aligned}\end{align} \]

Then, the nonlinear equilibrium equations are evaluated with the updated unknowns \(f(x)\). The procedure is repeated until convergence is reached.

Plugins may modify the update of the unknowns in the after_linear_solve hook, see IterationState. E.g., a backtracking line search (LinesearchPlugin) scales the Newton increment by a step length \(\alpha \in (0, 1]\), see Eq. (7).

(7)#\[x = x_n + \alpha\ dx\]

Examples

>>> import numpy as np
>>> import felupe as fem
>>>
>>> region = fem.RegionHexahedron(fem.Cube(n=6))
>>> field = fem.FieldContainer([fem.Field(region, dim=3)])
>>> boundaries, loadcase = fem.dof.uniaxial(
...     field, move=0.2, clamped=True, return_loadcase=True
... )
>>> solid = fem.SolidBody(umat=fem.NeoHooke(mu=1.0, bulk=2.0), field=field)
>>> res = fem.newtonraphson(items=[solid], **loadcase, verbose=2)  
 _______  _______  ___ ...
Newton-Raphson solver
=====================

| # | norm(fun) |  norm(dx) |
|---|-----------|-----------|
| 1 | 7.553e-02 | 1.898e+00 |
| 2 | 1.310e-03 | 5.091e-02 |
| 3 | 3.086e-07 | 6.698e-04 |
| 4 | ...e-14 | ...e-07 |

Converged in 4 iterations ...

Newton’s method had success

>>> print(res.success)
True

and 4 iterations were needed to converge within the specified tolerance.

>>> res.iterations
4

The norm of the objective function for all active degrees of freedom is lower than 3e-15.

>>> print(np.linalg.norm(res.fun[loadcase["dof1"]]) < 3e-15)
True

class felupe.tools.IterationState(iteration=None, fnorm=None, xnorm=None, success=None, tol=None, error=None, x=None, dx=None, fun=None, jac=None, alpha=None, updated=False)[source]#

A class to keep track of the state of an iteration during evaluation.

Parameters:

iteration (int or None, optional) – The (zero-based) index of the current iteration (default is None).
fnorm (float or None, optional) – The norm of the objective function of the current iteration (default is None).
xnorm (float or None, optional) – The norm of the increment of the current iteration (default is None).
success (bool or None, optional) – A flag if the current iteration has converged (default is None).
tol (float or None, optional) – The tolerance of the Newton-Raphson method (default is None).
error (bool or None, optional) – A flag if an error occured in the current iteration (default is None).
x (felupe.FieldContainer or ndarray or None, optional) – The unknowns (default is None). Before the update of an iteration, these are the unknowns of the last accepted iteration \(x_n\). After the update, these are the updated unknowns \(x\).
dx (ndarray or None, optional) – The increment of the unknowns of the current iteration (default is None).
fun (ndarray or None, optional) – The values of the objective function (the residuals), evaluated at the unknowns x (default is None).
jac (sparse matrix or ndarray or None, optional) – The Jacobian of the objective function of the current iteration, evaluated at the unknowns of the last accepted iteration \(x_n\) (default is None).
alpha (float or None, optional) – The step length which scales the Newton increment of the current iteration, e.g. of a line search (default is None). None means that the full Newton increment is applied.
updated (bool, optional) – A flag if the unknowns (and the values of the objective function) are already updated in the current iteration (default is False).

Notes

One state is created by newtonraphson() and its attributes are updated in-place during the iterations. Hence, a plugin may keep a reference to the state and attributes set by a plugin persist until they are updated by the Newton-Raphson method. The attributes alpha and updated are reset in the beginning of each iteration.

Attributes of the state in the hooks of an iteration.#
Hook	`x`	`fun`	`jac`	`dx`	`updated`
`before_iteration`	\(x_n\)	\(f(x_n)\)	previous	previous	False
`before_linear_solve`	\(x_n\)	\(f(x_n)\)	\(K(x_n)\)	previous	False
`after_linear_solve`	\(x_n\)	\(f(x_n)\)	\(K(x_n)\)	\(dx\)	False
`after_iteration`	\(x\)	\(f(x)\)	\(K(x_n)\)	\(dx\)	True

A plugin may modify the update of the unknowns in the after_linear_solve hook. Either the increment dx is replaced, e.g. by a scaled increment, which is then used by the Newton-Raphson method to update the unknowns x = update(x_n, dx) and to evaluate the objective function f(x). Or the plugin performs the update on its own, e.g. by the callables update and fun of the Context. Then, the plugin has to set the updated unknowns x, the values of the objective function fun evaluated at the updated unknowns, the applied increment dx and the flag updated=True. The Newton-Raphson method re-uses the values of the objective function and no additional assembly is performed. If the objective function is evaluated for more than one trial, the accepted trial must be evaluated last, because the (temporary) state variables of the items belong to the last evaluation.

Note

The evaluation of the objective function for a list of items links the fields of the items to the evaluated unknowns. If the unknowns of the last accepted iteration \(x_n\) are the field container of an item (e.g. in the first iteration of a job), the value arrays of \(x_n\) are replaced by the value arrays of the evaluated trial. A plugin which evaluates more than one trial must restore the value arrays of \(x_n\) before the next trial is created. If a plugin does not perform the update, the value arrays of \(x_n\) are restored by the Newton-Raphson method.

See also

felupe.Plugin: Base class for plugins.
felupe.tools.Context: A class to keep track of the context of a Job during evaluation.
felupe.LinesearchPlugin: A backtracking line search for the Newton-Raphson method.

class felupe.tools.NewtonResult(x, fun=None, jac=None, success=None, iterations=None, xnorms=None, fnorms=None)[source]#

A data class which represents the result found by Newton’s method. All parameters are available as attributes.

Parameters:

x (felupe.FieldContainer or ndarray) – Array or Field container with values at a solution found by Newton’s method.
fun (ndarray or None, optional) – Values of objective function (default is None).
jac (ndarray or None, optional) – Values of the Jacobian of the objective function (default is None).
success (bool or None, optional) – A boolean flag which is True if the solution converged (default is None).
iterations (int or None, optional) – Number of iterations until solution converged (default is None).
xnorms (array of float or None, optional) – List with norms of the values of the solution (default is None).
fnorms (float or None, optional) – List with norms of the objective function (default is None).

Notes

The objective function’s norm is relative based on the function values on the prescribed degrees of freedom \((\bullet)_0\). A small number \(\varepsilon\) is added to avoid numeric instabilities.

\[\text{norm}(\boldsymbol{f}) = \frac{||\boldsymbol{f}_1||} {\varepsilon + ||\boldsymbol{f}_0||}\]

class felupe.NewtonConvergenceError[source]#

An error which is raised if the Newton-Raphson method fails to converge.

This error is raised by newtonraphson() if the solution of the linearized equation system contains NaN values or if the maximum number of iterations is reached, by a LinesearchPlugin if no acceptable step length is found and by a CutbackPlugin if a failed substep is not recovered.

Notes

It is a subclass of ValueError. Hence, existing code which catches a ValueError is not affected. Plugins like the CutbackPlugin catch this error to recover a failed substep, while other errors, e.g. a ValueError of mismatching shapes of arrays, are raised without a recovery.

Examples

try:
    res = fem.newtonraphson(items=[solid], ...)
except fem.NewtonConvergenceError:
    ...  # e.g. reduce the load increment

class felupe.tools.EventDispatcher(plugins=None)[source]#

A class to dispatch events to plugins during evaluation.

Parameters:: plugins (list or None, optional) – The list of plugins.

add_plugin(plugin)[source]#: Add a plugin to the dispatcher and reconfigure it.

configure(plugins)[source]#: Configure the dispatcher with the given plugins and return a dict with hooks as keys.

trigger(hook, context, state)[source]#: Trigger a hook with context and current state to all registered functions.

class felupe.tools.Context(job=None, step=None, substep=None, items=None, dof1=None, dof0=None, ext0=None, fun=None, update=None, check=None, x0=None, solve=None)[source]#

A class to keep track of the context of a Job during evaluation.

Parameters:

job (felupe.Job or None, optional) – The job object.
step (felupe.Step or None, optional) – The step object.
substep (felupe.tools.NewtonResult or None, optional) – The result of the completed substep. The field container of the substep is available as substep.x.
items (list or None, optional) – The list of items of the Newton-Raphson method or of a step (default is None). Only available in the hooks of newtonraphson() and in the hooks of a substep, see generate_states().
dof1 (ndarray or None, optional) – The active degrees of freedom of the Newton-Raphson method (default is None). Only available in the hooks of newtonraphson().
dof0 (ndarray or None, optional) – The prescribed degrees of freedom of the Newton-Raphson method (default is None). Only available in the hooks of newtonraphson().
ext0 (ndarray or None, optional) – The external values of the prescribed degrees of freedom of the Newton-Raphson method (default is None). Only available in the hooks of newtonraphson().
fun (callable or None, optional) – A callable f = fun(x) which evaluates the objective function (the residuals) for given unknowns x. For a list of items, the vectors of all items are assembled and the fields of the items are linked to x. Additional arguments of the Newton-Raphson method are already bound. Default is None. Only available in the hooks of newtonraphson().
update (callable or None, optional) – The callable x = update(x, dx) of the Newton-Raphson method which updates the unknowns (default is None). Only available in the hooks of newtonraphson().
check (callable or None, optional) – A callable xnorm, fnorm, success = check(dx, x, f) which checks the convergence of the Newton-Raphson method for given unknowns x with the values of the objective function f (default is None). In contrast to the check-argument of newtonraphson(), the state variables of the items are not updated. Only available in the hooks of newtonraphson().
x0 (felupe.FieldContainer or None, optional) – The field container with the unknowns of a step, which is the starting point of the Newton-Raphson method of a substep (default is None). Only available in the hooks of a substep, see generate_states().
solve (callable or None, optional) – A callable res = solve(values) which updates the ramped items (and boundaries) of a step with a dict of values, updates the load case and evaluates the Newton-Raphson method, starting from the unknowns x0. It returns a NewtonResult and errors of the Newton-Raphson method are raised. On success, the unknowns x0 are linked to the result, i.e. the result is the starting point of the next evaluation. Default is None. Only available in the hooks of a substep, see generate_states().

See also

felupe.Plugin: Base class for plugins.
felupe.IterationState: A class to keep track of the state of an iteration.
felupe.SubstepState: A class to keep track of the state of a substep.

felupe.save(region, field, forces=None, gradient=None, filename='result.vtu', cell_data=None, point_data=None)[source]#

Write field-data to a VTU file.

Parameters:

region (Region) – The region to be saved.
field (FieldContainer) – The field container to be saved.
forces (ndarray, optional) – Array with reaction forces to be saved (default is None).
gradient (list of ndarray, optional) – The result of umat.gradient() with the first Piola-Kirchhoff stress tensor as the first list item (default is None).
filename (str, optional) – The filename for the results (default is “result.vtu”).
cell_data (dict or None, optional) – Additional dict with cell data (default is None).
point_data (dict or None, optional) – Additional dict with point data (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.uniaxial(
...     field, clamped=True, move=0.3, return_loadcase=False
... )
>>>
>>> umat = fem.NeoHooke(mu=1)
>>> solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
>>> step = fem.Step(items=[solid], boundaries=boundaries)
>>> job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate()

>>> fem.save(region, field, forces=job.res.fun, gradient=solid.results.stress)

felupe.topoints(values, region, average=True, mean=False)[source]#

Shift array of values located at quadrature points of cells to mesh-points.

Parameters:

values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points q of cells c.
region (Region) – A region used to shift the values to the mesh-points.
average (bool, optional) – A flag to return values averaged at mesh-points if True and the mesh of the region is not already disconnected or to return values on a disconnected mesh with discontinuous values at cell transitions if False (default is True).
mean (bool, optional) – A flag to shift the cell-means of the values to the mesh-points (default is False).

Returns:

Field.values – Array of values shifted to the mesh-points p.

Return type:

ndarray of shape (p, …)

Notes

If the number of quadrature-points per cell is greater than the number of points- per-cell, then the values are trimmed.

felupe.project(values, region, average=True, mean=False, dV=None, solver=<function spsolve>)[source]#

Project given values at quadrature-points to mesh-points.

Parameters:

values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points q of cells c.
region (Region) – A region used to project the values to the mesh-points.
average (bool, optional) – A flag to return values averaged at mesh-points if True and the mesh of the region is not already disconnected or to return values on a disconnected mesh with discontinuous values at cell transitions if False (default is True).
mean (bool, optional) – A flag to extrapolate the cell-means of the values to the mesh-points, see felupe.tools.extrapolate(). If True, dV and solver are ignored. Default is False.
dV (ndarray of shape (q, c) or None, optional) – Differential volumes located at the quadrature points of the cells. If None, the differential volumes are taken from the region (default is None).
solver (callable, optional) – A function for a sparse solver with signature x=solver(A, b). Default is scipy.sparse.linalg.spsolve().

Returns:

Field.values – Array of values projected to the mesh-points p.

Return type:

ndarray of shape (p, …)

Notes

The projection finds \(\hat{\boldsymbol{u}}\), located at mesh-points p, for values \(v\), evaluated at quadrature-points q, such that the variation of the minimization problem in Eq. (8) is solved.

(8)#\[ \begin{align}\begin{aligned}\Pi &= \int_V \frac{1}{2} (u - v) \cdot (u - v)\ dV \quad \rightarrow \quad \min\\\delta \Pi &= \int_V (u - v) \cdot \delta u\ dV = 0\end{aligned}\end{align} \]

This leads to assembled system-matrix and -vector(s) as shown in Eq. (9)

(9)#\[ \begin{align}\begin{aligned}\int_V v\ \cdot \delta u\ dV &\qquad \rightarrow \qquad \hat{\boldsymbol{A}}\\\int_V u\ \cdot \delta u\ dV &\qquad \rightarrow \qquad \hat{\boldsymbol{b}}\end{aligned}\end{align} \]

of an equation system to be solved, see Eq. (10). The right-arrows in Eq. (9) represent the assembly of the integral forms into the system matrix or vectors.

(10)#\[\hat{\boldsymbol{A}}\ \hat{\boldsymbol{u}} = \hat{\boldsymbol{b}}\]

The region, on which the values are projected to, may contain a mesh with a different (e.g. a disconnected) points- but same cells-array. The only requirement on the region is that it must use a compatible quadrature scheme. With the values at mesh-points at hand, new fields may be created to project between fields on different regions. If the array of differential volumes is not given, then the region must be created with Region(grad=True) to include the differential volumes.

Warning

Triangular element formulations require regions with quadrature schemes with higher integration orders, the default choices are not sufficient for projection. For RegionTriangle and RegionTetra a second-order scheme is used.

import felupe as fem

mesh = fem.Rectangle(n=2).triangulate().add_midpoints_edges()
quadrature = fem.TriangleQuadrature(order=5)
fem.RegionQuadraticTriangle(mesh, quadrature=quadrature)

For RegionQuadraticTriangle and RegionQuadraticTetra use a fifth-order scheme.

mesh = fem.Cube(n=2).triangulate().add_midpoints_edges()
quadrature = fem.TetrahedronQuadrature(order=5)
fem.RegionQuadraticTetra(mesh, quadrature=quadrature)

Examples

import felupe as fem

mesh = fem.Rectangle(n=2).triangulate()
region = fem.RegionTriangle(mesh)
field = fem.FieldAxisymmetric(region, dim=2)
values = field.extract()

region2 = region.copy()
region2.reload(quadrature=fem.TriangleQuadrature(order=2))
values_projected = project(values, region2)

felupe.tools.extrapolate(values, region, average=True, mean=False)[source]#

Extrapolate (and optionally average) an array with values at quadrature points to mesh-points.

Parameters:

values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points q of cells c.
region (Region) – A region used to extrapolate the values to the mesh-points.
average (bool, optional) – A flag to return values averaged at mesh-points if True and the mesh of the region is not already disconnected or to return values on a disconnected mesh with discontinuous values at cell transitions if False (default is True).
mean (bool, optional) – A flag to take the cell-means, averaged by the quadrature weights, of the values (default is False).

Returns:

Field.values – Array of values projected to the mesh-points p.

Return type:

ndarray of shape (p, …)

felupe.tools.force(field, forces, boundary)[source]#: Evaluate the force vector sum on points of a boundary.

felupe.tools.moment(field, forces, boundary, centerpoint=array([0., 0., 0.]))[source]#: Evaluate the moment vector sum on points of a boundary at a given center point.