Tools#
Optimization

Find a root of a real function using the NewtonRaphson method. 

A data class which represents the result found by Newton's method. 
Export of Results

Write fielddata to a VTU file. 
Convert CellData to PointData

Shift array of values located at quadrature points of cells to meshpoints. 

Project given values at quadraturepoints to meshpoints. 

Extrapolate (and optionally average) an array with values at quadrature points to meshpoints. 
ReactionForce and Moment

Evaluate the force vector sum on points of a 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=np.float64(1.4901161193847656e08), items=None, dof1=None, dof0=None, ext0=None, solver=<function spsolve>, verbose=None)[source]#
Find a root of a real function using the NewtonRaphson 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 vectorvalued objective function. Additional args and kwargs are passed. Function signature of a userdefined function has to be
fun = lambda x, *args, **kwargs: f
.jac (callable, optional) – Callable which assembles the matrixvalued Jacobian. Additional args and kwargs are passed. Function signature of a userdefined 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 functionsignature, then these arguments are passed tosolve
.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.490e8).
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) – 1darray of int with all active degrees of freedom (default is None).
dof0 (ndarray or None, optional) – 1darray of int with all prescribed degress of freedom (default is None).
ext0 (ndarray or None, optional) – Field values at meshpoints 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 usepypardiso.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. Iftqdm
is missing or verbose is 2, more detailed textbased 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 nottrue
, then logging is turned off.
 Returns:
The result object.
 Return type:
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 righthand 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 inital 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.
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) NewtonRhapson solver =====================  #  norm(fun)  norm(dx)    1  7.553e02  1.898e+00   2  1.310e03  5.091e02   3  3.086e07  6.698e04   4  2.255e14  1.527e07  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 3e15.
>>> np.linalg.norm(res.fun[loadcase["dof1"]]) 2.7384964752762237e15
 felupe.save(region, field, forces=None, gradient=None, filename='result.vtu', cell_data=None, point_data=None)[source]#
Write fielddata 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 PiolaKirchhoff 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 meshpoints.
 Parameters:
values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points
q
of cellsc
.region (Region) – A region used to shift the values to the meshpoints.
average (bool, optional) – A flag to return values averaged at meshpoints 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 cellmeans of the values to the meshpoints (default is False).
 Returns:
Field.values – Array of values shifted to the meshpoints
p
. Return type:
ndarray of shape (p, …)
Notes
If the number of quadraturepoints per cell is greater than the number of points percell, then the values are trimmed.
 felupe.project(values, region, average=True, mean=False, dV=None, solver=<function spsolve>)[source]#
Project given values at quadraturepoints to meshpoints.
 Parameters:
values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points
q
of cellsc
.region (Region) – A region used to project the values to the meshpoints.
average (bool, optional) – A flag to return values averaged at meshpoints 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 cellmeans of the values to the meshpoints, see
felupe.tools.extrapolate()
. If True,dV
andsolver
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 isscipy.sparse.linalg.spsolve()
.
 Returns:
Field.values – Array of values projected to the meshpoints
p
. Return type:
ndarray of shape (p, …)
Notes
The projection finds \(\hat{\boldsymbol{u}}\), located at meshpoints
p
, for values \(v\), evaluated at quadraturepointsq
, such that the variation of the minimization problem in Eq. (7) is solved.(7)#\[ \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 systemmatrix and vector(s) as shown in Eq. (8)
(8)#\[ \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. (9). The rightarrows in Eq. (8) represent the assembly of the integral forms into the system matrix or vectors.
(9)#\[\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 cellsarray. The only requirement on the region is that it must use a compatible quadrature scheme. With the values at meshpoints 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
andRegionTetra
a secondorder 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
andRegionQuadraticTetra
use a fifthorder 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 meshpoints.
 Parameters:
values (ndarray of shape (..., q, c)) – Array with values located at the quadrature points
q
of cellsc
.region (Region) – A region used to extrapolate the values to the meshpoints.
average (bool, optional) – A flag to return values averaged at meshpoints 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 cellmeans, averaged by the quadrature weights, of the values (default is False).
 Returns:
Field.values – Array of values projected to the meshpoints
p
. Return type:
ndarray of shape (p, …)