dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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

Attachment: signature.asc
Description: Digital signature