Transposed is only defined for rank 2 tensors

Linear Algebra
1

Hey,

i have been trying to implement a PDE for Fiber Orientation. The output should be a 3x3 Tensor so i defined a TensorFunctionSpace with shape 3x3 both for Orientation o and velocity u.
Nevertheless, when solving the nabla_grad(u).T FEniCs recognizes o as a higher rank tensor (above 2).
Why is that? Did i accidentally specified it somewhere?
Here is my code:

from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np

T = 5.0            # final time
num_steps = 100    # number of time steps
dt = T / num_steps # time step size
eps = 0.01         # diffusion coefficient

# Read mesh from file
mesh = UnitSquareMesh(8, 8, "left/right") #creating 3x3 tensor on 2d mesh

# Define function space for velocity
W = TensorFunctionSpace(mesh, 'P', 1, shape = (3,3))

V = TensorFunctionSpace(mesh, "Lagrange", 2, shape = (3,3))

#Define boundaries

inflow  = 'near(x[0], 0)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

#Define boundary conditions

u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2) 
o_D = np.zeros([3,3]) 
bco_inflow = DirichletBC(W, o_D, inflow)
bco_noslip = DirichletBC(W, o_D, walls)  Constant((0,0,0))
bco = [bco_noslip]

#Define Test-function

v = TestFunction(W)

# Define functions for velocity and orientations

u = TrialFunction(W)    #velocity
o = TrialFunction(V)    #orientation
o_n = Function(V)  #orientatin at next time step

#Define source terms

def symepsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) #symmetric part of gradient of velocity

def asymepsilon(u):
    return 0.5*(nabla_grad(u) - nabla_grad(u).T) #antisymmetric part of gradient of velocity

def fourth(o):
    return inner(o, o)  

#Define expressions used in variational form

Lamba = Constant(0.5) #assume Lenght to Thickness ratio is between 0 an 1
#eps = Constant(eps) #for symmetric and antisymmetric part of velocity
D = Constant(0.003) #tendency of fibers to adopt an isotropic orientation state
delta = Identity(dim = 2)

#Define variational Problem

F = ((-0.5*v*(inner(asymepsilon(u), skew(o))- inner(skew(o), asymepsilon(u)))) + \
            (0.5*v*Lambda*(inner(symepsilon(u), sym(o)) + inner(sym(o), symepsilon(u)) - 2*inner(symepsilon(u), fourth(o)))) + \
            (2*D*v*(delta - 3*(o))))

I get the error message:

Transposed is only defined for rank 2 tensors.

---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-6-ee03bdbfc5bd> in <module>
     69 
     70 F = ((-0.5*v*(inner(asymepsilon(u), skew(o))- inner(skew(o), asymepsilon(u)))) + \
---> 71             (0.5*v*Lambda*(inner(symepsilon(u), sym(o)) + inner(sym(o), symepsilon(u)) - 2*inner(symepsilon(u), fourth(o)))) + \
     72             (2*D*v*(delta - 3*(o))))
     73 #Weak form first without impulse bilance; stress tensor

<ipython-input-6-ee03bdbfc5bd> in asymepsilon(u)
     53 
     54 def asymepsilon(u):
---> 55     return 0.5*(nabla_grad(u) - nabla_grad(u).T) #antisymmetric part of gradient of velocity
     56 
     57 def fourth(o):

~/anaconda3/envs/fenicsproject/lib/python3.6/site-packages/ufl/exproperators.py in _transpose(self)
    350     """Transpose a rank-2 tensor expression. For more general transpose
    351     operations of higher order tensor expressions, use indexing and Tensor."""
--> 352     return Transposed(self)

I know this will probably be a pretty simple question but i am new to FEniCs and struggling a bit with the math.

Thanks in advance

2

Hi
I believe both the fiber orientation and velocity should be vectors and not a (3,3) tensors. Consider changing the definition of W and V to start with:

W = VectorFunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh,  "CG", 2)

As you can see from the error message, the gradient of a second order tensor would be a third order tensor and the symmetric part of the gradient of a third order tensor is meaningless.

Also in the definition of symepsilon, you can simply use sym(grad(u)). Similarly for the skew symmetric part skew(grad(u)) I believe

1 Like
3

Thanks for your reply.
The error message disappeared but now it has a problem with my Dirichlet BC.
If I define the FunctionSpaces am I still able to define a Tensor as a Dirichlet Boundary condition? The final goal of the simulation is to get the exact value for the Orientation Tensor (3x3) at every time in the simulation.

The Orientation Tensor at x = 0 and t=0 is supposed to be the unit tensor and the velocity at t=0 freely selectable.

#Define boundaries

inflow  = 'near(x[0], 0)'
walls   = 'near(x[1], 0) || near(x[1], 1)'

#Define boundary conditions

u_D = np.array([[1,1,1], [0,0,0], [0,0,0]])
o_D = np.array([[1,0,0], [0,1,0], [0,0,1]]) 
bco_inflow = DirichletBC(W, o_D, inflow)
bco_noslip = DirichletBC(W, o_D, walls) 
bco = [bco_noslip]

I am sure there is something wrong with the BCs itself.
The error message I get is:

*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  Illegal value rank (2), expecting (1).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  15b823f2c45d3036b6af931284d0f8e3c77b6845

Thanks in advance.

4

Hi, W is a VectorFunctionSpace, so your u_D and o_D should also be 3d vector.