### How do I get the polar decomposition of a field?

194
views
2
3 months ago by
I'm solving a problem in nonlinear elasticity.  As a part of this, I need to take the polar decomposition of the deformation gradient field $F = grad(y)$.  When I do this naively, I get the following error:

R,U=polar(F)
File "/usr/lib/python3/dist-packages/scipy/linalg/_decomp_polar.py", line 102, in polar
raise ValueError("a must be a 2-D array.")

My MWE that reproduces this error is as follows:
from fenics import *
from scipy.linalg import *

mesh = RectangleMesh(Point(-1,-1),Point(1,1), 4, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

def W1(y):
R,U=polar(F)
return tr(U*U) # test strain energy density

y  = Function(V)
Pi=W1(y)*dx



I think I see the problem -- that $F$ is a tensor field and not a tensor? -- but I don't know how to fix that.
Community: FEniCS Project
1
Thanks Marco.  The links are very useful.  While looking, I also found this:
The algorithm in Simo and Hughes book seems easy to implement within FENICS and stable.
written 3 months ago by K D
Hello Marco and others,

I followed the procedure at the first link that you pasted above.  It works great to find eigenvalues, and then I use the algorithm mentioned above from Simo and Hughes to construct $U$.

One minor issue is that finding the eigenvalues involves a square root.  The argument is always non-negative theoretically.  However, when it is very small, round-off errors can cause it become negative and gives me a "nan" result.  A quick fix is to simply use an absolute operator inside the square root.  I wonder if anybody has a cleaner fix (this seems a bit hacky to me).

I'll paste my final code snippet here once I get it going for the record.

written 3 months ago by K D
2
I was looking into this myself some time ago, and compared a few different approaches (with no clear winner).  My code for the algorithm from Simo and Hughes was
# Algorithm from Simo and Hughes, "Computational Inelasticity", Box 7.1.
# - Not friendly toward automatic differentiation, which removes a major
# - An advantage of this over Denman--Beavers is that you get the eigenvalues
#   at an intermediate step, which may be useful, e.g., for computing the
#   logarithmic/Hencky strain.
# - If the tolerance is too small, the resulting displacement starts to
#   deviate very significantly from the other two options.
TOL3 = Constant(1e-5) #DOLFIN_EPS
def getU(C):
I1 = tr(C)
I2 = 0.5*(I1**2 - tr(C*C))
I3 = det(C)
b = I2 - I1*I1/3.0
c = -2.0/27.0*I1**3 + I1*I2/3.0 - I3
m = 2.0*sqrt(-b/3.0)
n = 2.0*c/(m*b)
# apparently there is no atan2(?)
t = atan(sqrt(1.0-n**2)/n)/3.0
x = []
lam = []
for A in range(0,3):
# must guard against c<0 to avoid nan
xA = conditional(lt(abs(b),TOL3),
-pow(Max(c,Constant(0.0)),1.0/3.0),
m*cos(t+2.0*(A-1.0)*pi/3.0))
x += [xA,]
lam += [sqrt(xA + I1/3.0),]
i1 = lam[0] + lam[1] + lam[2]
i2 = lam[0]*lam[1] + lam[0]*lam[2] + lam[1]*lam[2]
i3 = lam[0]*lam[1]*lam[2]
D = i1*i2 - i3
U = (1.0/D)*(-C*C + (i1**2 - i2)*C + i1*i3*Identity(3))

return U​
Something else that I found more robust and amenable to automatic differentiation (albeit of questionable accuracy) is to approximate the matrix square root of $\mathbf{C}$ using a fixed number of Denman--Beavers iterations:
# Denman--Beavers iteration, see, e.g.,
#
#  https://en.wikipedia.org/wiki/Square_root_of_a_matrix#By_Denman%E2%80%93Beavers_iteration
#
# - Does not pose problems for automatic differentiation, but UFL code may blow
#   up for large number of iterations
# - There is no convergence check on how accurate the square root is.
def matSqrt(n,Y,Z=Identity(3)):
if(n==0):
return Y
Y1 = 0.5*(Y + inv(Z))
Z1 = 0.5*(Z + inv(Y))
return matSqrt(n-1,Y1,Z=Z1)​
For strain energies that depend on the eigenvalues of strain only, I've also tried using the following:
# Algorithm from Wikipedia to get eigenvalues:
#
# https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
#
# plus some trial-and-error epsilons that I added. This is only
# useful if you have a strain-energy that can be computed from eigenvalues
# alone, as it does not compute the eigenvectors, and convergence is shaky
# at best in problems without a mass matrix.
def eig3x3(A):
A11 = A[0,0]
A12 = A[0,1]
A13 = A[0,2]
A21 = A[1,0]
A22 = A[1,1]
A23 = A[1,2]
A31 = A[2,0]
A32 = A[2,1]
A33 = A[2,2]
p1 = A12*A12 + A13*A13 + A23*A23
q = tr(A)/3.0
p2 = (A11-q)**2 + (A22-q)**2 + (A33-q)**2 + 2.0*p1
p = sqrt(p2/6.0)
B = (1.0/p)*(A - q*I)
r = 0.5*det(B)
phi = conditional(gt(-1.0+DOLFIN_EPS, r), pi/3.0, \
conditional(gt(r,1.0-DOLFIN_EPS),Constant(0.0),\
acos(r)/3.0))
eig1 = conditional(gt(p1,DOLFIN_EPS),\
q + 2.0*p*cos(phi),A11)
eig3 = conditional(gt(p1,DOLFIN_EPS),\
q + 2.0*p*cos(phi + 2.0*pi/3.0),A33)
eig2 = 3.0*q - eig1 - eig3

return eig1, eig2, eig3​

However, taking the second derivative of the strain energy, to get a consistent tangent, is unstable, and it only seems to work reliably with implicit methods in dynamic problems, where the mass matrix stabilizes the Newton iteration.

written 3 months ago by David Kamensky

Thanks David!  Your codes are very helpful.  For now, I'm in 2D, so it looks friendlier to automatic differentiation, though I haven't tried yet.

PS.  For reasons unknown to me, FEniCS implements atan2 as atan_2.  See here:
https://fenicsproject.org/olddocs/ufl/1.5.0/ufl.html

def W1(y):
C=F.T*F

L1 = (tr(C) + sqrt(abs(tr(C)**2-4.*det(C))))/2.    # bigger eigenvalue of C
L2 = (tr(C) - sqrt(abs(tr(C)**2-4.*det(C))))/2.    # smaller eigenvalue of C
L3 = 1.0 #unit  eigenvalue of C
l1,l2,l3=L1**0.5,L2**0.5,L3**0.5 # eigenvalues of U
i1, i2, i3 = l1 + l2 + l3, l1*l2+l2*l3+l3*l1, l1*l2*l3 # invariants of U
D = i1*i2-i3 # see Simo and Hughes Box 7.1
U = (-C*C + (i1**2-i2)*C + i1*i3*Identity(2))/D # see Simo and Hughes Box 7.1
written 3 months ago by K D

1
3 months ago by
I'm unable to reproduce your problem, because the 'polar' function isn't define...:

NameError: global name 'polar' is not defined

Where are you getting that function? You might try to decompose it using other means, using
R = F.T*F​
but I'm not sure if UFL has a Hermitian sqrt function either...

Incidentally, the MR strain energy density doesn't require you to do this. You can easily find the rotation-independent tensors (Green or Lagrangian) which are related to F.T * F (and don't require the square root or the rotation matrix explicitly). The det(F) is then related to det(C) or det(E) trivially.

from fenics import *

mesh = RectangleMesh(Point(-1,-1),Point(1,1), 4, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

def W1(y):
#R,U=polar(F)
C = F.T*F
detF = sqrt(det(C))
return tr(C) - ln(detF) + (detF-1)**2 # strain energy density (compressible Mooney-Rivlin)

y  = Function(V)
Pi = W1(y)*dx​
Hello Mike,

Thank you for trying this out.  I had missed an import statement that is now fixed in the MWE.  It should work now.

About the energy: I am just using Mooney-Rivlin as a placeholder to put together the MWE.  The actual energy that I want to use requires $U$ and is a complex expression.
written 3 months ago by K D
Ah, I suspected it might be from a different package. Per my understanding, scipy is expecting basically a numpy array (i.e.: the value of your tensor field at a single point). So you could technically dump the tensor field to a numpy array, process it in scipy / numpy and reinsert it into your U and P functions.
Or as Marco says there may be 'better' solutions available.

Personally, I'd still recommend trying to formulate all of this in terms of C or E, or the strain tensor invariant for which you also don't need to decompose F... unless the rotation of your material really affects the answer...
written 3 months ago by Mike Welland
Hello Mike,

My material model is anisotropic, and I have quantities of the form $(n.U.n)^2$ where $n$ is a given material direction, and I could not find a way to write it in terms of $C$ or $E$.

I think I will try some version of Marco's suggestion because it seems possible entirely within the FENICS framework.

For reference however, can you point to how I would "dump the tensor field to a numpy array, process it in scipy / numpy and reinsert it into your U and P functions"?  Thanks.
written 3 months ago by K D