Calculate 4rank inverse tensor using ufl

Dear dolfinx users and developer.

First of all, I would like to express my deepest gratitude for your assistance and very good software.

Below is my code for inverse tensor of 4-rank (for ufl expression, put the tensor function in tensorinv. e.g., tensorinv(tensor) ).
But it takes a lot of time to calculate the inverse tensor.

If it is possible, would you let me know how to quickly derive inverse tensor of 4-rank?

Very thanks and I wish you have a nice day!

def flat(tensor):
    """
    Recursively flattens a ListTensor into a list of its elements.
    Expands 'Zero' objects to match their dimension.
    """
    from itertools import product
    flattened_list = []
    tensor_shape=tensor.ufl_shape
    index_combinations = product(*[range(d) for d in tensor_shape])
    for index in index_combinations:
        element = tensor[index]
        flattened_list.append(element)
    return flattened_list

def reshape(tensor, new_shape):
    """
    Reshapes a ufl.ListTensor to a new shape.
    
    Args:
        tensor (ufl.ListTensor): Input ListTensor to be reshaped.
        new_shape (tuple): New shape as a tuple of integers.
        
    Returns:
        ufl.ListTensor: Reshaped ListTensor.
    """
    # Flatten the original ListTensor into a single list

    flattened = flat(tensor)
    reshaped=np.reshape(flattened, new_shape).tolist()

    # Return a new ListTensor
    return reshaped

def determinant(matrix):
    n = len(matrix)
    if n == 1:
        return matrix[0][0]
    if n == 2:
        return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
    
    det = 0.0
    for col in range(n):
        sub_matrix = [[matrix[i][j] for j in range(n) if j != col] for i in range(1, n)]
        cofactor = (-1.0) ** col * matrix[0][col] * determinant(sub_matrix)
        det += cofactor
    return det

def cofactor_matrix(matrix):
    n = len(matrix)
    cofactor_mat = [[0] * n for _ in range(n)]
    for i in range(n):
        for j in range(n):
            sub_matrix = [
                [matrix[row][col] for col in range(n) if col != j]
                for row in range(n) if row != i
            ]
            cofactor_mat[i][j] = (-1) ** (i + j) * determinant(sub_matrix)
    return cofactor_mat

def transpose(matrix):
    return [list(row) for row in zip(*matrix)]

def inverse_matrix(matrix):
    det = determinant(matrix)
    if det == 0:
        return "there is no inverse matrix"
    
    cofactor_mat = cofactor_matrix(matrix)
    adjugate = transpose(cofactor_mat)
    
    n = len(matrix)
    inverse = [[adjugate[i][j] / det for j in range(n)] for i in range(n)]
    return ufl.as_tensor(inverse)

def tensorinv(tensor, ind=2):
    # tensor=ufl.as_tensor(np.arange(1, 3*3*3*3 + 1).reshape(3, 3, 3, 3).tolist())
    oldshape = tensor.ufl_shape
    prod = 1
    if ind > 0:
        invshape = oldshape[ind:] + oldshape[:ind]
        for k in oldshape[ind:]:
            prod *= k
    else:
        raise ValueError("Invalid ind argument.")
    tensor = reshape(tensor, (prod, prod))
    tensor_inverse = inverse_matrix(tensor)
    return ufl.as_tensor(reshape(tensor_inverse, (invshape)))

If there is no alternative, I would appreciate it if you could let me know how to implement a numerical local solver for inverting a 9x9 matrix (which represents the function values at the quadrature points)

I would probably use the dolfinx external operator by @a-latyshev GitHub - a-latyshev/dolfinx-external-operator: Extension of DOLFINx implementing the concept of external operator
For inverting the 9x9 matrix at quadrature points you could then use numpy or something along those lines.

1 Like

Very thanks dokken!! I will try what you recommended!