### How to Make a Custom Tensor Expression

197
views
1
5 months ago by
I am trying to make 2x2 tensor expression, but it is not working. How can I make a tensor expression?
from fenics import *
class MarkerTest(Expression):
def eval(self, values, x):
values[0,0] = 1.0
values[1,0] = 2.0
values[0,1] = 3.0
values[1,1] = 4.0
def value_shape(self):
return (2,2)

marker = MarkerTest(degree=1)
print(marker(0))​

Here is the type of error message encountered:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input> in <module>()
10
11 marker = MarkerTest(degree=1)
---> 12 print(marker(0))

/usr/lib/python3/dist-packages/dolfin/functions/expression.py in __call__(self, *args, **kwargs)
863
864         # The actual evaluation
--> 865         self.eval(values, x)
866
867         # If scalar return statement, return scalar value.

<ipython-input> in eval(self, values, x)
2 class MarkerTest(Expression):
3     def eval(self, values, x):
----> 4         values[0,0] = 1.0
5         values[1,0] = 2.0
6         values[0,1] = 3.0

IndexError: too many indices for array

Community: FEniCS Project

3
5 months ago by
As the error message says, there ar too many indices. For an array, there should be only one:
class MarkerTest(Expression):
def eval(self, values, x):
values[0] = 1.0
values[1] = 2.0
values[2] = 3.0
values[3] = 4.0
def value_shape(self):
return (2,2)​

So would the tensor be as follows?

\begin{bmatrix}
values[0] & values[2]  \\
values[1] & values[3] \\
\end{bmatrix}
written 5 months ago by Saro
It is the inverse, i.e.,
$\begin{bmatrix} \mathrm{values}[0] & \mathrm{values}[1] \\ \mathrm{values}[2] & \mathrm{values}[3] \end{bmatrix}$
You can check it with something like
mesh = UnitSquareMesh(1, 1)
n1 = Constant(("1.0", "0.0"))
n2 = Constant(("0.0", "1.0"))

print(assemble(dot(n1, marker*n1)*dx(domain=mesh)))
print(assemble(dot(n1, marker*n2)*dx(domain=mesh)))
print(assemble(dot(n2, marker*n1)*dx(domain=mesh)))
print(assemble(dot(n2, marker*n2)*dx(domain=mesh)))​
written 5 months ago by Adam Janecka
Thank you. That solves my problem.
written 5 months ago by Saro