Tools#
- 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 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.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 usepypardiso.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).
- 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\).
\[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)[source]#
Projection (and optionally averaging) of scalar or vectorial values at quadrature points to mesh-points.