DOLFIN

assemble error: All terms in form must have same rank

Asked by Sebastian Rau on 2012-05-30

Hello,

i am facing the "All terms in form must have the same rank" error while trying to assamble a form. Here is a minimal code:

meshsize = 32
mesh = UnitSquare(meshsize,meshsize)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
Q = VectorFunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V,Q,V])

v = TestFunction(W)

class State(object):
def __init__(self):
  self.uk = Function(W)
  self.u0 = Function(W)
  self.E = Function(V)
  self.Vext = Function(V)
  self.problem = None

def solve(self, u0, E, massn0):
  # self.E.interpolate(E)
  self.u0 = u0
  self.massn0 = massn0
  self.uk = Function(W, self.u0)
  self.du = TrialFunction(W)

  self.u00, self.u01, self.u02 = split(self.u0)
  self.uk0, self.uk1, self.uk2 = split(self.uk)
  self.du0, self.du1, self.du2 = split(self.du)
  v0, v1, v2 = split(v)

  self.J = (4.0/dt)*(self.uk0*self.du0 - self.du0*self.u00)*v0*dx \
   - 2.0*inner( div(outer((outer(self.uk1,self.uk1)), self.du0/self.uk0**3)), v1 )*dx \
   - 2.0*( self.uk0*self.du0*(self.uk2 + self.E*Vext[0]) + self.du0*lapl(self.uk0) + self.uk0*lapl(self.du0) )*div( v1 )*dx \
   - 2.0*inner( outer(grad(self.uk0*self.du0),( self.uk2 + self.E*Vext[0] - eps**2*lapl(self.uk0)/self.uk0 )), v1 )*dx \
   + eps**2*inner( outer(grad(self.uk0**2), (lapl(self.du0)*self.uk0 - lapl(self.uk0)*self.du0)/(self.uk0**2) ), v1 )*dx \
   - 2.0*self.uk0*self.du0*dx \
   - inner(self.du1, grad(v0))*dx + 1.0/dt*inner(self.uk1,v1)*dx + inner( div( 2*outer(outer(self.uk1, self.du1),(1/self.uk0**2)) ), v1)*dx \
   - self.uk0**2*self.du2*div(v1)*dx - inner( outer(grad(self.uk0**2),self.du2),v1 )*dx + lamda**2*inner( grad(self.du2),grad(v2) )*dx

self.du = Function(W)
self.A = assemble(self.J)

The error occurs in the last line, thus the mistake should be somewhere in the definition of J, but i can't find it! I also played a bit with Test and Trial Functions for self.du, but it doesnt help either...

Any help is appreciated!

Regards
Sebastian

Question information

Language:: English Edit question

Status:: Answered

For:: DOLFIN Edit question

Assignee:: No assignee Edit question

Last query:: 2012-05-30

Last reply:: 2012-05-30

Link existing bug

Revision history for this message

Martin Sandve Alnæs (martinal) said on 2012-05-30:

This is not a minimal code. Here's how to reduce it to find the error:
Set self.J to only the first ...*dx term and run. For example you can
do like this to quickly disable the other lines:

self.J = (4.0/dt)*(self.uk0*self.du0 - self.du0*self.u00)*v0*dx
"""\
   - 2.0*inner( div(outer((outer(self.uk1,self.uk1)), self.du0/self.uk0**3)), v1 )*dx \
   - 2.0*( self.uk0*self.du0*(self.uk2 + self.E*Vext[0]) + self.du0*lapl(self.uk0) + self.uk0*lapl(self.du0) )*div( v1 )*dx \
   - 2.0*inner( outer(grad(self.uk0*self.du0),( self.uk2 + self.E*Vext[0] - eps**2*lapl(self.uk0)/self.uk0 )), v1 )*dx \
   + eps**2*inner( outer(grad(self.uk0**2), (lapl(self.du0)*self.uk0 - lapl(self.uk0)*self.du0)/(self.uk0**2) ), v1 )*dx \
   - 2.0*self.uk0*self.du0*dx \
   - inner(self.du1, grad(v0))*dx + 1.0/dt*inner(self.uk1,v1)*dx + inner( div( 2*outer(outer(self.uk1, self.du1),(1/self.uk0**2)) ), v1)*dx \
   - self.uk0**2*self.du2*div(v1)*dx - inner( outer(grad(self.uk0**2),self.du2),v1 )*dx + lamda**2*inner( grad(self.du2),grad(v2) )*dx
"""

Until you get the error, add one term at a time and run. When you get the error, you've found the term. Then you can probably see the error yourself.

Revision history for this message

Kristian B. Ølgaard (k.b.oelgaard) said on 2012-05-30:

On 30 May 2012 14:25, Martin Sandve Alnæs
<email address hidden> wrote:
> Question #198874 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/198874
>
> Status: Open => Answered
>
> Martin Sandve Alnæs proposed the following answer:
> This is not a minimal code. Here's how to reduce it to find the error:
> Set self.J to only the first ...*dx term and run. For example you can
> do like this to quickly disable the other lines:
>
> self.J = (4.0/dt)*(self.uk0*self.du0 - self.du0*self.u00)*v0*dx
> """\
> - 2.0*inner( div(outer((outer(self.uk1,self.uk1)), self.du0/self.uk0**3)), v1 )*dx \
> - 2.0*( self.uk0*self.du0*(self.uk2 + self.E*Vext[0]) + self.du0*lapl(self.uk0) + self.uk0*lapl(self.du0) )*div( v1 )*dx \
> - 2.0*inner( outer(grad(self.uk0*self.du0),( self.uk2 + self.E*Vext[0] - eps**2*lapl(self.uk0)/self.uk0 )), v1 )*dx \
> + eps**2*inner( outer(grad(self.uk0**2), (lapl(self.du0)*self.uk0 - lapl(self.uk0)*self.du0)/(self.uk0**2) ), v1 )*dx \
> - 2.0*self.uk0*self.du0*dx \
> - inner(self.du1, grad(v0))*dx + 1.0/dt*inner(self.uk1,v1)*dx + inner( div( 2*outer(outer(self.uk1, self.du1),(1/self.uk0**2)) ), v1)*dx \
> - self.uk0**2*self.du2*div(v1)*dx - inner( outer(grad(self.uk0**2),self.du2),v1 )*dx + lamda**2*inner( grad(self.du2),grad(v2) )*dx
> """
>
> Until you get the error, add one term at a time and run. When you get
> the error, you've found the term. Then you can probably see the error
> yourself.