diff --git a/CHANGELOG.md b/CHANGELOG.md
index fbed6cfb..552e95ae 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,9 @@ All notable changes to this project will be documented in this file. The format
## [Unreleased]
+### Added
+- Add `math.solve_2d(A, b, solver=np.linalg.solve, **kwargs)` to be used in `newtonrhapson(solve=solve_2d, ...)` for two-dimensional unknowns. This is useful for local Newton-iterations related to viscoelastic evolution equations inside constitutive material formulations.
+
### Fixed
- Sort array of principal values in descending order before plotting in `Scene.plot("Principal Values of ...")`. This ensures that the labels are matching user-defined arrays of principal values.
diff --git a/docs/felupe/math.rst b/docs/felupe/math.rst
index 8b65c1fb..25af760e 100644
--- a/docs/felupe/math.rst
+++ b/docs/felupe/math.rst
@@ -56,7 +56,13 @@ This module contains math functions.
math.reshape
math.ravel
+** Equation system**
+
+.. autosummary::
+
+ math.solve_2d
+
**Detailed API Reference**
.. automodule:: felupe.math
- :members: linsteps, rotation_matrix, deformation_gradient, strain, strain_stretch_1d, extract, values, norm, interpolate, grad, identity, sym, dya, inv, det, dev,cof, eig, eigh, eigvals, eigvalsh, equivalent_von_mises, transpose, majortranspose, trace, cdya_ik, cdya_il, cdya, cross, dot, ddot, tovoigt, reshape, ravel
+ :members: linsteps, rotation_matrix, deformation_gradient, strain, strain_stretch_1d, extract, values, norm, interpolate, grad, identity, sym, dya, inv, det, dev,cof, eig, eigh, eigvals, eigvalsh, equivalent_von_mises, transpose, majortranspose, trace, cdya_ik, cdya_il, cdya, cross, dot, ddot, tovoigt, reshape, ravel, solve_2d
diff --git a/src/felupe/math/__init__.py b/src/felupe/math/__init__.py
index e74a13d8..9acf112a 100644
--- a/src/felupe/math/__init__.py
+++ b/src/felupe/math/__init__.py
@@ -12,6 +12,7 @@
values,
)
from ._math import linsteps
+from ._solve import solve_2d
from ._spatial import rotation_matrix
from ._tensor import (
cdya,
@@ -77,4 +78,5 @@
"tovoigt",
"trace",
"transpose",
+ "solve_2d",
]
diff --git a/src/felupe/math/_solve.py b/src/felupe/math/_solve.py
new file mode 100644
index 00000000..9b4a50cf
--- /dev/null
+++ b/src/felupe/math/_solve.py
@@ -0,0 +1,85 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe. If not, see .
+"""
+
+import numpy as np
+
+
+def solve_2d(A, b, solve=np.linalg.solve, **kwargs):
+ r"""Solve a linear equation system for two-dimensional unknowns with optional
+ elementwise-operating trailing axes.
+
+ Parameters
+ ----------
+ A : ndarray of shape (M, N, M, N, ...)
+ The left-hand side array of the equation system.
+ b : ndarray of shape (M, N, ...)
+ The right-hand side array of the equation system.
+ solve : callable, optional
+ A function with a compatible signature to :func:`numpy.linalg.solve`
+ (default is :func:`numpy.linalg.solve`).
+ **kwargs : dict, optional
+ Optional keyword arguments are passed to the callable ``solve(A, b, **kwargs)``.
+
+ Returns
+ -------
+ x : ndarray of shape (M, N, ...)
+ The solved array of unknowns.
+
+ Notes
+ -----
+ The first two axes of the rhs ``b`` and the first four axes of the lhs ``A`` are the
+ tensor dimensions and all remaining trailing axes are treated as batch dimensions.
+ This function finds :math:`x_{kl}` for Eq. :eq:`solve-system`.
+
+ .. math::
+ :label: solve-system
+
+ A_{ijkl} : x_{kl} = b_{ij}
+
+ Examples
+ --------
+ >>> import numpy as np
+ >>> import felupe as fem
+ >>>
+ >>> np.random.seed(855436)
+ >>>
+ >>> A = np.random.rand(3, 3, 3, 3, 2, 4)
+ >>> b = np.ones((3, 3, 2, 4))
+ >>>
+ >>> x = fem.math.solve_2d(A, b)
+
+ >>> x[..., 0, 0]
+ array([[ 1.1442917 , 2.14516919, 2.00237954],
+ [-0.69463749, -2.46685827, -7.21630899],
+ [ 4.44825615, 1.35899745, 1.08645703]])
+
+ >>> x.shape
+ (3, 3, 2, 4)
+
+ """
+
+ size = np.prod(b.shape[:2])
+ trax = b.shape[2:]
+
+ # flatten and reshape A to a 2d-matrix of shape (..., M * N, M * N) and
+ # b to a 1d-vector of shape (..., M * N)
+ b_1d = np.einsum("i...->...i", b.reshape(size, np.prod(trax)))
+ A_1d = np.einsum("ij...->...ij", A.reshape(size, size, np.prod(trax)))
+
+ # move the batch-dimensions to the back and reshape x
+ return np.einsum("i...->...i", solve(A_1d, b_1d, **kwargs)).reshape(b.shape)
diff --git a/tests/test_constitution_newton.py b/tests/test_constitution_newton.py
new file mode 100644
index 00000000..01304262
--- /dev/null
+++ b/tests/test_constitution_newton.py
@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+"""
+ _______ _______ ___ __ __ _______ _______
+| || || | | | | || || |
+| ___|| ___|| | | | | || _ || ___|
+| |___ | |___ | | | |_| || |_| || |___
+| ___|| ___|| |___ | || ___|| ___|
+| | | |___ | || || | | |___
+|___| |_______||_______||_______||___| |_______|
+
+This file is part of felupe.
+
+Felupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+Felupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Felupe. If not, see .
+
+"""
+import numpy as np
+import tensortrax as tr
+from tensortrax.math import trace
+from tensortrax.math.linalg import inv
+from tensortrax.math.special import dev, from_triu_1d, triu_1d
+
+import felupe as fem
+
+
+def test_visco_newton():
+ mesh = fem.Rectangle(n=3)
+ region = fem.RegionQuad(mesh)
+ field = fem.FieldContainer([fem.FieldPlaneStrain(region, dim=2)])
+ boundaries, loadcase = fem.dof.uniaxial(field, clamped=True)
+
+ @fem.constitution.isochoric_volumetric_split
+ def finite_strain_viscoelastic_newton(C, Cin, mu, eta, dtime):
+ "Finite Strain Viscoelastic material formulation."
+
+ def evolution(Ci, Cin, C, mu, eta, dtime):
+ "Viscoelastic evolution equation."
+ return mu / eta * dev(C @ inv(Ci)) @ Ci - (Ci - Cin) / dtime
+
+ # update of state variables by evolution equation
+ Cin = from_triu_1d(Cin, like=C)
+ Ci = fem.newtonrhapson(
+ x0=Cin,
+ fun=tr.function(evolution, ntrax=C.ntrax),
+ jac=tr.jacobian(evolution, ntrax=C.ntrax),
+ solve=fem.math.solve_2d,
+ args=(Cin, C.x, mu, eta, dtime),
+ verbose=0,
+ tol=1e-2,
+ ).x
+
+ # first invariant of elastic part of right Cauchy-Green deformation tensor
+ I1 = trace(C @ inv(Ci))
+
+ # strain energy function and state variable
+ return mu / 2 * (I1 - 3), triu_1d(Ci)
+
+ umat = fem.Hyperelastic(
+ finite_strain_viscoelastic_newton, mu=1.0, eta=1.0, dtime=0.1, nstatevars=6
+ )
+ solid = fem.SolidBodyNearlyIncompressible(umat, field, bulk=5000)
+ solid.results.statevars[[0, 3, 5]] += 1
+
+ move = fem.math.linsteps([0, 0.3], num=3)
+ step = fem.Step(
+ items=[solid], ramp={boundaries["move"]: move}, boundaries=boundaries
+ )
+ job = fem.CharacteristicCurve(steps=[step], boundary=boundaries["move"]).evaluate()
+ assert np.all([norm[-1] < 1e-8 for norm in job.fnorms])
+
+
+if __name__ == "__main__":
+ test_visco_newton()