Hi all,
I've recently been having some issues when trying to solve with
forms where I attach multiple discrete functions that share the same
mesh but not necessarily the same function_space (i.e. element).
This used to work under 0.8.1 and I suspect it's a bug.
A reasonably simple example of the problem can be stated as follows.
let V, p be discrete Velocity and pressure functions generated from a
Stoke's solve with Taylor-hood mixed elements. The two solutions
share the same mesh, but are in different function spaces (V is P2, p
is P1).
I now want to form a new vector valued function
v = V - grad(p)
by projecting V and grad(p) onto an appropriate function space. (If
you're wondering, this is a toy version of a deformable porous media
problem)
Taking what seems to be the obvious approach
----------------------------------------------
# read in Velocity and pressure functions (from a Taylor-Hood Stokes
problem)
p = Function("pressure.xml.gz")
V = Function("velocity.xml.gz")
# extract the mesh
mesh = p.function_space().mesh()
# define a new function space for v on the same mesh
BDM = FunctionSpace(mesh, "BDM", 1)
# Define projection problem
w = TestFunction(BDM)
u = TrialFunction(BDM)
a = dot(w,u)*dx
L = dot(w,V)*dx - dot(w,grad(p))*dx
# Compute solution
problem = VariationalProblem(a, L)
v = problem.solve()
-----------------------------------------------------------
generates the following error from solve (JIT/FFC seems to work fine)
Solving linear variational problem
Traceback (most recent call last):
File "test.py", line 44, in <module>
v = problem.solve()
File "/Users/mspieg/mspieg/hg_repositories/FEniCS/srcFEniCS/dolfin/
local/lib/python2.5/site-packages/dolfin/variationalproblem.py", line
52, in solve
cpp.VariationalProblem.solve(self, u)
File "/Users/mspieg/mspieg/hg_repositories/FEniCS/srcFEniCS/dolfin/
local/lib/python2.5/site-packages/dolfin/cpp.py", line 10742, in solve
return _cpp.VariationalProblem_solve(*args)
RuntimeError: *** Error: Unable to extract mesh from form
(nonmatching meshes for function spaces).
If I just project the pressure part of v i.e.
L = -dot(w,grad(p))
everything works fine.
Similarly if I just project the velocity part (and extract the mesh
from V rather than p) i.e.
mesh = V.function_space().mesh()
...
L = dot(w,V)
that also works fine.
If i use the mesh from p to project V, things break and vice-versa
(and yet it should be the same mesh. And they certainly are if you
write them out and diff the .xml files)
I can obviously project each part separately and add them together,
but that seems like a hack and I don't see why this shouldn't work.
I also see this behaviour in my c++ codes as well and have tracked it
to problems where I try to attach multiple discrete functions on
different function spaces as coefficients to a form. Like I said,
this used to work in 0.8.1, but breaks now.
In general, I don't see why you can't have arbitrary discrete
functions in a linear form as long as they can be evaluated during
assembly.
apologies for the long e-mail but all help greatly appreciated
cheers
marc
p.s. I've attached a test script with the above problem and variations
p.p.s versions should all be up to date development versions
Dolfin 0.9.1 dev
ufc 1.1.1 dev
FFC 0.6.1 dev
instant 0.9.6 dev (changset 283)

----------------------------------------------------
Marc Spiegelman
Lamont-Doherty Earth Observatory
Dept. of Applied Physics/Applied Math
Columbia University
http://www.ldeo.columbia.edu/~mspieg
tel: 845 704 2323 (SkypeIn)
----------------------------------------------------
