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

194

views

2

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:

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.

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

1

I'm unable to reproduce your problem, because the 'polar' function isn't define...:

Where are you getting that function? You might try to decompose it using other means, using

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.

`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)
#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.

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...

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.

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.

^{T}F with its spectal decomposition.https://www.allanswered.com/post/omlpn/

The algorithm in Simo and Hughes book seems easy to implement within FENICS and stable.

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.

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:

For strain energies that depend on the eigenvalues of strain only, I've also tried using the following:

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.

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