Solving a 0D PDE system as BC for a 1D problem

Hi everyone,
I need to solve iteratively some coupled equations at the boundaries of my 1D model (surface/bulk diffusion model).

I know how to write the variational formulation for the bulk problem and for the boundary, but I don’t know how to “connect” the two. Especially, as there is diffusion in the BC equations, I need to rely on the existing mesh to evaluate the gradient.

It’s also unclear to me which functionspace to use : the concentration is 1D but I’ll be using some other unknowns that are going to be 0D.

Does anyone has any experience with this ?

Hi, you could define a mixed function space involving a 1D unknown u and scalar unknown v for the boundary using a “Real” element, something like:

from dolfin import *

mesh = UnitIntervalMesh(10)
Ve = FiniteElement("CG", mesh.ufl_cell(), 1)
Re = FiniteElement("Real", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]))
w = Function(W)
u, v = split(w)

Then it is possible to write things like u.dx(0)+v in a form

Thank you for the quick reply, so far it seems to work!

I have another question regarding u.dx(0). What exactly does this do ? It would make more sense to me to write something like u.ds(1) or u.ds(2) depending on the surface I’m working on (but I realise this does not work).

No, u.dx(0) is \partial u/\partial x.
But you can mark the extremities of the interval as 1 and 2 using something such as

def left(x, on_boundary):
    return near(x[0], 0)
def right(x, on_boundary):
    return near(x[0], 1)
facets = MeshFunction("size_t", mesh, 0) # 0 being the topological dimension of the facet (points here)
AutoSubDomain(left).mark(facets, 1)
AutoSubDomain(right).mark(facets, 2)
ds = Measure("ds", subdomain_data=facets)

and do a form such as expr*ds(1) or expr*ds(2)

Thank you for your answer, that is what I meant by ds(1) and ds(2) although I realize now I didn’t include the facets marking.

Yet another question. Now that my functionspace and my variational formulation are working, what kind of boundary condition am I supposed to include? I have a hard time grasping the concept for a scalar problem but I guess that not setting any is not the right way either (as the solver “finishes” with 0 iterations, although something is bound to happen as I have non-zero initial values for my scalar quantities).

I am sorry but it is hard to answer without knowing the actual problem you want to solve. Besides, you are speaking about boundary conditions and initial conditions which are not the same.
For boundary conditions, I imagine that they are included in your variational formulation since they are coupled to boundary unknowns.
For initial conditions, you just need to assign non-zero initial values to your mixed function (possibly 0 for the u component and a non-zero scalar to the v component), you can look at FunctionAssigner for instance.
https://fenicsproject.org/qa/9444/interpolation-on-mixed-functionspace-shapes-not-matching/?show=9447#a9447
https://fenicsproject.org/qa/8876/interpolate-on-mixedfunctionspace/