Partial integration

Errors
1

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)
2

See, for example https://fenicsproject.org/qa/11306/computing-projection-of-correlation-matrix/.

3

Hi, thanks for your answer. I’m new to Fenics so your example was unfortunately not very enlightening. Can you maybe explain why “f*dx(1)” does not work?

4

See Form language — Unified Form Language (UFL) 2021.1.0 documentation

6

Ok I wrote an UserExpression for rho as follows. It works.

class f_x(UserExpression):
    def __init__(self, x):
        UserExpression.__init__(self)
        self.x = x
    def eval(self, value, y):
        value[0] = f((self.x[0],y[1]))
    def value_shape(self):
        return ()

class Rho(UserExpression):
    def eval(self, value, x):
        fx = f_x(x)
        value[0] = assemble(fx*dx)
    def value_shape(self):
        return ()

rho = Rho()

But now the solve function is stuck on “Solving nonlinear variational problem.” for any time I’ve let the program run. My thinking is that this implementation is very slow. Do you think there’s any possible quicker implementation?