How do I get the polar decomposition of a field?

3 months ago by
K D  
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:

File "/usr/lib/python3/dist-packages/scipy/linalg/", 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):
    F = grad(y)
    return tr(U*U) # test strain energy density
y  = Function(V)

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
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  
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
#   advantage of using FEniCS.
# - 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),
        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.,
# - 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)):
        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:
# 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, \
    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:

def W1(y):
    F = grad(y)
    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 Answer

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):
    F = grad(y)
    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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »