dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

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.

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.

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