Hybrid FEM formulation (with some non-spatial DOFs) with FeniCS

Question

I would like to implement a hybrid finite element formulation, that has, in addition to normal (mixed) functionspace, several additional degrees of freedom that are related to the boundary condition. These additional DOFs form their own block(s) in FEM matrix and vectors, that are connected to the ordinary DOFs that correspond to (vector) field on the mesh. I would like to know, if this is at all possible currently with FeniCS?

For additional details of the hybrid method, see https://jyx.jyu.fi/dspace/handle/123456789/36645

Answer 1

You can use Real spaces for for the ODE dofs. Even though the Real spaces were meant to be used as, for example, a global Lagrange multiplier, they can be used to formulate global ODEs which can be coupled to any spatial entity in your domain.

Have in mind that the Real dofs are not well integrated into DOLFIN. They are by default added to each element matrix and assemble with many Real dofs can therefore take unnecessary long time. This can be avoided using dedicated block structures for these dofs, but that is not supported by DOLFIN (yet). You need to use external softwares like cbc.block for that purposes.

Comments

Thanks for quick answer. Can you suggest me some simple example code? Which classes should I use?

---

The following simple script solves a diffusion equation coupled to a simple growth ode on the boundary:

from dolfin import *

mesh = UnitSquareMesh(10,10)
right = AutoSubDomain(lambda x, on_boundary: on_boundary and near(x[0], 1.0) and 0.25 <= x[1] <= 0.75)
left = AutoSubDomain(lambda x, on_boundary: on_boundary and near(x[0], 0.0))
boundaries = FacetFunction("size_t", mesh, 0)
right.mark(boundaries, 1)
left.mark(boundaries, 2)
dss = ds[boundaries]

# Function spaces
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
VR = V*R

# Solution functions
U = Function(VR)
U0 = Function(VR)
u, p = split(U)
u0, p0 = split(U0)

v, q = TestFunctions(VR)

dt = Constant(0.1)
t_stop = 10.0

bc = DirichletBC(VR.sub(0), Constant(0), boundaries, 2)

# Area needed to scale "volume" integral
area = Constant(1.0)
roof = Constant(100)

# Diffusion equation with p-depend flux condition
F = ((u-u0)*v/dt + inner(grad(u), grad(v)))*dx - p*v*dss(1)


# Simple growth equation for the p variable
F += ((p-p0)/dt - (roof-u))/area*q*dx

# Initial conditions
t = 0.0
U.vector()[:] = 1.0
U0.vector()[:] = 1.0
u_plot = Function(V)
while t < t_stop:
    solve(F==0, U, bc)
    U0.assign(U)

    u_plot.assign(U.split(deepcopy=1)[0])
    plot(u_plot)
    t += float(dt)

plot(u_plot, interactive = True)

---

Ok I believe there has been some misunderstanding.

I wish I could have a number of real valued DOFs that are not in relation with any spatial node. If I understand correctly, in your example, there is a real-valued DOF for each spatial node (that's why it's related to mesh). What I need, is a real dof that is not related to mesh. I.e. something like this:

V = FunctionSpace(mesh, "CG", 1)
R1 = SingleRealDOF() #no mesh parameter here
R2 = SingleRealDOF()
...
VR = MixedFunctionSpace([V,R1,R2,R3,...RN])
U = Function(VR)
u, r1,r2,r3...rn = split(U)
...

and the relation to other (spatial) DOFS is then built in weak form as follows:

F += somecoefficientsur1ds(1) + someothercoefficientsur2ds(1) ...

(my actual application of interest is found in the paper I cited, eq. (25)).
Then, the after solving with solve(F==0,U,bc), we would obtain field u and single real value r as follows:

u,r1,r2,r3,...,rn = split(U). realdofs r1...rn would not be a field over mesh, but just a single real values, such that

print type(r1)
< type 'float' >

But there doesn't seem to be any way to do like this? Or is there?

---

I wish I could have a number of real valued DOFs that are not in relation with any spatial node.

That's exactly what the 'R' function space is, i.e. one global DOF. Dependence on a mesh is only an implementational issue.

---

Jan is correct. However it is important to notice that a geometrical aspect is introduced to the weak formulation whenever you integrate a global dof. That is why I have to divide by the area of the whole domain in:

F += ((p-p0)/dt - (roof-u))/area*q*dx

as dx is an area integral over the whole domain.

If your mesh is a 3D mesh and is very big, you can choose an integration domain which is much smaller, for example some few triangles on the boundary. The above line could then be:

area_of_subdomain_2 = some_value
F += ((p-p0)/dt - (roof-u))/area_of_subdomain_2*q*ds(2)

---

I have now been able to solve my problem, thanks! As you (Johan) suggested, I have used cbc.block and it makes assemble a lot faster. However, it's pretty slow still. Block structure in my case is as follows:

  [[G,     Hhat, Khat],
  [Htilde, H,   0 ],
  [Ktilde, 0,   K ]]

So assemble times for *hat and *tilde matrice blocks are pretty slow. Still some ideas that could help me get the assemble function quicker?

Or if there is no other way I could look forward to edit the dolfin package itself, to make it faster. Any estimate how complicated this could be?