v11.2.0

[11.2.0] - 2026-10-13

Added

• Add mesh.extrude(mesh, path, other_mesh=None) and Mesh.extrude(path, other_mesh=None) to extrude a line- or quad-mesh along a given path. The section is transported by rotation-minimizing frames and its initial orientation w.r.t. the tangent of the path is preserved. Optionally, the section is interpolated between two meshes.
• Add constitution.ThirdMediumContactMixed, a mixed-formulation for a stabilized third medium contact.
• Add LinesearchPlugin(admissible=True, residual=True, c=1e-4, max_halvings=10, criteria=None), a backtracking line search for the Newton-Raphson method, which is enabled by Job(steps, plugins=[LinesearchPlugin()]) or newtonraphson(..., plugins=[LinesearchPlugin()]). The Newton increment is scaled by a step length alpha (starting with one) which is halved until all acceptance criteria are fulfilled: positive determinants of the deformation gradients of all items (evaluated first, without assembly) and a sufficient decrease of the norm of the residuals (Armijo-type, the change of the prescribed degrees of freedom is included by its linearization). Custom criteria criterion(trial) -> bool get a plugins.LinesearchTrial. The residuals of the accepted trial are re-used and a NewtonConvergenceError is raised if the maximum number of halvings is exceeded. The accepted step lengths are available in LinesearchPlugin.alphas and IterationState.alpha, and reduced step lengths are shown by the ProgressPlugin.
• Add the keyword-argument plugins=None to newtonraphson(), to use plugins with hooks for the Newton-Raphson method without a job.
• Add the attributes x, dx, fun, jac, alpha and updated to IterationState and the attributes items, dof1, dof0, ext0, fun, update and check to Context in the hooks of newtonraphson(). In the after_linear_solve hook, plugins may replace the Newton increment or perform the update of the unknowns on their own, see IterationState.
• Add CutbackPlugin(factor=0.5, max_cutbacks=10, growth=2.0, exceptions=(NewtonConvergenceError, ArithmeticError, numpy.linalg.LinAlgError)), a cutback of the increment of failed substeps, which is enabled by Job(steps, plugins=[CutbackPlugin()]). If the Newton-Raphson method of a substep fails, the unknowns and the items are restored to a checkpoint and the substep is subdivided into increments with linearly interpolated values of the ramp. The increment is multiplied by factor after each failed attempt and by growth after each converged increment (limited by the remaining part of the substep). The increment must not become smaller than factor**max_cutbacks, i.e. only reductions below the smallest increment reached so far count towards max_cutbacks (failed attempts with increased increments do not count). The states of the items, e.g. the internal fields of a SolidBodyNearlyIncompressible, are restored by their methods checkpoint() and restore() and the state variables of the materials are only updated for converged increments. The load factors of the converged increments and the number of cutbacks are available in CutbackPlugin.load_factors and CutbackPlugin.cutbacks, and recovered substeps are shown by the ProgressPlugin. A NewtonConvergenceError, caused by the error of the last attempt, is raised if the increment would become smaller than the smallest increment. A cutback may be combined with a LinesearchPlugin.
• Add NewtonConvergenceError, a subclass of ValueError, which is raised by newtonraphson() if the solution contains NaN values or if the maximum number of iterations is reached, by the LinesearchPlugin if no acceptable step length is found and by the CutbackPlugin if a failed substep is not recovered. Errors of other types, e.g. a ValueError of mismatching shapes of arrays, do not lead to a cutback by the CutbackPlugin with its default exceptions. The sparse solvers scipy.sparse.linalg.spsolve() and pypardiso.spsolve() do not raise an error for a singular matrix (the former returns NaN values, the latter a perturbed solution), hence no further error types of the solvers are included.
• Add the hooks before_substep and after_failed_substep to Plugin, which are triggered in Step.generate_states(). In the after_failed_substep hook, plugins may recover a failed substep, see SubstepState.
• Add SubstepState, a JobState with the additional attributes values, result and load_factors. One state is created for each substep and it is passed to the hooks before_substep, after_failed_substep and after_substep. Add the attributes x0 and solve to Context in the hooks of a substep. The callable solve(values) updates the ramped items with given values, evaluates the Newton-Raphson method and links the unknowns x0 to the result.
• Add Step.generate_states(stepnumber=None, time=None, **kwargs), which yields the SubstepState of each completed substep. Step.generate() still yields the results of the substeps.
• Add the attribute error to JobState.
• Add the methods checkpoint() and restore(checkpoint) to thermal.TimeStep (the old time) and to ContactRigidPlane (the state of the frictional contact).
• Add view.select_surface_points() and Mesh.select_surface_points() to interactively select points by surface patches. A left click cycles the patch under the cursor from unselected to selected, to excluded (excluded_color="darkred") and back to unselected. Excluded patches have priority, i.e. the points shared by selected and excluded patches are not selected. Patches may be selected or excluded initially by the coordinates of points next to them, select_surface_points(selected=None, excluded=None), same as left clicks on these patches. The window shows the axes and uses the default view of Mesh.plot(). The surface of quadratic (or Lagrange) cells is subdivided by select_surface_points(nonlinear_subdivision=1), same as in Mesh.plot(), i.e. it is plotted with smooth shading and curved edges of the faces, and the patches are evaluated on the subdivided surface. The returned point ids are always those of the faces of the mesh, including the mid-edge points. The edges of the mesh are hidden by default (show_edges=False), the patches are separated by their borders.
• Add view.select_edge_points() and Mesh.select_edge_points() to interactively select points by edge patches. The edges are the borders of the surface patches, i.e. the edges between faces with a kink angle above angle and the boundary edges of open surfaces, e.g. of 2d-meshes. An edge patch is a chain of connected edges, which is split at points with a kink angle above angle and at points with more than two edges. Edges which are hidden behind the surface are not selected. Edge patches are selected and excluded (excluded_color="red", initially by selected=None, excluded=None) like the surface patches of view.select_surface_points(). The edges of the mesh on the surface are hidden by default (show_edges=False).
• Add Boundary(..., select_points=None) to interactively select points for the boundary condition by surface patches (select_points="surfaces") or by edge patches (select_points="edges").
• Add interactive selection methods for a boundary condition, Boundary.select_surface_points() and Boundary.select_edge_points() to apply a new mask.
• Add CylindricalConstraint(field, points, center=(0, 0, 0), axis=(0, 0, 1), displacement=(0, 0, 0), skip=(False, False, False), multiplier=1e3), a penalty-based constraint of the displacements of points in cylindrical coordinates (radial, circumferential and axial components w.r.t. a given axis), e.g. to fix the outer surface of a cylinder in radial direction only. The constraint is geometrically exact: the circumferential displacement is the arc length on the undeformed radius, i.e. a constrained circumferential direction rotates the radial and circumferential directions by the prescribed angle, and a radially constrained point with a skipped circumferential direction slides on the cylinder (like in a rigid sleeve). The prescribed displacements may be ramped in a Step by CylindricalConstraint.update(displacement). 2d-meshes are supported in polar coordinates.
• Add CartesianConstraint(field, points, rotation=None, displacement=(0, 0, 0), skip=(False, False, False), multiplier=1e3), a penalty-based constraint of the displacements of points in a fixed Cartesian coordinate system, given by the rotation matrix of the global base vectors (like math.rotation_matrix()) with the (orthogonal and right-handed) base vectors as columns, e.g. for supports on inclined planes where the points slide in a plane which is not aligned to the global axes. The gaps are linear in the displacements, i.e. the stiffness matrix is constant. The prescribed displacements may be ramped in a Step by CartesianConstraint.update(displacement). The CartesianConstraint and the CylindricalConstraint share a common base class for the selection of the points, the skipped directions, the prescribed displacements, the assembly of the penalty potential and the plot of the constrained points.
• Add rotations to the MultiPointConstraint(field, points, centerpoint, skip=(False, False, False), multiplier=1e3, rotationpoint=None). The values of an additional rotation-point are the rotations (in rad) of the connected points about the center-point, i.e. the points are moved rigidly by the displacements of the center-point and by a rotation, which may be prescribed or free, e.g. for a twisted end face of a beam or for the rigid core of a coated sphere which rolls on a plane. Without a rotation-point, the rotations are locked (unchanged). The rotation matrix is parametrized incrementally, R = exp([ψ - ψ_n]×) R_n, with the rotation-vector ψ_n and the rotation matrix R_n of the last converged substep as state variables. Hence, the parametrization is not singular for (accumulated) rotations of 2π and the values of the rotation-point accumulate the increments of the rotation, i.e. a rotation about a fixed axis may be ramped beyond one revolution, while the rotation of a path with varying axes is path-dependent. The residuals and the symmetric stiffness matrix are consistent, with the derivatives of Rodrigues' formula in closed form (power series for small angles). The rotation matrix and the moment w.r.t. the deformed center-point are available as results.rotation and results.moment. For 2d-meshes, the first value of the rotation-point is the rotation about the z-axis. The rotation-point may be located at the center-point. Skipped axes are not supported with a rotation-point. The positional arguments of the MultiPointConstraint are unchanged.
• Add the methods checkpoint() and restore(checkpoint) to MultiPointConstraint (the state variables of the rotation), which are used by the CutbackPlugin.
• Add the keyword-arguments multiplier_update="increment" and symmetric=False to ContactRigidPlane. The stiffness-based multipliers (with items) are evaluated once per increment, based on the stiffness of the last converged increment, and kept constant during the Newton-iterations (multiplier_update="increment"). Before, they were re-evaluated in each assembly without being linearized, which slowed down or prevented the convergence of the Newton-Raphson method. The former behaviour is available by multiplier_update="iteration" and constant multipliers, evaluated only on the first evaluation, by multiplier_update="once". With symmetric=True, the symmetric part of the (non-symmetric) stiffness matrix of frictional contact is returned, e.g. for symmetric solvers.
• Add the keyword-argument elastic_slip=None to ContactRigidPlane, the allowable elastic slip of sticking points. If given, the tangential multiplier is proportional to the normal force, λ_t = μ f_n / elastic_slip, i.e. a sticking point starts to slide at an elastic slip of elastic_slip, independent of the contact pressure (similar to the penalty friction of Abaqus, Ansys and Marc). The scale factor multiplier_tangential is not used and the derivative of the tangential multiplier is included in the stiffness matrix. This improves the convergence at the transition from partial stick to gross sliding, e.g. for soft (rubber-like) bodies with large tangential displacements.
• Add ContactRigidCylinder(field, points, centerpoint, radius, axis=None, inside=False, friction=0.0, multiplier=1.0, multiplier_tangential=0.1, items=None, multiplier_update="increment", symmetric=False, elastic_slip=None, curvature="converged"), a node-to-surface contact with a rigid cylinder (a rigid circle for 2d-meshes), e.g. for cylindrical indenters or rollers (inside=False) and rigid tubes (inside=True). The cylinder is connected to a center-point on its axis and does not rotate. The contact formulation is the same as for a ContactRigidPlane, but the gaps and the outward normal vectors depend on the positions of the points. The elastic slip is measured in the current tangent plane and the stiffness matrix includes the derivatives of the normal vectors (curvatures) in the normal and in the tangential contributions. The curvature term of the normal stiffness is evaluated with the gaps of the last converged increment by default (curvature="converged"), i.e. bounded by the converged normal forces, because the gaps of the trial state of large increments may lead to an indefinite stiffness matrix. The consistent linearization is available by curvature="trial". The cylinder is plotted by ContactRigidCylinder.plot().
• Add mechanics.ContactRigidBase, the base class of ContactRigidPlane and ContactRigidCylinder with the normal contact, the Coulomb friction (return mapping from the converged state), the stiffness-based multipliers, the checkpoints and the assembly. A subclass defines the geometry of the rigid surface by the abstract method _kinematics(dx), which returns the gaps, the outward normal vectors and their derivatives w.r.t. the gap vectors (or None for flat surfaces).
• Add the evaluation of constitutive materials in chunks of cells. ConstitutiveMaterial.gradient(x, *args, out=None, **kwargs) and ConstitutiveMaterial.hessian(x, *args, out=None, **kwargs) evaluate the methods _gradient(x) and _hessian(x) of a derived material in chunks of cells, i.e. the arrays of the input list are sliced along their last axis. The intermediate arrays of a chunk fit into the CPU caches. The number of evaluation points (quadrature points times cells) per chunk is given by the attribute ConstitutiveMaterial.chunksize (default is 8192, None evaluates all cells at once) and the chunks are evaluated by a pool of threads if the attribute ConstitutiveMaterial.parallel is True (default is False). Output arrays (out) are handled by the base class and results which are broadcasted along the cells are taken from the first chunk. Materials which implement the methods gradient() and hessian() are unchanged.
• Add the how-to guide Evaluate in parallel, which explains the parallel evaluation of the assembly (Job.evaluate(parallel=True)) and of the constitutive materials (e.g. NeoHooke(..., parallel=True)), both in chunks of cells, and the thread-safety of user-defined materials with _gradient() and _hessian().

Changed

• The state variables of ContactRigidPlane (the gap vectors, the reference gap vectors of the frictional contact and the flags of active contact points) are now only updated after a converged increment, like the state variables of a SolidBody. They are stored in results.statevars (converged) and results._statevars (trial). The forces and the stiffness matrix are functions of the displacements within an increment, evaluated by a return mapping which starts from the converged state. Before, the reference gap vectors were moved in each evaluation of the vector, e.g. also by trial evaluations of a LinesearchPlugin and by failed Newton-iterations, and the stiffness matrix depended on the order of the evaluations of the vector and the matrix. The attributes results.dx_ref, results.active and results.slip now hold the trial results of the last evaluation. The checkpoint of a ContactRigidPlane contains the state variables and the multipliers.
• The stiffness matrix of ContactRigidPlane is now consistent for frictional contact. For sliding points, the derivative of the direction of the friction force is evaluated with the trial (elastic) slip of the converged reference, and the non-symmetric coupling of the friction force to the normal force is included. For points which come into contact within an increment, the reference gap vector is located at the linearly interpolated point of first contact (zero gap), with its derivative in the stiffness matrix. Before, the reference was the gap vector of the first Newton-iteration with an active contact.
• The stiffness-based multipliers of ContactRigidPlane are evaluated by the mean stiffness of the diagonal blocks of the stiffness matrices of the items at all points in potential contact (before: the active points). Hence, they are defined for the first increment of a contact and they do not depend on the number of active points. The mean tangential stiffness is divided by the number of tangential directions (one in 2d, before: two). ContactRigidPlane.contact_multipliers() takes no arguments.
• Assemble the vector and the matrix of ContactRigidPlane from vectorized COO-arrays, instead of setting items of LIL-matrices. Before, the blocks of the points in contact were assembled from dense arrays of shape (n * dim, n * dim), i.e. the time and the memory scaled quadratically with the number of points in contact n (e.g. 31 s for 5000 points). Now, the assembly scales linearly (e.g. 14 ms for 5000 points).
• newtonraphson() creates one IterationState which is updated in-place during all iterations. Before, a new state was created after each iteration and the hooks before_newton and before_iteration got an empty state or the state of the previous iteration. The before_newton hook is now triggered after the unknowns are determined and before the initial assembly of the objective function.
• The check for NaN values in the Newton increment is now performed before the after_linear_solve hook is triggered.
• Copy StateNearlyIncompressible.u from field[0].values, instead of using it directly.
• Step.generate() catches errors of the Newton-Raphson method of a substep and triggers the after_failed_substep hook of the plugins. If the substep is not recovered by a plugin, the error is raised (unchanged).
• The unknowns x0 are now linked to the result of each completed substep by Step.generate() (before by Job.evaluate()). Hence, also a standalone Step.generate() starts each substep from the result of the previous substep.
• The hook after_substep gets the SubstepState of the substep (a subclass of JobState with the same attributes stepnumber, substepnumber and time). The ProgressPlugin shows the increments of subdivided substeps in this hook.
• The hook after_job is now also triggered if the evaluation of a job fails. The error is available as JobState.error and it is raised after all plugins are called.
• FieldContainer.checkpoint() shares the regions, meshes, mesh containers and indices of the fields with the field container, only the field containers, the fields and their values are copied. Before, the whole field container was deep-copied, including the regions with the gradients of the shape functions, e.g. 84 MB instead of 0.7 MB for a checkpoint of 8000 hexahedrons. This also applies to SolidBody.checkpoint(), SolidBodyNearlyIncompressible.checkpoint() and the CutbackPlugin, which creates a checkpoint before each substep. Changes of the region or the mesh after a checkpoint is created are visible in the field container of the checkpoint, the restore of the values is not affected. The field container of the checkpoint holds the copied top-level field container "field.x0" as x0 (before, x0 was copied twice).
• Change the perturbation in eigh and eigvalsh with new helpers in the JAX backend. The helper is enabled in becker(), storakers(), extended_tube(), ogden() and morph().
• Change the rate tensor in morph() to a symmetric tensor.
• Evaluate the bilinear form of the gradients of the test and the trial field in IntegralFormCartesian.integrate(parallel=False) in chunks of cells with a fixed contraction path. NumPy's default path limits the intermediate arrays to the size of the largest input array and hence, contracts the last three operands at once for many element types. The chunks are sized to fit into the CPU caches (assembly._cartesian.CHUNKSIZE_BYTES, 2 MiB, about the L2 cache of one core). This speeds up the integration of the stiffness matrix by about 1.2-2.5x (linear tetrahedrons), 3-5x (linear hexahedrons) and up to 25x for a single Lagrange hexahedron of order 4. The results are equal up to rounding errors.
• Assemble bilinear forms in IntegralFormCartesian.assemble() with a cached sparsity pattern of the test and the trial field. The values are summed up by numpy.bincount() into the data array of the CSR matrix, instead of converting all (duplicate) values from the COO to the CSR format on each assembly. The pattern is created on the first assembly (from the unique pairs of points per cell) and it is weakly referenced by the indices of the fields, i.e. it is not copied with the fields and it is released together with them. Fields with degrees of freedom which are not numbered point-wise fall back to the COO-assembly. This speeds up the assembly of the stiffness matrix by about 4-5x.
• Evaluate the COO-indices IntegralFormCartesian.indices of a bilinear form on demand. These are only required if the form is not assembled with a sparsity pattern.
• Return the assembled sparse matrix (vector) of a single block in IntegralForm.assemble(block=True) directly, without the block-assembly by scipy.sparse.bmat() (scipy.sparse.vstack()), which copied the matrix of a single field.
• Evaluate the linear form of the gradient of the test field in IntegralFormCartesian.integrate() by a single contraction (optimize=False). This is about 1.1-11x faster than NumPy's default contraction path, e.g. for the internal force vector.
• The update of a field container in newtonraphson() shares the regions, meshes and indices of the fields with the field container of the last iteration (like FieldContainer.checkpoint()), only the field containers, the fields and their values are copied. Before, the whole field container was deep-copied in each iteration, including the regions with the gradients of the shape functions. The arithmetic operators of a field container, e.g. x + dx, are not changed.
• Evaluate the chunks of the bilinear form of the gradients of the test and the trial field in IntegralFormCartesian.integrate(parallel=True) by a pool of threads (one per usable CPU of the process, e.g. the CPUs of a job on a shared HPC-node), instead of einsumt with NumPy's default contraction path. With a multi-threaded BLAS-library, limit its number of threads to avoid an oversubscription of the CPUs.
• Use the dtype from the field values for the statevars in SolidBody and SolidBodyNearlyIncompressible and ContactRigidPlane.
• Use the dtype from the field values for SparsityPattern.assemble().
• Keep the dtype in solve() and newtonraphson() when items are used.
• Assemble the vectors and matrices of MultiPointConstraint and MultiPointContact from vectorized COO-arrays, instead of setting items of LIL-matrices. Before, the sparse identity matrix of the connected points was converted to a dense array of shape (n, n) in the assembly of the matrix, i.e. the time and the memory scaled quadratically with the number of points n (e.g. 5.5 s for 10k points, 80 GiB for 100k points). Now, the assembly scales linearly, e.g. 3 ms for 10k points and 26 ms for 100k points. The results are unchanged.
• Move the evaluation of numpy.einsum() in chunks from assembly._cartesian to math._einsum.einsum_chunks(), together with math._einsum.cpu_count() (both are not part of the public API). The operands are now sliced by the chunked index of their subscripts (before: by the length of their last axis), i.e. operands without this index are not sliced, also if their last axis has the length of the cells. Optionally, the chunks are taken along another free index of the output (idx, default is the last axis of the output), the output array is created if it is not given (out=None) and the items are divided into one chunk per thread by default (chunksize=None).
• Evaluate the contractions of IntegralFormCartesian.integrate(parallel=True) (the linear forms and the bilinear forms without the gradients of both fields) and of assembly.expression.BasisField(parallel=True) by math._einsum.einsum_chunks() in chunks of cells (one per thread), instead of einsumt. The chunks are written directly into the output array, without concatenating them. This speeds up the parallel evaluation of the bases by about 2x and avoids the overhead of the threads if there are less cells than threads, e.g. for a single cell of high order. The serial evaluation is unchanged. The assembly does not use einsumt anymore.
• math._einsum.einsum_chunks() passes a single chunk in serial (parallel=False, chunksize=None) to numpy.einsum() without further checks, i.e. its result has the memory layout chosen by NumPy. The output array of more chunks is created in C-order, which is faster for the sparse assembly. The output array of the bilinear form of the gradients of both fields in IntegralFormCartesian.integrate() is always created in C-order, also for a single chunk.
• Evaluate the contractions of the math helpers math.dot(), math.ddot(), math.dddot(), math.cdya_ik(), math.cdya_il() and math.cdya() and of constitution.AreaChange with parallel=True by math._einsum.einsum_chunks() in chunks of cells (one per thread), instead of einsumt. The serial evaluation (parallel=False) is unchanged. FElupe does not use einsumt anymore.
• math._einsum.einsum_chunks() re-uses a pool of threads (one per usable CPU), which is created on its first use. Before, a new pool was created for each evaluation. A forked child process creates its own pool. This reduces the overhead of the threads for small arrays, e.g. the parallel assembly of the vector and the matrix of a mixed u/p/J-formulation with 64 hexahedrons is about 1.3-2x faster. Optional keyword-arguments are passed to numpy.einsum() and contractions without a free last index of the output, e.g. with a scalar output, are evaluated as a single chunk (instead of raising an error).
• The methods of a material decorated by constitutive_material() take precedence over the methods of ConstitutiveMaterial, i.e. the order of the base classes of the derived class is reversed. Hence, a class which is already derived from ConstitutiveMaterial may be decorated (before, this raised a TypeError).
• Move the pool of threads and cpu_count() from math._einsum to math._threads. Add math._threads.map_threaded(), which evaluates calls from a thread of the pool in serial, e.g. a math-function with parallel=True within a chunk of a constitutive material. Otherwise, a task of the pool which waits for further tasks of the same pool may block all threads of the pool.
• Evaluate the gradient and the hessian of NeoHooke and NeoHookeCompressible in chunks of cells, i.e. the materials implement _gradient() and _hessian() of ConstitutiveMaterial. This speeds up the evaluation of the hessian by about 1.4x and the assembly of the stiffness matrix of a SolidBody by about 1.2-1.4x (27k hexahedrons). With parallel=True, the chunks are evaluated by a pool of threads (before: the math-functions within the material), which speeds up the evaluation of the hessian by about 3x (2 CPUs). The output arrays are handled by ConstitutiveMaterial, i.e. the materials don't re-use output arrays on their own. The strain energy density functions are evaluated without threads.
• Evaluate the gradients and the hessians of the composed materials ThreeFieldVariation, NearlyIncompressible, ThirdMediumContactMixed, OgdenRoxburgh, CompositeMaterial and MaterialStrain in chunks of cells, i.e. they implement _gradient() and _hessian() of ConstitutiveMaterial and the (inner) material is evaluated for each chunk. This speeds up the evaluation of the hessians by about 1.2-2x (27k hexahedrons), e.g. ThreeFieldVariation(NeoHooke()) by about 2x. The composed materials don't pass output arrays to the (inner) material anymore.
• The argument parallel of ThreeFieldVariation, NearlyIncompressible and ThirdMediumContactMixed defaults to None, i.e. the flag is taken from the material, e.g. ThreeFieldVariation(NeoHooke(mu=1, bulk=5000, parallel=True)) evaluates the chunks of cells by a pool of threads. Add the argument parallel=None to OgdenRoxburgh. A CompositeMaterial evaluates the chunks by threads if any of its materials is parallel.
• The intermediate results of ThreeFieldVariation are local variables of a chunk (before: attributes of the material) and MaterialStrain does not change its keyword-argument tangent during the evaluation, i.e. both materials are thread-safe. Remove the unused fourth-order projection tensor in the hessian of ThreeFieldVariation.
• Evaluate the gradients and the hessians of the tensortrax-based materials Hyperelastic and MaterialAD in chunks of cells by ConstitutiveMaterial, instead of the chunks of tensortrax. This includes the conversions of the stresses and the elasticity tensors, e.g. of the second Piola-Kirchhoff stress tensor and its tangent, which speeds up the evaluation of the hessian of Hyperelastic(neo_hooke) by about 2x (27k hexahedrons). With parallel=True, the chunks are evaluated by the pool of threads of FElupe; tensortrax is evaluated without threads. The JAX-based materials are unchanged.
• With parallel=True, ConstitutiveMaterial evaluates a probe of two cells to create the output arrays and all chunks by the pool of threads. Before, the first chunk was evaluated in serial, which limited the speedup of the threads (e.g. by about 12 % for eight chunks and two threads).
• Evaluate the stress tensor of LinearElasticTensorNotation in chunks of cells, i.e. parallel=True evaluates the chunks by a pool of threads, like for all other materials. The elasticity tensor is constant and it is evaluated without threads (before: parallel=True evaluated the math-functions of the constant elasticity tensor by threads and the stress tensor without threads).
• Run the examples of the docstrings as doctests. The configuration of pytest is moved from tox.ini to pyproject.toml ([tool.pytest.ini_options]), because pytest ignored the section [tool:pytest] of tox.ini, i.e. the option --doctest-modules was never applied. The tests are collected from tests and src/felupe with --import-mode=importlib (the doctests of an installed package are evaluated for the installed package) and the doctest-flags NORMALIZE_WHITESPACE and ELLIPSIS. A conftest.py renders the plots off-screen and evaluates each doctest in a temporary directory, i.e. the files written by the examples don't end up in the working directory.
• Replace black, isort and flake8 by ruff, the linter and formatter configured in pyproject.toml. The rules of pycodestyle and pyflakes (E4, E7, E9, F), mutable default arguments (B006), sorted imports (I) and whitespace (W) are checked by pre-commit and for each pull request in the workflow lint.yml. The code is formatted with ruff, unused imports and trailing whitespace are removed and the default arguments of Circle(centerpoint, sections), Mesh.mirror(normal, centerpoint), mesh.mirror(normal, centerpoint), Mesh.add_runouts(values, centerpoint), mesh.runouts(values, centerpoint) and constitution.jax.vmap2(in_axes, out_axes, methods) are tuples instead of lists (the values are unchanged).

Fixed

• Fix the stiffness-based multipliers of ContactRigidPlane, which were evaluated at the degrees of freedom of the wrong points. The indices of the active points within points were used as point ids of the mesh.
• Fix the coupling blocks of the center-point in the stiffness matrix of ContactRigidPlane, which were transposed. This only affected non-symmetric blocks of sliding points.
• Fix the XDMFWriterPlugin and the AnimationWriterPlugin, which closed their file (or plotter) in the hook after_iteration if the Newton-Raphson method reached the maximum number of iterations or had NaN-valued norms. A substep which was recovered afterwards (e.g. by a cutback) could not be written. Both plugins now close their file only in after_job, which is also triggered if the evaluation fails. Before, the file was not closed for other errors (e.g. NaN values in the solution).
• Fix SolidBodyNearlyIncompressible.restore(): the internal fields (pressure and volume ratio) of the checkpoint are now restored exactly. Before, they were overwritten by the re-evaluation of the results, which updated the internal fields with the deformation gradient of the last evaluation (e.g. of a failed Newton iteration). With restore_state=False, the internal fields are re-initialized by the deformation of the restored field.
• Fix SolidBody.restore() for state variables which are initialized after the checkpoint was created, e.g. the temperature of a thermal.SolidBodyThermal. A copy of the state variables of the checkpoint is now assigned (was written in-place before, which failed for different shapes).
• Fix the assembly of bilinear forms with plane-strain fields and additional vector-valued fields in IntegralFormCartesian. Only the axes of two-dimensional fields are trimmed, e.g. a coupling block of shape (3, 3, 9) is now trimmed to (2, 2, 9) instead of (2, 2, 2).
• Fix the assembly of mixed-field formulations in IntegralFormAxisymmetric: the gradient flags of non-axisymmetric fields are now passed to the linear and bilinear forms, empty (None) blocks of bilinear forms are supported and an axisymmetric field may be coupled with a vector-valued field.
• Fix re-used output arrays in NeoHooke.gradient(x, out=P) and NeoHooke.hessian(x, out=A4): the arrays are now reset with fill(0.0). Before, P was not reset at all (wrong, accumulated stresses for NeoHooke(mu=None, bulk=...)) and A4 was multiplied by zero, which kept NaN and inf values in the output array.
• Fix CharacteristicCurve.res by using a property. This was always None before, because it was not linked to / updated by CharacteristicCurve._curve.res. The attributes x and y are not affected, because these initially empty lists are updated in-place.
• Fix repeated evaluations of a Job or a CharacteristicCurve. Before, Job.evaluate() added its built-in plugins, i.e. the ProgressPlugin and the XDMFWriterPlugin, to the dispatcher of the job on each call. Hence, the built-in plugins of all previous evaluations were triggered again, e.g. the progress was shown multiple times and the result file of a previous evaluation was overwritten. Now, Job.evaluate() creates a new dispatcher for each evaluation with the plugins of the job and the built-in plugins of this evaluation. The dispatcher of the job, Job.dispatcher, is not modified.
• Fix a wrong contraction in saint_venant_kirchhoff_orthotropic() (tensortrax backend).
• Fix CharacteristicCurve, when additional plugins are passed.
• Fix the matrix in ContactRigidPlane for center points, which are only partly constrained.
• Fix Boundary.apply_mask() and Boundary.plot() for boundary conditions with a dof-based mask (with Boundary.skip=None). Before, a point-based mask applied afterwards, e.g. by Boundary.select_surface_points(), and the plot raised a TypeError. Now, all axes of the points of a point-based mask applied after a dof-based mask are prescribed and the prescribed directions are plotted by the (dof-based) mask.
• Fix the parameters in the title of ConstitutiveMaterial.plot(), which raised a TypeError for parameters with a value of None (e.g. tangent of MaterialStrain) and for 0-d arrays (e.g. E of LinearElastic1D). Parameters with a value of None are not shown, 0-d arrays are shown like scalars and strings are shown in quotes.
• Fix the examples of the docstrings which failed as doctests, e.g. continuation lines with >>> instead of ..., undefined variables (loadcase, mesh), a missing hess=True of the region for the hessian of a basis, a missing dim=3 of the displacement field in ViewField and ViewSolid, an outdated output of mesh.Triangle and a None-valued block of the hessian of ThreeFieldVariation in the example of FormItem. Outputs of NumPy scalars are printed, i.e. they don't depend on the version of NumPy.

Removed

• Remove the optional dependency einsumt and the extra felupe[parallel]. The threaded evaluations with parallel=True use math._einsum.einsum_chunks() and require no additional packages. An installation with the extra felupe[parallel] still works, but the installer warns that the extra is not provided.

What's Changed

• Add mesh extrusion along paths with rotation-minimizing frames by @adtzlr in #1118
• Add ThirdMediumContactMixed by @adtzlr in #1119
• Fix re-used output arrays in NeoHooke by @adtzlr in #1120
• Add LinesearchPlugin, a backtracking line search for newtonraphson() by @adtzlr in #1121
• Add CutbackPlugin, a cutback of the increment of failed substeps by @adtzlr in #1123
• Share the geometry in FieldContainer.checkpoint() by @adtzlr in #1125
• Fix CharacteristicCurve.res attribute by using a @property by @adtzlr in #1127
• Add NewtonConvergenceError to recover only convergence failures by cutbacks by @adtzlr in #1129
• Use a fresh dispatcher for each Job.evaluate() by @adtzlr in #1130
• Change the rate L_G in morph to a symmetric tensor by @adtzlr in #1133
• Fix CharacteristicCurve with additional plugins passed by @adtzlr in #1135
• Enhance the performance of IntegralForm.integrate() by @adtzlr in #1137
• Use a cached SparsityPattern for the assembly of bilinear forms by @adtzlr in #1138
• Use lazy-evaluation for indices by @adtzlr in #1139
• No bmat() for the assembly of single-block sparse matrices / vectors by @adtzlr in #1140
• Enhance performance of vector integration by @adtzlr in #1141
• The updated field container shares the geometry with the old one by @adtzlr in #1142
• Re-add parallel-keyword to IntegralForm.integrate() by @adtzlr in #1143
• Fix MORPH material implementations by @adtzlr in #1144
• Don't change the dtype in SolidBody.assemble methods by @adtzlr in #1145
• Don't change the dtype in SolidBodyNearlyIncompressible.assemble methods by @adtzlr in #1146
• Fix centerpoint contribution in ContactRigidPlane by @adtzlr in #1147
• Add Mesh.select_surface_points() by @adtzlr in #1148
• Add Mesh.select_edge_points() by @adtzlr in #1149
• Exclude surface patches in select_surface_points() by @adtzlr in #1150
• Exclude edge patches in select_edge_points() by @adtzlr in #1151
• Add examples and initial selections to select_surface_points() and select_edge_points() by @adtzlr in #1153
• Add nonlinear_subdivision to select_surface_points() and hide mesh edges by default by @adtzlr in #1154
• Add interactive selection methods for Boundary by @adtzlr in #1155
• Add CylindricalConstraint by @adtzlr in #1156
• Vectorize the assembly of MultiPointConstraint and MultiPointContact by @adtzlr in #1157
• Add CartesianConstraint by @adtzlr in #1158
• Add rotations to MultiPointConstraint by @adtzlr in #1159
• Fix the state, the linearization and the multipliers of ContactRigidPlane by @adtzlr in #1160
• Add ContactRigidCylinder and a pressure-dependent tangential multiplier for rigid contacts by @adtzlr in #1161
• Move einsum_chunks() to felupe.math by @adtzlr in #1162
• Use einsum_chunks() in the assembly instead of einsumt by @adtzlr in #1163
• Use einsum_chunks() in the math helpers instead of einsumt by @adtzlr in #1164
• Remove the optional dependency einsumt by @adtzlr in #1165
• Evaluate constitutive materials in chunks of cells by @adtzlr in #1166
• Evaluate NeoHooke and NeoHookeCompressible in chunks of cells by @adtzlr in #1167
• Evaluate composed materials in chunks of cells by @adtzlr in #1168
• Evaluate the tensortrax materials in chunks of cells by FElupe by @adtzlr in #1169
• Document the parallel evaluation and make parallel consistent by @adtzlr in #1170
• Run the docstring examples as doctests by @adtzlr in #1171
• Replace black, isort and flake8 by ruff by @adtzlr in #1172
• Release 11.2.0 by @adtzlr in #1173

Full Changelog: v11.1.3...v11.2.0