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[, sym, mode])	Shift of scalar or tensorial values at quadrature points to mesh-points.
`project`(values, region[, average, mean, dV, ...])	Project given values at quadrature-points projected 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=True)[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, optional) – Verbosity level: False or 0 for no logging, True or 1 for a progress bar and 2 for a text-based logging output (default is True).If the environmental variable FELUPE_VERBOSE is set and its value is false, then this argument is ignored and 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’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.7482611016095555e-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"])
>>> job.evaluate()

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

felupe.topoints(values, region, sym=True, mode='tensor')[source]#: Shift of scalar or tensorial values at quadrature points to mesh-points.

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

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

valuesndarray of shape (…, q, c)
Array with values located at the quadrature points q of cells c.

regionRegion
A region used to project the values to the mesh-points.

averagebool, 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).

meanbool, optional
A flag to take the cell-means of the values. Quadrature weights are ignored. Default is False.

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

solvercallable, optional
A function for a sparse solver with signature x=solver(A, b). Default is scipy.sparse.linalg.spsolve().

Field.valuesndarray of shape (p, …)
Array of values projected to the mesh-points p.

The projection finds \(\hat{oldsymbol{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)#\[\Pi &= \int_V\]

rac{1}{2} (v - u)^2 delta udV quad ightarrow quad min

delta Pi &= int_V (v - u) delta udV = 0

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

(3)#\[&\int_V v\ \delta u\ dV \qquad\]

ightarrow qquad hat{oldsymbol{A}}

&int_V udelta udV qquad

ightarrow qquad hat{oldsymbol{b}}

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

(4)#\[\hat{oldsymbol{A}}\ \hat{oldsymbol{u}} = \hat{oldsymbol{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)
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.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.