Inner product of the strain field with a fourth order tensor

Hello,
How can I perform the inner product of the strain field with a fourth order tensor?

from dolfin import *
import numpy

mesh = Mesh('mesh.xml')

# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)

xi2 = numpy.array([[1,0,0], [0,1,0], [0,0,1]])
ns = numpy.array([2**(-1/2),2**(-1/2),0])
ndn = numpy.zeros((3,3))
n4 = numpy.zeros((3,3,3,3))

def epsilon(u):
    return sym(grad(u))

for i in range(2):
    for j in range(2):
      ndn[i,j] = ns[i]*ns[j]

for i1 in range(0,2):
    for i2 in range(0,2):
        for i3 in range(0,2):
            for i4 in range(0,2):
                    n4[i1,i2,i3,i4] = ndn[i1,i3]*xi2[i2,i4] + ndn[i1,i4]*xi2[i2,i3] + ndn[i2,i4]*xi2[i1,i3] + ndn[i2,i3]*xi2[i1,i4]

prod = inner(n4,epsilon(u))

An inner product has to be between something with the same shape. A fourth order tensor cannot be in an inner product with a second order tensor. I would suggest using dot. See Form language — Unified Form Language (UFL) 2021.1.0 documentation

1 Like

When using dot, it gives the following error:

Traceback (most recent call last):
File “product.py”, line 30, in
prod = dot(n4,epsilon(u))
File “/usr/lib/python3/dist-packages/ufl/operators.py”, line 172, in dot
a = as_ufl(a)
File “/usr/lib/python3/dist-packages/ufl/constantvalue.py”, line 437, in as_ufl
raise UFLValueError(“Invalid type conversion: %s can not be converted”
ufl.log.UFLValueError: Invalid type conversion:

How can I convert the array to the form that UFL can understand?

You need to use as_tensor(n4)

1 Like

Thank you for the answer. I have another question.
When I use a vector expression statement like

n = Expression(("-(3*a3*x[0]*x[0]+2*a2*x[0]+a1)/sqrt(pow((3*a3*x[0]*x[0]+2*a2*x[0]+a1),2)+1)","1.0/sqrt(pow((3*a3*x[0]*x[0]+2*a2*x[0]+a1),2)+1)"),a3=p0[0],a2=p0[1],a1=p0[2],degree=0)

Can I use it in tensor algebra operations? Also, how can I display the values of the expression statements and tensors to see if it works correctly or not?

I am not sure what you mean by using tensor algebra operations in this setting.
The code you are writing has to be “C++” code, so that one can generate efficient code for it.

To visualize an expression you should interpolate it into a suitable function space (in your case DG, 0, as you are using degree=0) and visualize it in for instance paraview

1 Like

I mean such as outer(n,n) or inner(n,n). Also, how can I display the elements of the tensor such as as_tensor(n4) in the previous code?

No, that is not supported in legacy dolfin.

Again, project the as_tensor into a suitable function space.

1 Like

Can’t I use functions that are solutions to the variational problem like displacement in outer or inner operations? So, how can I create a function over the mesh that its values are assigned externally to use in outer or inner operations?

You can use them in inner and outer operations, But an Expression does not support the inner and outer operations (at least as far as I am aware).

inner, outer, etc can thus be used to assemble scalar values, vectors or matrices.

Without an example illustrating what you would like to do with the inner/outer products, i.e. What do you want to do with

or the equivalent corrected expression

1 Like

Let’s say I have an array:

coor = V.tabulate_dof_coordinates()
length = len(coor)
n = numpy.zeros((length,2))
for i in range(length):
    n[i,0] = coor[i,0]
    n[i,1] = coor[i,1]

How can I make this array a function for a vector function space over the mesh and use it in inner or outer operations?

See for instance Application of point forces, mapping vertex indices to corresponding DOFs - #2 by dokken

1 Like