dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

Hello Marc!

This is a known bug. The two read in meshes are of course the same but during 
assemble, we extract one mesh from the same form. To be able to check that 
the mesh are the same we just compare the mesh pointers. Line 79 in 
dolfin/fem/Form.cpp. This is really not what we want, as your case will break 
this. 

We have previously discussed this issue, can't find the email now. Then we 
thought of including some easily compiled mesh signature that could be used 
in the comparison.

I do not know why this worked before the FunctionSpace rework.

A possible work-around for you before we have worked this out could be:

**********************
# read in Velocity and pressure functions (from a Taylor-Hood Stokes problem)
p = Function("pressure.xml.gz")
v0 = Function("velocity.xml.gz")

# extract the mesh
mesh = p.function_space().mesh()

V = VectorFunctionSpace(mesh,"CG",2)

# Then you interpolate
v = interpolate(v0,V)
# or just
v = Function(V)
v.vector().assign(v0.vector())

**********************

If this work: ;)
We had a drag-lift demo that had the same problem. Before we released 9.0 we 
removed the part that did not work. Should we re-add it for the instructive 
cause? 

Johan

On Wednesday 04 March 2009 05:00:02 mspieg wrote:
> 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)
> ----------------------------------------------------