Finding eigenvalues of a 3x3 tensor in dolfinx

dolfinx
1

Hello

I have been using dolfin to find the eigenvalues of a 2x2 ufl tensor and it’s been working great. I now need to use dolfinx for a 3x3 tensor as the calculation involves complex numbers.

In dolfin I use the following piece of code to determine whether I have the right eigenvalues in 2D:

from dolfin import *

def e1(A): 
     return ((tr(A) + sqrt(tr(A)**2-4*det(A)))/2)

def e2(A): 
     return ((tr(A) - sqrt(tr(A)**2-4*det(A)))/2)

S = FunctionSpace(mesh, 'CG', 2)
Ca = Expression((('5','2'),('2', '1')), degree=0)
C = project(e1(Ca), S)
vertex_values = C.compute_vertex_values(mesh)
for i, x in enumerate(mesh.coordinates()):
    print(x, C(x))

Here is the code I have written in dolfinx to find the eigenvalues of a 3D tensor. Can someone tell me where I’m going wrong?

from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import numpy as np
import ufl

C = ufl.as_matrix([[-2, -4, 2], [-2, 1, 2], [4, 2, 5]])

a = -1
b = ufl.tr(C)
c = -0.5*((ufl.tr(C))**2-ufl.tr(C*C))
d = ufl.det(C)

del_0 = b**2 - 3*a*c
del_1 = 2*(b**3) - 9*a*b*c + 27*(a**2)*d

C_pos = (0.5*(del_1+ufl.sqrt (del_1**2-4*del_0**3) ))**(1/3)

prim = 0.5*(-1+ufl.sqrt(-3))

x0 = (-1/3*a)*(b+C_pos+(del_0/C_pos))
x1 = (-1/3*a)*(b+prim*C_pos+(del_0/(prim*C_pos)))
x2 = (-1/3*a)*(b+(prim**2)*C_pos+(del_0/((prim**2)*C_pos)))

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = fem.FunctionSpace(domain, ("CG", 1))
uD = fem.Function(V, dtype=np.complex128)

uD.interpolate(fem.Expression(ufl.real(x0), V.element.interpolation_points, dtype=np.complex128))

The corresponding error message is:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [2], in <cell line: 32>()
     29 V = fem.FunctionSpace(domain, ("CG", 1))
     30 uD = fem.Function(V, dtype=np.complex128)
---> 32 uD.interpolate(fem.Expression(ufl.real(x0), V.element.interpolation_points, dtype=np.complex128))

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/function.py:98, in Expression.__init__(self, ufl_expression, X, form_compiler_params, jit_params, dtype)
     68 def __init__(self, ufl_expression: ufl.core.expr.Expr, X: np.ndarray,
     69              form_compiler_params: dict = {}, jit_params: dict = {},
     70              dtype=PETSc.ScalarType):
     71     """Create DOLFINx Expression.
     72 
     73     Represents a mathematical expression evaluated at a pre-defined
   (...)
     95 
     96     """
---> 98     assert X.ndim < 3
     99     num_points = X.shape[0] if X.ndim == 2 else 1
    100     _X = np.reshape(X, (num_points, -1))

AttributeError: 'builtin_function_or_method' object has no attribute 'ndim'

Thanks in advance,

Will.

2

Hi,
Following my question and Dokken’s answer in another topic, it looks like you are missing a pair of parenthesis in

uD.interpolate(fem.Expression(ufl.real(x0), V.element.interpolation_points, dtype=np.complex128))

Just write V.element.interpolation_points() and it should work.

3

Hi Pamgur,

That change results in the following new error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [8], in <cell line: 29>()
     26 V = fem.FunctionSpace(domain, ("CG", 1))
     27 uD = fem.Function(V, dtype=np.complex128)
---> 29 uD.interpolate(fem.Expression(ufl.real(x0), V.element.interpolation_points(), dtype=np.complex128))

File /usr/local/dolfinx-complex/lib/python3.10/dist-packages/dolfinx/fem/function.py:102, in Expression.__init__(self, ufl_expression, X, form_compiler_params, jit_params, dtype)
     99 num_points = X.shape[0] if X.ndim == 2 else 1
    100 _X = np.reshape(X, (num_points, -1))
--> 102 mesh = ufl_expression.ufl_domain().ufl_cargo()
    104 # Compile UFL expression with JIT
    105 if dtype == np.float32:

AttributeError: 'NoneType' object has no attribute 'ufl_cargo'
4

I managed to resolve the problem myself. If anyone in the future needs to find eigenvalues of 3x3 tensors, here is my code:

from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import numpy as np
import ufl
from petsc4py.PETSc import ScalarType

domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [1, 1, 1]], [5, 5, 5], mesh.CellType.hexahedron)
C = ufl.as_matrix([[fem.Constant(domain, ScalarType(-2)), fem.Constant(domain, ScalarType(-4)), fem.Constant(domain, ScalarType(2))],
                   [fem.Constant(domain, ScalarType(-2)), fem.Constant(domain, ScalarType(1)), fem.Constant(domain, ScalarType(2))],
                   [fem.Constant(domain, ScalarType(4)), fem.Constant(domain, ScalarType(2)), fem.Constant(domain, ScalarType(5))]])
a = -1
b = ufl.tr(C)
c = -0.5*((ufl.tr(C))**2-ufl.tr(C*C))
d = ufl.det(C)

del_0 = b**2 - 3*a*c
del_1 = 2*(b**3) - 9*a*b*c + 27*(a**2)*d

C_pos = (0.5*(del_1+ufl.sqrt (del_1**2-4*del_0**3) ))**(1/3)

V = fem.FunctionSpace(domain, ("CG", 1))

prim = 0.5*(-1+fem.Constant(domain, ScalarType(ufl.sqrt(-3))))

x0 = (-1/3*a)*(b+C_pos+(del_0/C_pos))
x1 = (-1/3*a)*(b+prim*C_pos+(del_0/(prim*C_pos)))
x2 = (-1/3*a)*(b+(prim**2)*C_pos+(del_0/((prim**2)*C_pos)))

domain = mesh.create_unit_square(MPI.COMM_WORLD, 5, 5, mesh.CellType.quadrilateral)
expr = fem.Expression(ufl.real(x2), V.element.interpolation_points(), dtype=np.complex128)
uD = fem.Function(V)
uD.interpolate(expr)
print(uD.x.array)
2 Likes