How to do polar decomposition (RU decomposition) of the deformation gradient tensor in FEniCS?

General
1

In the below MWE, I have a unit cube. I have defined a displacement field U1 = Expression(("0.0","cos(theta) - sin(theta)","sin(theta) + cos(theta)"), theta = pi/3, degree=0). I projected it on a VectorFunctionSpace defined on the domain. Next, I define the deformation gradient tensor as F = Identity(3) + grad(disp). I want to do a polar decomposition as R,U = polar(F) where polar() is imported from scipy.linalg. This gives me the error

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

Is there a way to do the polar decomposition of the deformation gradient tensor pointwise on the entire domain?

Here is the MWE:

from dolfin import *  
from scipy.linalg import polar   

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}
               


# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
T = TensorFunctionSpace(mesh, "Lagrange", 1)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

vtk1 = File("R_ref.pvd")
vtk2 = File("U_ref.pvd")
              
U1 = Expression(("0.0",
                "cos(theta) - sin(theta)",
                "sin(theta) + cos(theta)"), theta = pi/3, degree=0)
                
disp = project(U1, V)

F = Identity(3) + grad(disp)

R,U = polar(F)

R_ref = project(R, T)
U_ref = project(U, T)

R_ref.rename('R_ref','R_ref')
T_ref.rename('T_ref','T_ref')

vtk1<<R_ref
vtk2<<U_ref

Any suggestion or comments would be extremely helpful.

Thanks in advance!

2

F is a symbolic expression. It is not an array of values.
You would have to implement the polar decomposition in UFL. ie. take the logic from: scipy/scipy/linalg/_decomp_polar.py at v1.15.2 · scipy/scipy · GitHub
and implement it with ufl operators, similar to:
Calculating local eigenvalues and eigenvectors for use in weak form
(There are analytical formulas for svd of 2x2 matrices at): Singular value decomposition - Wikipedia

1 Like
3

I am not sure I understand. So are you saying that we need to run a for loop over the vertices of a tetrahedral mesh, convert F (in the MWE I posted originally) into an array locally at each vertex, do the SVD, and then reconfigure it as a matrix.

I am sorry this is a trivial thing that I am not getting. Would it be possible to give me an MWE on how to do the SVD or RU decomposition of the deformation gradient tensor?