ufl_legacy.log.UFLException: determinant_expr not implemented for shape (4, 4)

Linear Algebra
1

I want to assemble the matrices of my weak form, I’m doing compressible flow and plastic fluid structure interaction, the variables are (p,u,T), which have 4 variables, then I write a function to assemble the weak form

from fenics import *
import numpy as np
from ufl_legacy import transpose

from Compute_DC import compute_DC_constants
from compute_SUPG import compute_tau_supg

from fenics import *
import numpy as np

def assemble_weak_form(V_fluid, Y, Y_n, dt, A0_f, A1_f, A2_f, A1_sp, A2_sp,
K11, K12, K21, K22, A0_s, A1_s, A2_s,
markers, beta=1.0):
“”"
Assemble the weak form for FSI problem with 4x4 matrices
“”"
# Test and trial functions
W = TestFunction(V_fluid)
dY = TrialFunction(V_fluid) # For bilinear form

# Get mesh info
mesh = V_fluid.mesh()
n = FacetNormal(mesh)
h = CellDiameter(mesh)

# Define measures
dx = Measure('dx', domain=mesh)
dx_s = Measure('dx', subdomain_data=markers, subdomain_id=2)
dS = Measure('dS', domain=mesh)

# Initialize bilinear and linear forms
a = 0  # Bilinear form
L = 0  # Linear form

# Get gradients
grad_dY = grad(dY)
grad_W = grad(W)
grad_Y = grad(Y)

# 1. Mass terms
# LHS (bilinear form)
a += inner(W, dot(A0_f, dY)) * dx  # Whole domain
a -= inner(W, dot(A0_f, dY)) * dx_s  # Subtract solid domain
a += inner(W, dot(A0_s, dY)) * dx_s  # Add solid domain mass

# RHS (linear form)
L += (1/dt)*inner(W, dot(A0_f, Y_n)) * dx  # Whole domain
L -= (1/dt)*inner(W, dot(A0_f, Y_n)) * dx_s  # Subtract solid domain
L += (1/dt)*inner(W, dot(A0_s, Y_n)) * dx_s  # Add solid domain mass

# 2. Convective terms
# LHS
conv_term = dot(A1_f, grad_dY[:, 0]) + dot(A2_f, grad_dY[:, 1])
a += inner(W, conv_term) * dx  # Whole domain
a -= inner(W, conv_term) * dx_s  # Subtract solid domain

# 3. Diffusive terms
# LHS
F_diff_x = dot(K11, grad_dY[:, 0]) + dot(K12, grad_dY[:, 1])
F_diff_y = dot(K21, grad_dY[:, 0]) + dot(K22, grad_dY[:, 1])

a -= inner(grad_W[:, 0], F_diff_x) * dx
a -= inner(grad_W[:, 1], F_diff_y) * dx
a += inner(grad_W[:, 0], F_diff_x) * dx_s
a += inner(grad_W[:, 1], F_diff_y) * dx_s

# 4. SUPG stabilization
# Compute residual for current solution (RHS)
F_diff_x_n = dot(K11, grad_Y[:, 0]) + dot(K12, grad_Y[:, 1])
F_diff_y_n = dot(K21, grad_Y[:, 0]) + dot(K22, grad_Y[:, 1])
div_F_diff = grad(F_diff_x_n)[:, 0] + grad(F_diff_y_n)[:, 1]

Res_n = dot(A0_f, (Y - Y_n) / dt) + \
        dot(A1_f, grad_Y[:, 0]) + dot(A2_f, grad_Y[:, 1]) - \
        div_F_diff

# LHS SUPG terms
tau_SUPG = compute_tau_supg(dt, A1_sp, A2_sp, K11, K12, K21, K22, mesh)

# Linearized residual for LHS
Res_lin = dot(A0_f, dY / dt) + conv_term - (grad(F_diff_x)[:, 0] + grad(F_diff_y)[:, 1])

B_st = inner(dot(transpose(A1_f), grad_W[:, 0]) +
             dot(transpose(A2_f), grad_W[:, 1]),
             dot(tau_SUPG, Res_lin)) * dx
B_st_s = inner(dot(transpose(A1_f), grad_W[:, 0]) +
               dot(transpose(A2_f), grad_W[:, 1]),
               dot(tau_SUPG, Res_lin)) * dx_s

a += B_st - B_st_s

# 5. DC stabilization
K_DC = compute_DC_constants(h, Res_n, Y)
a += inner(grad_W, dot(K_DC, grad_dY)) * dx
a -= inner(grad_W, dot(K_DC, grad_dY)) * dx_s

# 6. Solid domain terms
a += inner(W, dot(A1_s, grad_dY[:, 0]) + dot(A2_s, grad_dY[:, 1])) * dx_s

# 7. Interface terms
a += beta * h * inner(jump(dot(grad_W, n)), jump(dot(grad(dY), n))) * dS

# Assemble matrices
A = assemble(a)
b = assemble(L)

return A, b

here, the K11 and other K_ij matrices are 44, and all matrices here are 44, then when I want to assemble the matrices, it says that

determinant_expr not implemented for shape (4, 4).
Traceback (most recent call last):
  File "/home/a/Desktop/FluidStructure/main.py", line 258, in <module>
    a, L = assemble_weak_form(V_fluid, solution_fluid, solution_fluid_n, dt,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/a/Desktop/FluidStructure/WeakForm.py", line 105, in assemble_weak_form
    A = assemble(a)
        ^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
    return Form(form,
           ^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
                                            ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
                ^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/compute_form_data.py", line 253, in compute_form_data
    form = apply_algebra_lowering(form)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_algebra_lowering.py", line 175, in apply_algebra_lowering
    return map_integrand_dags(LowerCompoundAlgebra(), expr)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 46, in map_integrand_dags
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 27, in map_integrands
    mapped_integrals = [map_integrands(function, itg, only_integral_type)
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 35, in map_integrands
    return itg.reconstruct(function(itg.integrand()))
                           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/map_integrands.py", line 46, in <lambda>
    return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 99, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/apply_algebra_lowering.py", line 135, in inverse
    return inverse_expr(A)
           ^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/compound_expressions.py", line 143, in inverse_expr
    return adj_expr(A) / determinant_expr(A)
                         ^^^^^^^^^^^^^^^^^^^
  File "/usr/lib/python3/dist-packages/ufl_legacy/compound_expressions.py", line 101, in determinant_expr
    error("determinant_expr not implemented for shape %s." % (sh,))
  File "/usr/lib/python3/dist-packages/ufl_legacy/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl_legacy.log.UFLException: determinant_expr not implemented for shape (4, 4).

how to solve it

2

You would need to implement the determinant by hand. I would do it by lowering the determinant to the sum of 3x3 determinants, and use ufl.det on those. Ref Form language — Unified Form Language (UFL) 2021.1.0 documentation (section Taking components of tensors)

3

Could you please tell me how to transform this weak form to matrices? I encounter this error while doing this, should I change the source code?

4

Your code is not properly formatted, and I do not have the bandwidth to reformat it and reduce it to a minimal example.

Please reduce it to a code where you have as few terms as possible and still reproduce the error.

From what I can tell, calling assemble(a) and remove term by term from a should tell you what term is problematic.