dolfin team mailing list archive

[Martin Sandve Alnæs <martinal@xxxxxxxxx>]

On Wed, Mar 4, 2009 at 7:47 AM, Johan Hake <hake@xxxxxxxxx> wrote:
> 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.

If I understand this right, the actual problem isn't the pointer comparison,
but that the mesh is read twice. That shouldn't be necessary to do.

Martin


> 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)
>> ----------------------------------------------------
>
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>