### Error due to Inverse

118
views
0
3 months ago by
Hello everyone,

In my code, I'm trying to get the inverse of an (6,6) tensor with  $inv\left(A\right)$inv(A) and then the kernel gives the error:

"The kernel appears to have died. It will restart automatically."

And Then I got an error message as:

adj_expr not implemented for dimension 6.​
What can be the problem ?

Here is the sample code:

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import warnings
%matplotlib inline

E=210000 #Young's Modulus
Poisson=0.3 #Poisson's Ratio
mu = E/(2*(1+Poisson))
kappa   = E/(3*(1-2*Poisson))
lame=Poisson*E/((1+Poisson)*(1-2*Poisson))

L = 20; W = 10 #Box Mesh, Length and Width

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W),10, 5,5)
V_tensor=VectorFunctionSpace(mesh,"CG",1,dim=6)
Fth = TensorFunctionSpace(mesh, "CG",1, shape=(6,6),symmetry=True)

I = as_vector([1,1,1,0,0,0])
I = project(I,V_tensor)

C_e = as_matrix([[2*mu+lame, lame, lame ,0,0,0],
[lame, 2*mu+lame, lame,0,0,0],
[lame, lame, 2*mu+lame,0,0,0],
[0, 0, 0,mu,0,0],
[0, 0, 0,0,mu,0],
[0, 0, 0,0,0,mu]])

P = project(C_e,Fth)

a=project(dot(inv(P),I),V_tensor)​

Regards,
Community: FEniCS Project
Please, provide a minimum working example that reproduces the error.
written 3 months ago by Eleonora Piersanti
1
Updated the code....
written 3 months ago by Christian
The code you added does not reproduce the error you mentioned. I get a different error:
NameError: name 'inv' is not defined
​
Please be sure that the code you are providing is really minimal (for example remove %matplotlib) and reproduces the error that you refer to so that we can help you fixing it! :)
written 3 months ago by Eleonora Piersanti
I tried it before typing here actually. It gives the following error for me : (v 2017.1.0)

adj_expr not implemented for dimension 6.​​
written 3 months ago by Christian
Ok I performed the block matrix method and now getting following error.. I do not know the meaning of the number given in the error as 15625. I decreased the mesh up to nearly 8-9 elements but the integration point number did not change in the error.
WARNING: The number of integration points for each cell will be: 15625
Consider using the option 'quadrature_degree' to reduce the number of points
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py in jit(ufl_object, form_compiler_parameters, mpi_comm)
141     try:
--> 142         result = ffc.jit(ufl_object, parameters=p)
143     except Exception as e:

/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py in jit(ufl_object, parameters, indirect)
203     # Inspect cache and generate+build if necessary
--> 204     module = jit_build(ufl_object, module_name, parameters)
205

/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py in jit_build(ufl_object, module_name, parameters)
119                                     params=params,
--> 120                                     generate=jit_generate)
121     return module

/usr/local/lib/python3.5/dist-packages/dijitso/jit.py in jit(jitable, name, params, generate, send, receive, wait)
159             # 1) Generate source code
--> 160             header, source, dependencies = generate(jitable, name, signature, params["generator"])
161             # Ensure we got unicode from generate

/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py in jit_generate(ufl_object, module_name, signature, parameters)
65     code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
---> 66             prefix=module_name, parameters=parameters, jit=True)
67

/usr/local/lib/python3.5/dist-packages/ffc/compiler.py in compile_form(forms, object_names, prefix, parameters, jit)
140     return compile_ufl_objects(forms, "form", object_names,
--> 141                                prefix, parameters, jit)
142

/usr/local/lib/python3.5/dist-packages/ffc/compiler.py in compile_ufl_objects(ufl_objects, kind, object_names, prefix, parameters, jit)
187     cpu_time = time()
--> 188     ir = compute_ir(analysis, prefix, parameters, jit)
189     _print_timing(2, time() - cpu_time)

/usr/local/lib/python3.5/dist-packages/ffc/representation.py in compute_ir(analysis, prefix, parameters, jit)
189     irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
--> 190            for (form_id, fd) in enumerate(form_datas)]
191     ir_integrals = list(chain(*irs))

/usr/local/lib/python3.5/dist-packages/ffc/representation.py in <listcomp>(.0)
189     irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
--> 190            for (form_id, fd) in enumerate(form_datas)]
191     ir_integrals = list(chain(*irs))

/usr/local/lib/python3.5/dist-packages/ffc/representation.py in _compute_integral_ir(form_data, form_id, prefix, element_numbers, classnames, parameters, jit)
452                                    classnames,
--> 453                                    parameters)
454

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/uflacsrepresentation.py in compute_integral_ir(itg_data, form_data, form_id, element_numbers, classnames, parameters)
--> 139                                 parameters)
140     ir.update(uflacs_ir)

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/build_uflacs_ir.py in build_uflacs_ir(cell, integral_type, entitytype, integrands, tensor_shape, coefficient_numbering, quadrature_rules, parameters)
352                 ir["unique_tables"], p["enable_table_zero_compression"],
--> 353                 rtol=p["table_rtol"], atol=p["table_atol"])
354

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py in build_optimized_tables(num_points, quadrature_rules, cell, integral_type, entitytype, modified_terminals, existing_tables, compress_zeros, rtol, atol)
565     # Analyze tables for properties useful for optimization
--> 566     unique_table_ttypes = analyse_table_types(unique_tables, rtol=rtol, atol=atol)
567

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py in analyse_table_types(unique_tables, rtol, atol)
544     return { uname: analyse_table_type(table, rtol=rtol, atol=atol)
--> 545              for uname, table in unique_tables.items() }
546

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py in <dictcomp>(.0)
544     return { uname: analyse_table_type(table, rtol=rtol, atol=atol)
--> 545              for uname, table in unique_tables.items() }
546

/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py in analyse_table_type(table, rtol, atol)
517         ttype = "ones"
--> 518     elif is_quadrature_table(table, rtol=rtol, atol=atol):
519         # Identity matrix mapping points to dofs (separately on each entity)

493     num_entities, num_points, num_dofs = table.shape
--> 494     I = numpy.eye(num_points)
495     return (num_points == num_dofs

/usr/lib/python3/dist-packages/numpy/lib/twodim_base.py in eye(N, M, k, dtype)
232         M = N
--> 233     m = zeros((N, M), dtype=dtype)
234     if k >= M:

MemoryError:

During handling of the above exception, another exception occurred:

RuntimeError                              Traceback (most recent call last)
<ipython-input-1-4bc524112f77> in <module>()
356         print ("incremental stress is",inc_stress(px,py,pz))
357         #Stress n+1 from C_ep
--> 358         strain333=project(dot(inv_C(Block11,Block22,Block33,Block44),inc_stress),V_tensor)
359         print("strain33333333 is",strain3333(px,py,pz))
360         stress2 = project(inc_stress+stress1,V_tensor)

/usr/lib/python3/dist-packages/dolfin/fem/projection.py in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
140     # Assemble linear system
141     A, b = assemble_system(a, L, bcs=bcs,
--> 142                            form_compiler_parameters=form_compiler_parameters)
143
144     # Solve linear system for projection

/usr/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble_system(A_form, b_form, bcs, x0, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend)
349     # Create dolfin Form objects referencing all data needed by assembler
350     A_dolfin_form = _create_dolfin_form(A_form, form_compiler_parameters)
--> 351     b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
352
353     # Create tensors

/usr/lib/python3/dist-packages/dolfin/fem/assembling.py in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
65         return Form(form,
66                     form_compiler_parameters=form_compiler_parameters,
---> 67                     function_spaces=function_spaces)
68     else:
69         raise TypeError("Invalid form type %s" % (type(form),))

/usr/lib/python3/dist-packages/dolfin/fem/form.py in __init__(self, form, form_compiler_parameters, function_spaces)
87         # Jit the module, this will take some time...
88         jit_result = jit(form, form_compiler_parameters,
---> 89                          mpi_comm=mesh.mpi_comm())
90         if jit_result is None:
91             cpp.dolfin_error("form.py",

/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py in mpi_jit(*args, **kwargs)
68         # Just call JIT compiler when running in serial
69         if MPI.size(mpi_comm) == 1:
---> 70             return local_jit(*args, **kwargs)
71
72         # Default status (0 == ok, 1 == fail)

/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py in jit(ufl_object, form_compiler_parameters, mpi_comm)
145         cpp.dolfin_error("jit.py",
146                          "perform just-in-time compilation of form",
--> 147                          "ffc.jit failed with message:\n%s" % (tb_text,))
148
149     # Unpack result (ffc.jit returns different tuples based on input type)

2837
2838     """
-> 2839     return _common.dolfin_error(location, task, reason)
2840
2841 def deprecation(feature: 'std::string', version_deprecated: 'std::string', message: 'std::string') -> "void":

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to perform just-in-time compilation of form.
*** Reason:  ffc.jit failed with message:
Traceback (most recent call last):
File "/usr/lib/python3/dist-packages/dolfin/compilemodules/jit.py", line 142, in jit
result = ffc.jit(ufl_object, parameters=p)
File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 204, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 120, in jit_build
generate=jit_generate)
File "/usr/local/lib/python3.5/dist-packages/dijitso/jit.py", line 160, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/local/lib/python3.5/dist-packages/ffc/jitcompiler.py", line 66, in jit_generate
prefix=module_name, parameters=parameters, jit=True)
File "/usr/local/lib/python3.5/dist-packages/ffc/compiler.py", line 141, in compile_form
prefix, parameters, jit)
File "/usr/local/lib/python3.5/dist-packages/ffc/compiler.py", line 188, in compile_ufl_objects
ir = compute_ir(analysis, prefix, parameters, jit)
File "/usr/local/lib/python3.5/dist-packages/ffc/representation.py", line 190, in compute_ir
for (form_id, fd) in enumerate(form_datas)]
File "/usr/local/lib/python3.5/dist-packages/ffc/representation.py", line 190, in <listcomp>
for (form_id, fd) in enumerate(form_datas)]
File "/usr/local/lib/python3.5/dist-packages/ffc/representation.py", line 453, in _compute_integral_ir
parameters)
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/uflacsrepresentation.py", line 139, in compute_integral_ir
parameters)
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/build_uflacs_ir.py", line 353, in build_uflacs_ir
rtol=p["table_rtol"], atol=p["table_atol"])
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py", line 566, in build_optimized_tables
unique_table_ttypes = analyse_table_types(unique_tables, rtol=rtol, atol=atol)
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py", line 545, in analyse_table_types
for uname, table in unique_tables.items() }
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py", line 545, in <dictcomp>
for uname, table in unique_tables.items() }
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py", line 518, in analyse_table_type
File "/usr/local/lib/python3.5/dist-packages/ffc/uflacs/elementtables.py", line 494, in is_quadrature_table
I = numpy.eye(num_points)
File "/usr/lib/python3/dist-packages/numpy/lib/twodim_base.py", line 233, in eye
m = zeros((N, M), dtype=dtype)
MemoryError
.
*** Where:   This error was encountered inside jit.py.
*** Process: 0
***
*** DOLFIN version: 2017.1.0
*** Git changeset:  7215d3707edbf33a51e7b4bc3dfc8d67922651a2​
written 3 months ago by Christian
written 3 months ago by Jan Blechta
Integration is performed by quadrature. Your integrand probably has a terribly high polynomial degree. No surprise when using Cramer's rule which blows up quickly. You possibly want to decrease the quadrature order - search the reference manual for quadrature degree.
written 3 months ago by Jan Blechta
Right, it is incredibly huge.. I was using quad degree as 3, the error did not change when I decrease quad degree. Resultant inverse matrix also in (6,6) shape but includes thousands of polynomial degree on it. Is not there a method to simplyfy it ?
written 3 months ago by Christian
You must have done a mistake. When using quadrature degree 3 there can't be 15625 quadrature points.
written 3 months ago by Jan Blechta
The inverse function is in shape (6,6) and goes like below : (It is only small portion of the function, it did not include all...)
ListTensor(ListTensor(Indexed(Inverse(Sum(ListTensor(ListTensor(Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(0), FixedIndex(0)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(0), FixedIndex(1)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(0), FixedIndex(2))))), ListTensor(Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(1), FixedIndex(0)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(1), FixedIndex(1)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(1), FixedIndex(2))))), ListTensor(Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(2), FixedIndex(0)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(2), FixedIndex(1)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(2), FixedIndex(2)))))), ComponentTensor(Product(IntValue(-1), Indexed(ComponentTensor(IndexSum(Product(Indexed(ListTensor(ListTensor(Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 383)), MultiIndex((FixedIndex(3), FixedIndex(0)))), Indexed(Conditional(LE(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), FiniteElement('Discontinuous Lagrange', tetrahedron, 0)), 244), FloatValue(1e-10)), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), TensorElement(FiniteElement('Lagrange', tetrahedron, 1), shape=(6, 6), symmetry={(3, 2): (2, 3), (4, 3): (3, 4), (3, 0): (0, 3), (5, 2): (2, 5), (3, 1): (1, 3), (5, 4): (4, 5), (2, 1): (1, 2), (2, 0): (0, 2), (5, 0): (0, 5), (5, 1): (1, 5), (4, 2): (2, 4), (1, 0): (0, 1), (5, 3): (3, 5), (4, 1): (1, 4), (4, 0): (0, 4)})), 368), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 2), ​
written 3 months ago by Christian
I don't understand what's going on. Maybe you could start a new question. It seems to me that your original problem has been resolved.
written 3 months ago by Jan Blechta

3
3 months ago by
The error message is pretty self-explaining. The inverse is computed using Cramer's rule. That's not implemented for 6x6 matrix in UFL. The implementation exists for 2x2, 3x4 and 4x4 matrix: https://bitbucket.org/fenics-project/ufl/src/085616b79ffb8fa0b9cf20b682e9f2928a8c7b44/ufl/compound_expressions.py#compound_expressions.py-156.

Considering that the matrix involved is block-diagonal with blocks 3x3, you can implement the inversion for your matrix by inverting the 3x3 blocks.
Is there any procedure or predefined function to convert it to block-diagonal with 3x3 blocks ?
written 3 months ago by Christian
1
Something along the lines of

block = as_tensor([[P[i,j] for i in range(3)] for j in range(3)])​
written 3 months ago by Jan Blechta