Is there a way in Fenics to compute a partial integral? I need \rho(x)=\int f(x,y)dy, but I get a TypeError: “unsupported operand type(s) for *: ‘Expression’ and ‘Form’”. The problem is that the equation for c only involves the variable x but depends on \rho and not on f. In the following implementation I get the TypeError.

dolfin import * # Sub domain for Dirichlet boundary condition class NofluxBoundary(SubDomain): def inside(self, x, on_boundary): return bool((x[0] < DOLFIN_EPS or x[0] > (1.0 - DOLFIN_EPS)) \ and on_boundary) # Sub domain for Periodic boundary conditional class PeriodicBoundaryC(SubDomain): # Left boundary is "target domain" G def inside(self, x, on_boundary): return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary) # Map right boundary (H) to left boundary (G) def map(self, x, y): y[1] = x[1] - 1.0 y[0] = x[0] class PeriodicBoundary(SubDomain): # Left boundary is "target domain" G def inside(self, x, on_boundary): return bool(x[1] < DOLFIN_EPS and x[1] > -DOLFIN_EPS and on_boundary) # Create mesh and finite element mesh = UnitSquareMesh(32, 32) pbc = PeriodicBoundaryC() nofluxb = NofluxBoundary() pb = PeriodicBoundary() V = FunctionSpace(mesh, 'CG', 1, constrained_domain=pbc) Q = FiniteElement("CG", triangle, 1) domains = MeshFunction("size_t", mesh, mesh.topology().dim()) domains.set_all(0) boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) boundaries.set_all(0) pb.mark(boundaries, 1) nofluxb.mark(boundaries, 2) ds = Measure('ds', domain=mesh, subdomain_data=boundaries) W = FunctionSpace(mesh,MixedElement([Q,Q]), constrained_domain=pbc) #system u = Function(W) u.interpolate(Expression(("cos(0.5*pi*x[0])*cos(0.5*pi*x[0])", "cos(0.5*pi*x[0])*cos(0.5*pi*x[0])"), degree=2)) f, c = split(u) v, w = TestFunction(W) r = Expression('x[0]',degree=2) cospsi = Expression('cos(2*pi*x[1])',degree=2) sinpsi = Expression('sin(2*pi*x[1])',degree=2) #system B = r*cospsi*grad(f)[0]*v*dx - sinpsi*grad(f)[1]*v*dx \ + r*sinpsi*grad(c)[0]*f*grad(v)[1]*dx \ - r*grad(f)[0]*v*ds(1) \ + r*grad(c)[0]*grad(w)[0]*dx -r*(f*dx(1))*w*dx - \ - grad(c)[0]*w*ds(1) solve(B == 0, f)