Form has no parts with arity 2.

407
views
2
8 months ago by
Dear all,

I was working on a bigger problem where I had a nonlinear system of PDE's. More on that can be found in How to manage discontinous coefficients.

I was using the regular Newton nonlinear solver. However, as it was mentioned on the referenced post, I could probably be better using a fixed iteration method.
Thus I wanted to assembly the equations to manually solve. However after assembling I am getting this output (it is not exactly an error, but an output):

"Form has no parts with arity 2."

My knowlodge of math is still basic thus I would really like to now why do I get this error. In my interpretation the form is the integrand that I constructed to represent the weak form. The arity concept for me is totally new, I found that it could be defined as the number of arguments of the function.
However, I thought that if I had terms that had all the related test and trial functions I would be in good shape.
For this problem, I made a small MWE where I think of solving the following abstract problem (I do not know if it is well posed at all, I just wanted to give an example).

I tried to put the following code and when tried to assemble had the same output.

---------------------------------------------------------------------------------------------
MWE example:
from fenics import *

mesh = IntervalMesh(100, 0, 1)

P1 = FiniteElement('CG', interval, 2)
V = FunctionSpace(mesh, MixedElement([P1, P1]))
u = Function(V)
v_1, v_2 = TestFunctions(V)
u_1, u_2 = split(u)
f_1 = Constant(1)
f_2 = Constant(5)
g_1 = Constant(1)
g_2 = Constant(5)

f_1*v_1*dx  + f_2*v_2*dx + g_1*v_1*ds + g_2*v_2*ds

a = lhs(F)
A = assemble(a)

L = rhs(F)
b = assemble(L)​

I know it is not a complicated question, I am just trying to learn so any awnser that is didatic will be of great value.
Thank you very much for the help, and if I am being too vague or lazy, please just tell me and point me the right direction to understand and solve this doubt.
All the best, Murilo

Community: FEniCS Project

4
8 months ago by
"lhs()" and "rhs()" split the form into bilinear and linear forms, respectively.  You have not used any trial functions, and so there is no linear form to split up.  You have started forming a problem for use with the nonlinear solver.

This'll work though:
from fenics import *

mesh = IntervalMesh(100, 0, 1)

P1 = FiniteElement('CG', interval, 2)
V  = FunctionSpace(mesh, MixedElement([P1, P1]))
v_1, v_2 = TestFunctions(V)
u_1, u_2 = TrialFunctions(V)
f_1 = Constant(1)
f_2 = Constant(5)
g_1 = Constant(1)
g_2 = Constant(5)

f_1*v_1*dx  + f_2*v_2*dx + g_1*v_1*ds + g_2*v_2*ds

a = lhs(F)
A = assemble(a)

L = rhs(F)
b = assemble(L)​

u = Function(V)
solve(A, u.vector(), b)

print u.vector().max(), u.vector().min()

edit:
Note, however, that you have a term
$\int_{\Omega} u_1 u_2\ \mathrm{d}\Omega$
that does not include any test functions and so I would expect something bad to happen here...

Here's an experiment:

from fenics import *

mesh = UnitIntervalMesh(2)

P1 = FiniteElement('CG', interval, 2)
V  = FunctionSpace(mesh, P1)

u = TrialFunction(V)
v = TestFunction(V)

# NOTE: using lhs(F) in the assembles below to avoid RuntimeError.
# these first two are equivalent, as they are scalars :
F = inner(u, u)*dx
A = assemble(lhs(F))   # Form has no parts with arity 2.

F = u * u * dx
A = assemble(lhs(F))   # Form has no parts with arity 2.

F = v * v * dx
A = assemble(lhs(F))   # Form has no parts with arity 2.

# these are OK :

# vector :
F = u*dx
a = assemble(F)
print 'one:\t', a.array()

# vector :
F = v*dx
a = assemble(F)
print 'two:\t', a.array()

# matrix :
F = inner(u, v)*dx
a = assemble(F)
print 'three:\n', a.array()

# different matrix :
a = assemble(F)
print 'four:\n', a.array()

# functions can be assembled too :
f = Function(V)
f.vector()[:] = 1.0

# scalar :
F = f * dx
a = assemble(F)
print 'five:\t', a

# vector :
F = f * v * dx
a = assemble(F)
print 'six:\t', a.array()

# matrix :
F = f * u * v * dx
a = assemble(F)
print 'seven:\n', a.array()
That is perfect pf4d!
However, one small last doubt, to solve the system, should I call a:

u = Function(V)
solve(A, u.vector(), b)​
? I mean, would this method solve for u_1 and u_2?
Thank you very much for the help!
written 8 months ago by Murilo Moreira
1
[ [ [ [ [ [ [ this is correct.  ] ] ] ] ] ]
written 8 months ago by pf4d
I was going to say that I tested and got it right! Thank you very much, and with the edit, it is way more clear, you rock!
written 8 months ago by Murilo Moreira
I was looking forward to correct the term that you rightly described as wrong because it was not multiplied by any test function. I tried to add a multiplication with v_1, however now I have the following:
"Multiplying expressions with overlapping form argument number 1, argument is v_1."
I need to search on this error now.
Did you ever saw any error like this?
Thank you very much!
written 8 months ago by Murilo Moreira