Constructing 4th order tensor

1

Dear all,

I am a beginner to the FEniCS. I wish to construct the 4th order tensor in 2D as shown in the answer of the link here . I think ufl can be helpful for me. But I don’t know how to proceed. Help is appreciated.

from dolfin import *
from fenics import *
from ufl import *

mesh = UnitSquareMesh(2,2)
i,j,k,l = Indices(4)

After this I wish to construct the tensor as shown in the link pasted above.
Also, 2nd part of my question, if I have to modify the certain entry of the higher order tensor, how can I? Let’s say I have 3rd order constant tensor in 2D of name ‘K’. And I have to modify the entry corresponding to (1,2,2) and assign it a value, let’s say 2. I can’t perform K(1,2,2) = 2. It throws up an error. How should I perform this assignment?

Dolfin version: 2017.2.0
OS : Ubuntu

Thanks for the help.
Regards,
Peter

2

Is the following useful?

from dolfin import *
from fenics import *
from ufl import *

mesh = UnitSquareMesh(2,2)
i,j,k,l = indices(4)

# Tensor from StackExchange thread:
d = mesh.geometry().dim()
delta = Identity(d)
I = as_tensor(0.5*(delta[i,k]*delta[j,l] + delta[i,l]*delta[j,k]),(i,j,k,l))

# Kludgy way to change components:  Fill up a nested list with the components,
# change elements of the list, then use as_tensor to turn the list back
# into a tensor.
componentList = []
for i in range(0,d):
    componentList += [[],]
    for j in range(0,d):
        componentList[i] += [[],]
        for k in range(0,d):
            componentList[i][j] += [[],]
            for l in range(0,d):
                componentList[i][j][k] += [I[i,j,k,l],]

# Allowed to assign items in a list:
componentList[0][1][1][1] = 2.0

# Can convert lists into tensors:
I = as_tensor(componentList)
2 Likes
3

Thanks David Kamensky! This is what I wanted.
I pasted the following code below

print shape(I)
A = as_tensor([[1.0,2.0],[0.5,0.5]])
B = as_tensor(I[i,j,k,l]*A[k,l])
print (shape(B))

As far as the shape of I is concerned, as expected, it prints

(2,2,2,2)

But when I take the inner product with A, it outputs

()

I expect it to provide (2,2), because I am performing double contraction with tensor of order 2. Is there anything I am missing?

Thanks in advance.

4

The following modification gives the expected shape:

from dolfin import *
from ufl import *

mesh = UnitSquareMesh(2,2)
i,j,k,l = indices(4)

# Tensor from StackExchange thread:
d = mesh.geometry().dim()
delta = Identity(d)
I = as_tensor(0.5*(delta[i,k]*delta[j,l] + delta[i,l]*delta[j,k]),(i,j,k,l))

print(shape(I))
A = as_tensor([[1.0,2.0],[0.5,0.5]])

####### Added tuple of free indices as second argument to as_tensor() ######
B = as_tensor(I[i,j,k,l]*A[k,l],(i,j))

print(shape(B))
3 Likes
5

Thanks for the reply. But when I do it using lists and then use as_tensor to turn the list back into a tensor, it fails.

from dolfin import *
from ufl import *

mesh = UnitSquareMesh(2,2)
i, j, k,l = indices(4)

d = mesh.geometry().dim()
delta = Identity(d)
I = as_tensor(0.5*(delta[i,k]*delta[j,l] + delta[i,l]*delta[j,k]),(i,j,k,l))

componentList =

for i in range(0,d):
componentList += [,]

for j in range(0,d):
    componentList[i] += [[],]

    for k in range(0,d):
        componentList[i][j] += [[],]

        for l in range(0,d):
            componentList[i][j][k] += [I[i,j,k,l],]

componentList[0][1][1][1] = 2.0

I = as_tensor(componentList)
print shape(I)
A = as_tensor([[1.0,2.0],[0.5,0.5]])
B = as_tensor(I[i,j,k,l]*A[k,l],(i,j))
print (shape(B))

Sorry for the inappropriate indentation in the code. I could not find the proper formatting option. :sob:
Thanks again in advance.

6

This is because i, j, k, and l are no longer UFL indices after being used as iteration variables in the Python loops. You can re-define them by running

i, j, k, l = indices(4)

again before using them as indices to define B.

To format multi-line code blocks, you can type three backticks in a row before and after the code, like

```
line1
line2
etc.
```

(For more comprehensive formatting info, look up “Markdown”.)

2 Likes
7

Excellent! Thanks! It works now :smile: . Thanks for the pointing at “Markdown”.

8

Hi Kamensky,

thank you for your share. So, indices() function needs ufl import.
I just tried a part of the code you posted, using dolfin and ufl import.

from dolfin import *
from ufl import *

mesh = UnitSquareMesh(2,2)
i,j,k,l = indices(4)

# Tensor from StackExchange thread:
d = mesh.geometry().dim()
delta = Identity(d)
I = as_tensor(0.5*(delta[i,k]*delta[j,l] + delta[i,l]*delta[j,k]),(i,j,k,l))
f = Constant(0.0)

The indices() function is work, but the Constant function is error:

---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
<ipython-input-3-c87f29c3fa26> in <module>
      9 delta = Identity(d)
     10 I = as_tensor(0.5*(delta[i,k]*delta[j,l] + delta[i,l]*delta[j,k]),(i,j,k,l))
---> 11 f = Constant(0.0)

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/coefficient.py in Constant(domain, count)
    123 def Constant(domain, count=None):
    124     """UFL value: Represents a globally constant scalar valued coefficient."""
--> 125     domain = as_domain(domain)
    126     element = FiniteElement("Real", domain.ufl_cell(), 0)
    127     fs = FunctionSpace(domain, element)

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/domain.py in as_domain(domain)
    294         # closer to the user interface?
    295         # TODO: Make this configurable to be an error from the dolfin side?
--> 296         cell = as_cell(domain)
    297         return default_domain(cell)
    298 

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/cell.py in as_cell(cell)
    327         return TensorProductCell(cell)
    328     else:
--> 329         error("Invalid cell %s." % cell)

~/anaconda3/envs/fenicsproject/lib/python3.7/site-packages/ufl/log.py in error(self, *message)
    170         "Write error message and raise an exception."
    171         self._log.error(*message)
--> 172         raise self._exception_type(self._format_raw(*message))
    173 
    174     def begin(self, *message):

UFLException: Invalid cell 0.0.

Do I miss somthing, or is this a bug?

The code works well if I just import: from ufl import indices.

Thank you.

9

It looks like ufl.Constant and dolfin.Constant are defined differently. This sort of situation is a limitation of wildcard (i.e. *) imports, and why some people discourage using wildcard imports altogether. However, I find it convenient and aesthetically-pleasing in the case of dolfin, since the code looks more like mathematical notation without “dolfin.” everywhere (and the alternative of individually importing every mathematical function you need would get tedious).

2 Likes