Neumann BC special function

variational formulation
1

Hi, I have a Neumann boundary condition expressed by an integral function (using scipy.quad) of coordinates, however I am not able to define it as a part of the variational form, since the function accept only float argument. Can anybody help?

Here’s an MWE:

import dolfin

import numpy as np

from scipy import integrate

mesh = dolfin.UnitIntervalMesh(10)

# Marking of the facet markers:
boundary = dolfin.MeshFunction("size_t",
                               mesh,
                               mesh.topology().dim()-1,
                               0) # this enforces default 0 value

# Mark boundaries
for f in dolfin.facets(mesh):

    # Mark inner boundary:
    if dolfin.near(f.midpoint()[0], 0.0):

        boundary[f] = 1 # left

    # Mark outer boundary:
    elif dolfin.near(f.midpoint()[0], 1.0):

        boundary[f] = 2 # right

T = dolfin.FunctionSpace(mesh, dolfin.FiniteElement("CG", "interval", 1))

ds = dolfin.Measure("ds", subdomain_data = boundary)

n = dolfin.FacetNormal(mesh)

x = dolfin.SpatialCoordinate(mesh)

q_int = lambda z, u: np.exp(-z*u)

q = lambda z: integrate.quad(lambda u: q_int(z, u), 0.0, np.inf)[0]

_theta, theta_ = dolfin.TrialFunction(T), dolfin.TestFunction(T)

a = dolfin.inner(dolfin.grad(_theta), dolfin.grad(theta_))*dolfin.dx

L = q(x[0])*theta_*ds(1)

bc = dolfin.DirichletBC(T, 1.0, boundary, 2)

theta = dolfin.Function(T)

dolfin.solve(a == L, theta, bc)