dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Wednesday 04 March 2009 09:23:10 Garth N. Wells wrote:
> Johan Hake wrote:
> > On Wednesday 04 March 2009 08:09:00 Martin Sandve Alnæs wrote:
> >> 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.
> >
> > Well, that is an another related bug, which Anders mention in his answer.
> > When the issue of reading in Funcitons from file is solved this problem
> > will probably not be an issue anymore.
>
> This will be resolved by not reading in Functions (at least not
> directly). In the future, it will be required to read the mesh and the
> vector (and possibly the dof map), and the the user should supply the
> finite element.

Ok.

> > However I still think it would be neat, maybe not feasable, to compare
> > the equality of two meshes based on some signature, not only a pointer.
>
> We should try to avoid using signatures - it introduces various
> difficulties.

Ok.

Johan

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