dolfin team mailing list archive

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

On Wed, Mar 4, 2009 at 8:35 AM, Anders Logg <logg@xxxxxxxxx> wrote:
> On Wed, Mar 04, 2009 at 08:09:00AM +0100, 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.
>>
>> Martin
>
> The reason it even works at all (to read in the two functions) is that
> DOLFIN happens to have precompiled elements for the cases in question.
> If it had been cubics, it wouldn't have worked.
>
> For this reason, we've discussed removing the precompiled elements.
> One would then (in C++) first need to create the FunctionSpace and
> then read the degrees of freedom for the Function. The problem with
> this is that the mesh needs to be read first, so instead of just
>
>  Function u("function.xml");
>
> we get
>
>  Mesh mesh("mesh.xml");
>  MyFunctionSpace V(mesh);
>  Function u(V);
>  File file("udofs.xml");
>  file >> u;
>
> We could add a constructor to Function to simplify this to
>
>  Mesh mesh("mesh.xml");
>  MyFunctionSpace V(mesh);
>  Function u(V, "udofs.xml");
>
> The file udofs.xml would contain all the dofs of u stored as a Vector.

This looks good. In any but the most  trivial of applications, you'll need to
make the mesh and V anyway. We should have both approaches though,
because it's natural to load a file into some existing Function object.

Martin


> In Python, it would be possible to handle storing of general functions
> to file. The XML file could contain everything needed to dynamically
> instantiate the function space at run-time (no need for precompiled
> elements). But then the interfaces would be very different in C++ and
> Python for reading/writing Functions.
>
> --
> Anders
>
>
>>
>> > 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
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkmuLzsACgkQTuwUCDsYZdF9+ACfTbJGul/Bb7xCynDwWLM64542
> gFUAoIWhAVIz5Mllw8hHvRwJ5P/yCzZT
> =aUBt
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>
>