Tools#

Optimization

`newtonrhapson`([x0, fun, jac, solve, ...])	Find a root of a real function using the Newton-Raphson method.
`tools.NewtonResult`(x[, fun, jac, success, ...])	A data class which represents the result found by Newton's method.

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

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||}\]

felupe.newtonrhapson(x0=None, fun=<function fun>, jac=<function jac>, solve=<function solve>, maxiter=16, update=<function update>, check=<function check>, args=(), kwargs={}, tol=1.4901161193847656e-08, items=None, dof1=None, dof0=None, ext0=None, solver=<function spsolve>, verbose=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 the check the result.
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 degress 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.

Returns:

The result object.

Return type:

felupe.tools.NewtonResult

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\).

\[f(x_0) = 0\]

The linearization is given by

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

with the Jacobian, evaluated at given unknowns \(x_0\),

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

and is rearranged to an equation system with left- and right-hand sides.

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

After a solution is found, the unknowns are updated.

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

Repeated evaluations lead to an incrementally updated solution of \(x\). Herein \(x_n\) refer to the inital unknowns whereas \(x\) are the updated unknowns (the subscript \((\bullet)_{n+1}\) is dropped for readability).

(1)#\[ \begin{align}\begin{aligned}dx &= \text{solve} \left( K(x_n), -f(x_n) \right)\\ 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.

Examples

>>> 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)
>>> solid = fem.SolidBody(umat=fem.NeoHooke(mu=1.0, bulk=2.0), field=field)
>>> res = fem.newtonrhapson(items=[solid], **loadcase)  

Newton-Rhapson 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 | 2.255e-14 | 1.527e-07 |

Converged in 4 iterations ...

Newton’s method had success

>>> 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.

>>> np.linalg.norm(res.fun[loadcase["dof1"]])
2.7384964752762237e-15

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, loadcase = fem.dof.uniaxial(field, clamped=True, move=0.3)
>>>
>>> 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. (2) is solved.

(2)#\[ \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. (3)

(3)#\[ \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. (4). The right-arrows in Eq. (3) represent the assembly of the integral forms into the system matrix or vectors.

(4)#\[\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.