How to Perform Inner Product for Tensor Coefficients?

[Dreams-Cindy]

Dear Developers,

I am currently using Shenfun to solve the following equation on a disk. Using polar coordinates, the degrees of freedom in the radial and angular directions are set to N. The dimensions of both u and f are $2 \times N \times N$.

The equation is given by:

where C is a known fourth-order tensor with dimensions $2 \times 2 \times 2 \times 2$.

The way the fourth-order tensor $C_{ijkl}$ acts on the second-order tensor $\nabla u$ is:

I would like to know how to perform the inner product for a known tensor coefficient C.

It is clear that directly multiplying by C is incorrect:

inner(v,−div(C*∇u)).

I would appreciate guidance on how to resolve this problem. Thank you!

[mikaem]

Hi

Using a fourth order tensor like this has never been implemented before in Shenfun as far as I know. It is possible, but in order to do it you have to understand how shenfun expressions are stored under the hood. Basically, you need to create an Expr, using your tensor coefficient C to modify the scales of all terms. A naive implementation for the polar coordinates would be like shown below. Please test it, because I'm not 100% sure it is correct:-)

from shenfun import *
import sympy as sp 
import numpy as np

# Define polar coordinates using angle along first axis and radius second
theta, r = psi = sp.symbols('x,y', real=True, positive=True)
rv = (r*sp.cos(theta), r*sp.sin(theta))

N = 32
F = FunctionSpace(N, 'F', dtype='d')
L = FunctionSpace(N, 'L', bc=(0, 0), domain=(0, 1))
T = TensorProductSpace(comm, (F, L), axes=(1, 0), coordinates=(psi, rv))
V = VectorSpace(T)
u = TrialFunction(V)
du = grad(u)
C = np.arange(16).reshape((2, 2, 2, 2))

terms, scales, indices = [], [], []
for i in range(2):
     for j in range(2):
         f = []
         for k in range(2):
             for l in range(2):
                 f.append(C[i, j, k, l] * du[k][l])
         g = sum(f[1:], f[0])
         terms += g.terms()
         scales += g.scales()
         indices += g.indices()
Cgradu = Expr(du._basis, terms, scales, indices)

[Dreams-Cindy]

Thank you very much for your prompt reply. I will try to test it.

[Dreams-Cindy]

Thanks to the code you provided, I am now able to implement the inner product for tensor coefficients. I greatly appreciate your prompt assistance!