Assemble forms from mixed space formulation

I have solved the Helmholtz equation with complex “wave number” in dolfin without problems by splitting it in its real and imaginary parts and using mixed function spaces.

The problem along with the code is already shown here, but I copy here the code, just in case:

from dolfin import  RectangleMesh, MeshFunction, SubDomain, Measure, FunctionSpace, MixedFunctionSpace, Constant
from dolfin import TrialFunctions, TestFunctions, inner, grad, dx, Function, solve, plot, Point, near, Expression
from dolfin import assemble

#geometrical constants and mesh generation
Lx = 0.01
aspect_ratio = 100
Ly = aspect_ratio*Lx
Nx = 5
Ny = aspect_ratio*Nx 

mesh = RectangleMesh( Point(0,0), Point(Lx,Ly), Nx, Ny )

dim_edge = 1
def_value = 0

boundary_markers = MeshFunction("size_t", mesh, dim_edge, def_value)

tol = 1.0E-14

front = 1
back = 2

class Gamma_front(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0, tol)

GF = Gamma_front()
GF.mark(boundary_markers, front)

class Gamma_back(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], Lx, tol)

GB = Gamma_back()
GB.mark(boundary_markers, back)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

#Physical constants
kappa = 200
h_conv = 15 
epsilon = 0.8
sigma = 5.67E-8
T_air = 298
rho = 2700
c= 900
omega = 0.8

alpha = 0.4

I = 600*(1 + 0.0j)

sx = -15*Lx
sy = 0.5*Ly

Q_R = Expression("alpha*I_R/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_R = I.real, sx=sx, sy=sy)
Q_I = Expression("alpha*I_I/(2*pi)*(x[0]-sx)/((x[0]-sx)*(x[0]-sx) + (x[1]-sy)*(x[1]-sy))", degree=2, alpha = alpha, I_I = I.imag, sx=sx, sy=sy)

# Funtion spaces 
P2 = FunctionSpace(mesh, "Lagrange", 1 )
H = MixedFunctionSpace( P2, P2  ) 
Phi_R, Phi_I = TrialFunctions(H)
Psi_R, Psi_I = TestFunctions(H)

# variational forms
a =  kappa*( inner( grad(Phi_R), grad(Psi_R)) + inner( grad(Phi_I), grad(Psi_I) ))*dx     \
                                          - rho*omega*c*( Phi_R*Psi_I -  Phi_I*Psi_R)*dx  \
                                          + h_eff*( Phi_R*Psi_R + Phi_I*Psi_I)*ds(front)  \
                                          + h_eff*( Phi_R*Psi_R + Phi_I*Psi_I)*ds(back)
L = ( Q_R*Psi_R + Q_I*Psi_I )*ds(front)

#problem solution
T = Function(H)
solve(a==L, T)
T_R, T_I = T.split()

However, now I was interested in solving similar problems with different right hand sides, so I tried to compute the assembled matrix of the bilinear form using:

from dolfin import assemble
A = assemble(a)
b = assemble(L)

and I get the message error: ArityMismatch: Adding expressions with non-matching form arguments ('v_0^0', 'v_1^1') vs ('v_0^1', 'v_1^0').

After performing a look on the internet I saw that this can be caused by the forms not having all the trial and test functions of the mixed space, but I don’t understand how to fix it.

Thank you very much

Are you sure that you want to use MixedFunctionSpace, and not
H = FunctionSpace(mesh, MixedElement([P2, P2]))?
As far as im aware, MixedFunctionSpace is used for @cdaversin’s mixed dimensional FE implementation, and not to created a MixedElement on the same meshes

I’m sure I was not understanding what I was doing. I thought that both things were equivalent. I will rewrite it as you suggest. Is there any place where I can read about the difference between both?


Rewritten, works now!

See for instance;