User Expression broken in FEniCS 2018.1

[clason]

FEniCS 2018.1 changed how UserExpressions work (due to the switch to pybind11), see, e.g., https://www.allanswered.com/post/jrmxq/expression-with-cpp-code-behaviour-changed-with-2018-1-0-dev0pybind11/, https://www.allanswered.com/post/jrmxq/expression-with-cpp-code-behaviour-changed-with-2018-1-0-dev0pybind11/. It seems that this also impacts how Expressions are propagated in Julia: The following line from the README.md no longer works with FEniCS.jl (although it does work in Python and using PyCall!)

u_D = Expression("1+x[0]*x[0]+2*x[1]*x[1]", degree=2)

instead giving

ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:101 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'RuntimeError'>
RuntimeError('Must supply C++ code to Expression. You may want to use UserExpression',)
  File "/opt/miniconda/envs/fenicsproject/lib/python3.6/site-packages/dolfin/function/expression.py", line 369, in __init__
    raise RuntimeError("Must supply C++ code to Expression. You may want to use UserExpression")

Stacktrace:
 [1] pyerr_check at /home/clason/.julia/packages/PyCall/WcrLS/src/exception.jl:60 [inlined]
 [2] pyerr_check at /home/clason/.julia/packages/PyCall/WcrLS/src/exception.jl:64 [inlined]
 [3] macro expansion at /home/clason/.julia/packages/PyCall/WcrLS/src/exception.jl:84 [inlined]
 [4] __pycall!(::PyCall.PyObject, ::Ptr{PyCall.PyObject_struct}, ::PyCall.PyObject, ::PyCall.PyObject) at /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:101
 [5] _pycall! at /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:62 [inlined]
 [6] macro expansion at /home/clason/.julia/packages/PyCall/WcrLS/src/exception.jl:84 [inlined]
 [7] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{}, ::Int64, ::PyCall.PyObject) at /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:28
 [8] _pycall!(::PyCall.PyObject, ::PyCall.PyObject, ::Tuple{}, ::Base.Iterators.Pairs{Symbol,Any,NTuple{8,Symbol},NamedTuple{(:cppcode, :element, :cell, :domain, :degree, :name, :label, :mpi_comm),Tuple{String,Nothing,Nothing,Nothing,Int64,Nothing,Nothing,Nothing}}}) at /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:16
 [9] #call#89 at /home/clason/.julia/packages/PyCall/WcrLS/src/pyfncall.jl:147 [inlined]
 [10] PyObject at ./none:0 [inlined]
 [11] #Expression#13(::Nothing, ::Nothing, ::Nothing, ::Int64, ::Nothing, ::Nothing, ::Nothing, ::Type, ::String) at /home/clason/.julia/packages/FEniCS/MA1xs/src/jfem.jl:74
 [12] (::getfield(Core, Symbol("#kw#Type")))(::NamedTuple{(:degree,),Tuple{Int64}}, ::Type{Expression}, ::String) at ./none:0
 [13] top-level scope at none:0

It would also help to prominently note in the README which FEniCS version this package is compatible with (and maybe test the version when using FEniCS and warn if it's not the right one).

[tungli]

I had the same problem. Actually, for me there were two problems. One problem was related to the Conda version of dolfin where mpi_comm was missing from keyword arguments. The other is that cppcode is not a keyword. Therefore, line 74 in jfem.jl is wrong I think.

[ysimillides]

Hi, Ill check the cppcode line in jfem.jl over the weekend and let you know.

[EvergreenTree]

Same issue here.

Calling FFC just-in-time (JIT) compiler, this may take some time.
PyError ($(Expr(:escape, :(ccall(#= /Users/evergreen/.julia/packages/PyCall/0jMpb/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'RuntimeError'>
RuntimeError('Must supply C++ code to Expression. You may want to use UserExpression',)
  File "/anaconda3/lib/python3.6/site-packages/dolfin/function/expression.py", line 369, in __init__
    raise RuntimeError("Must supply C++ code to Expression. You may want to use UserExpression")


Stacktrace:
 [1] pyerr_check at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/exception.jl:60 [inlined]
 [2] pyerr_check at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/exception.jl:64 [inlined]
 [3] macro expansion at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/exception.jl:84 [inlined]
 [4] __pycall!(::PyObject, ::Ptr{PyCall.PyObject_struct}, ::PyObject, ::PyObject) at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/pyfncall.jl:44
 [5] _pycall!(::PyObject, ::PyObject, ::Tuple{}, ::Int64, ::PyObject) at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/pyfncall.jl:29
 [6] _pycall!(::PyObject, ::PyObject, ::Tuple{}, ::Base.Iterators.Pairs{Symbol,Any,NTuple{8,Symbol},NamedTuple{(:cppcode, :element, :cell, :domain, :degree, :name, :label, :mpi_comm),Tuple{String,Nothing,Nothing,Nothing,Int64,Nothing,Nothing,Nothing}}}) at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/pyfncall.jl:11
 [7] #call#89 at /Users/evergreen/.julia/packages/PyCall/0jMpb/src/pyfncall.jl:89 [inlined]
 [8] PyObject at ./none:0 [inlined]
 [9] #Expression#13(::Nothing, ::Nothing, ::Nothing, ::Int64, ::Nothing, ::Nothing, ::Nothing, ::Type, ::String) at /Users/evergreen/.julia/packages/FEniCS/MA1xs/src/jfem.jl:74
 [10] (::getfield(Core, Symbol("#kw#Type")))(::NamedTuple{(:degree,),Tuple{Int64}}, ::Type{Expression}, ::String) at ./none:0
 [11] top-level scope at In[5]:3

[EvergreenTree]

Seems that this project is no longer maintained? I decided to use python instead. Alternatively, one can simply use pycall to run a python version of fenics inside Julia. However, this results in a pyobject return, which is not compatible with FEniCS.jl.

[ysimillides]

@EvergreenTree I'll try to have a look into your issue in depth tomorrow when I get into the office. I was under the impression this issue was fixed with the code above. Are you on the latest version (i.e. master) of fenics.jl?

[EvergreenTree]

@ysimillides Thanks for your reply. I tried reinstalling both with Pkg.add("FEniCS") and with cloning this git, but got exactly the same error. I suspect that this is the same issue as @clason has said, i.e. the change of Expression function into UserExpression during one FEniCS update.

[EvergreenTree]

@ysimillides After looking into your code (line 74 of jfem.jl), I found that this works Expression(fenics.Expression("1+x[0]*x[0]+2*x[1]*x[1]", degree=2)) However this line in the tutorial still doesn't work:

> U = FEniCS.Function(V)
MethodError: no constructors have been defined for Function

Stacktrace:
 [1] top-level scope at In[159]:12

> methods(FEniCS.Function)
0 methods for generic function Type:

[EvergreenTree]

Yet this works: U = FEniCS.FeFunction(V) Now everything works fine for me. Hope it helps.

[EvergreenTree]

@ysimillides Anyway I think it's a bug. Could you please update this fix to the master or let me do this? I wish it would be helpful to FEniCS users. Thank you!

[clason]

@EvergreenTree I use this snippet to extract sparse Julia matrices from FEniCS: https://github.com/JuliaDiffEq/FEniCS.jl/issues/55#issuecomment-414434584 For PyCall 1.90, the snippet has to be slightly adapted (and a workaround seems to be no longer needed)

using SparseArrays,PyCall   

# create a demo matrix
fe = pyimport("fenics")
mesh = fe.UnitSquareMesh(64,64)
V = fe.FunctionSpace(mesh,"CG",1) 
u,v = fe.TrialFunction(V),fe.TestFunction(V)
a = fe.dot(fe.grad(u),fe.grad(v))*fe.dx + fe.Dx(u,1)*v*fe.dx # non-symmetric form
A_fe = fe.assemble(a) # wrapper for dolfin matrix

# convert to PETSc matrix, accessible via petsc4py
A_py = fe.as_backend_type(A_fe).mat()

# convert to Julia matrix -- note that we feed the CSR structure to a CSC constructor
m,n = A_py.getSize()
indptr,indices,data = A_py.getValuesCSR()
A_jl = SparseMatrixCSC(m,n,indptr.+1,indices.+1,data)
# transpose to get correct structure (collect lazy wrapper)
A_jl = sparse(A_jl') 

(Like you, I switched to direct PyCall access to FEniCS as this project seemed... less than fully maintained.)