dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Wednesday 18 February 2009 21:58:29 A Navaei wrote:
> 2009/2/18 Johan Hake <hake@xxxxxxxxx>:
> > On Wednesday 18 February 2009 21:22:15 A Navaei wrote:
> >> 2009/2/18 Johan Hake <hake@xxxxxxxxx>:
> >> > On Wednesday 18 February 2009 20:56:59 A Navaei wrote:
> >> >> 2009/2/18 Anders Logg <logg@xxxxxxxxx>:
> >> >> > On Wed, Feb 18, 2009 at 02:04:35PM +0000, A Navaei wrote:
> >> >> >> 2009/2/17 Anders Logg <logg@xxxxxxxxx>:
> >> >> >> > On Tue, Feb 17, 2009 at 03:28:08PM +0000, Garth N. Wells wrote:
> >> >> >> >> A Navaei wrote:
> >> >> >> >> > The following minimal test for Function in c++ reveals some
> >> >> >> >> > bugs. I guess this example can help me with dealing with the
> >> >> >> >> > current issues of ImageFunction.
> >> >> >> >> >
> >> >> >> >> > (1) interpolate.py does not work when a Function is created
> >> >> >> >> > in c++ and wrapped (see comment [2]). It seems that the bug
> >> >> >> >> > is originated from the copy constructor (see comment [3])
> >> >> >> >> >
> >> >> >> >> > (2) In order to perform the interpolation, why is it
> >> >> >> >> > necessary to create another Function and then copy it?
> >> >> >> >> >
> >> >> >> >> > (3) Signature checkes seem not working properly (see comment
> >> >> >> >> > [1]). The signature-based assignments are error-prone anyway,
> >> >> >> >> > why the object-oriented approach is not used?
> >> >> >> >>
> >> >> >> >> Signatures are used to permit reading/writing Functions to a
> >> >> >> >> file. They are indeed error prone, so I believe that we reached
> >> >> >> >> a consensus a short while ago that we would remove pre-compiled
> >> >> >> >> elements.
> >> >> >> >>
> >> >> >> >> Garth
> >> >> >> >
> >> >> >> > Instead of signatures, I'd recommend that you define a simple
> >> >> >> > form file for each of the different types of FunctionSpace you
> >> >> >> > need, for example:
> >> >> >> >
> >> >> >> >  element = FiniteElement("CG", "triangle", 1)
> >> >> >> >
> >> >> >> >  v = TestFunction(element)
> >> >> >> >  u = TrialFunction(element)
> >> >> >> >  a = v*u*dx
> >> >> >> >
> >> >> >> > If you put this in a file named My.form and compile it with FFC
> >> >> >> > using -l dolfin, you will get a class named MyFunctionSpace that
> >> >> >> > you can then instantiate using just a mesh:
> >> >> >> >
> >> >> >> >  MyFunctionSpace V(mesh);
> >> >> >> >
> >> >> >> > Create one form file for each of the different types of
> >> >> >> > FunctionSpace that you need, name the files to something
> >> >> >> > suitable and use the generated code. That way you won't need to
> >> >> >> > worry about signatures, dofmaps and finite elements.
> >> >> >>
> >> >> >> Effectively, I've been using the very same method all this time,
> >> >> >> it does not work.
> >> >> >
> >> >> > Yes, it does. It's used in about 20 of the demos.
> >> >> >
> >> >> >> The copy constructor fix never worked. I've been trying to explain
> >> >> >> this in many different ways, but the right attention was never
> >> >> >> paid to this. Let's see if the sandbox example can convince you
> >> >> >> this time.
> >> >> >>
> >> >> >> A Function instance still cannot be returned by reference (or
> >> >> >> value). Returning as shared_ptr seems to work initially, but
> >> >> >> eventually it generates segmentation fault -- see attached.
> >> >> >
> >> >> > Yes, it can. There's absolutely no problem to return a Function by
> >> >> > reference. See the updated sandbox demo.
> >> >> >
> >> >> > The only problem is when you want to copy a Function which is only
> >> >> > defined in terms of an eval() operator. In those cases the Function
> >> >> > cannot be copied.
> >> >> >
> >> >> > If you do the following:
> >> >> >
> >> >> > class MyFunction : public Function
> >> >> > {
> >> >> > public:
> >> >> >
> >> >> >  MyFunction(const FunctionSpace& V) : Function(V) {};
> >> >> >
> >> >> >  void eval(double* values, const double* x) const
> >> >> >  {
> >> >> >    values[0] = sin(x[0]);
> >> >> >  }
> >> >> > };
> >> >> >
> >> >> > MyFunction f(V);
> >> >> > Function g = f;
> >> >> >
> >> >> > Do you then expect g to return sin(x)? It would be possible to
> >> >> > implement this but it would require g to keep a pointer to f so
> >> >> > that the eval() in g may call the eval() in f.
> >> >>
> >> >> Yes, we eventually want to get the image data in eval() and obviously
> >> >> the work around:
> >> >>
> >> >>  Function g(V);
> >> >>  g.vector();
> >> >>  _f = g;
> >> >
> >> > I haven't fully followed the thread, so I might say somthing that has
> >> > already been mentioned, but shouldn't the approach used in
> >> > interpolate.py work?
> >>
> >> I initially used the approach in interpolate.py, which didn't work:
> >>
> >>  Function g(V);
> >>  _f = g;
> >
> > Which interpolate.py are you talking about?
>
> dolfin/site-packages/interpolate.py

Ok, but here we do nothing like:

Function g(V);
_f = g;

> > The c++(ish) equivalent of interpolate.py would be:
> >
> >  Function interpolate(&Function g, V)
> >  {
> >     Function f(V);
> >     g.interpolate(f.vector(),V);
> >     return f;
> >  }
> >
> > Here f is a discrete function which will be returned by value.
>
> Just to mention that the discussion is not really about interpolation.
> It's about having a Function instance inside a class initialised by
> the input FunctionSpace, see the below attachment.

Yes but if you want to return a Function by value you need to copy it. 

To be able to copy it, the function need to be a discrete one, i.e., the 
vector need to be initialized. This overrides your eval function, as the 
values in the vector is used to evaluate the function. (This is what happens 
to you)

To be able to use the information given by the overloaded eval function you 
need to interpolate the values to a vector. This vector needs to be attached 
to the function you want to return, also making it a discrete one that you 
want. 

The point here is that you need to make an interpolation before you return the 
Function by value. You do not do that in the attached code. I added this in 
my previous post below.

> >>  Then Anders added the, not so obvious, addition of initialising the
> >> vector:
> >>
> >>  Function g(V);
> >>  g.vector();
> >>  _f = g;
> >>
> >> But this doesn't help with eval().
> >
> > Correct as no interpolation is done, just an initialization of the
> > vector.
>
> Interpolation is done, see the original attachment:
> http://www.fenics.org/pipermail/dolfin-dev/attachments/20090218/d53d755d/at
>tachment-0001.cpp

Yes but _after_ you have returned the function from FunctionContainer, right? 
You really need to do it before it is returned, see above.

There you try to interpolate the Function onto it self. The functions vector 
is initilized, otherwise you could not return it by value from 
FunctionContainer, but it is initialized with zeros making the interpolation

  f2.interpolate(f2.vector(), f2.function_space());

pointless, as you are interpolating a Function f2, using its eval, which only 
returns zero, as it is a discrete function with a zero vector.

Johan

>
> -Ali
>
> >> Something is going on in the python code, or the wrapper, which c++
> >> lacks.
> >
> > No this is essentially what happens in interpolate.py
> >
> > def interpolate(v, V):
> >    #...
> >    # Compute interpolation
> >    Pv = Function(V)
> >    v.interpolate(Pv.vector(), V)
> >
> >    return Pv
> >
> > which is exactly the same I suggest above for C++.
> >
> >> How did you end up with using g instead of directly using f?
> >
> > I want to interpolate the function g, using its eval, to the _f's
> > _vector, which defines the coefficients in the discrete function in the
> > FunctionSpace V. This is just how you do it and which also interpolate.py
> > does.
> >
> > Johan
> >
> >> -Ali
> >>
> >> >  class FunctionContainer
> >> >  {
> >> >  public:
> >> >    FunctionContainer(const FunctionSpace& V):_f(V)
> >> >    {
> >> >      message("assigning function");
> >> >      MyFunction g(V);
> >> >      g.interpolate(_f.vector(),V);
> >> >    };
> >> >
> >> >    const Function& get_function()
> >> >    {
> >> >      message("returning");
> >> >      return _f;
> >> >    };
> >> >  protected:
> >> >    Function _f;
> >> >  };
> >> >
> >> > Here g's vector is never initialized and the eval function is used to
> >> > interpolate to the vector of _f.
> >> >
> >> > Johan
> >> >
> >> >> does not call eval(). I don't think if it is possible to do this
> >> >> without amending the Function class? The denial of changing the
> >> >> visibility of the member variables to protected is making this
> >> >> unnecessarily more and more complicated. I am attaching the updated
> >> >> sandbox test.
> >> >>
> >> >>
> >> >> -Ali
> >> >>
> >> >> > --
> >> >> > Anders
> >> >> >
> >> >> >> -Ali
> >> >> >>
> >> >> >> > -----BEGIN PGP SIGNATURE-----
> >> >> >> > Version: GnuPG v1.4.9 (GNU/Linux)
> >> >> >> >
> >> >> >> > iEYEARECAAYFAkma2rQACgkQTuwUCDsYZdHp4ACfSbCXc2FAulzIdDsKvhz/6EGV
> >> >> >> > aY4An0eyftGV3hxR3L25M9LPu3X7KFg+
> >> >> >> > =z1cY
> >> >> >> > -----END PGP SIGNATURE-----
> >> >> >> >
> >> >> >> > _______________________________________________
> >> >> >> > DOLFIN-dev mailing list
> >> >> >> > DOLFIN-dev@xxxxxxxxxx
> >> >> >> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> >> >> >>
> >> >> >> // Place for random tests
> >> >> >>
> >> >> >> #include <dolfin.h>
> >> >> >> #include "Poisson.h"
> >> >> >>
> >> >> >> using namespace dolfin;
> >> >> >>
> >> >> >> class MyFunction : public Function
> >> >> >> {
> >> >> >> public:
> >> >> >>
> >> >> >>   MyFunction(const FunctionSpace& V) : Function(V) {};
> >> >> >>
> >> >> >>   void eval(double* values, const double* x) const
> >> >> >>   {
> >> >> >>     message("Calling eval");
> >> >> >>     double dx = x[0] - 0.5;
> >> >> >>     double dy = x[1] - 0.5;
> >> >> >>     values[0] = 500.0*exp(-(dx*dx + dy*dy) / 0.02);
> >> >> >>   }
> >> >> >> };
> >> >> >>
> >> >> >> class FunctionContainer
> >> >> >> {
> >> >> >> public:
> >> >> >>   FunctionContainer(const FunctionSpace& V)
> >> >> >>   {
> >> >> >>     _f = Function(V);
> >> >> >>   };
> >> >> >>
> >> >> >>   const Function& get_function()
> >> >> >>   {
> >> >> >>     return _f;
> >> >> >>   };
> >> >> >> protected:
> >> >> >>   Function _f;
> >> >> >> };
> >> >> >>
> >> >> >>
> >> >> >> int main()
> >> >> >> {
> >> >> >>   UnitSquare mesh(2, 2);
> >> >> >>   PoissonFunctionSpace V(mesh);
> >> >> >>   MyFunction f(V);
> >> >> >>   Vector x;
> >> >> >>
> >> >> >>   message("Interpolating to another vector");
> >> >> >>   f.interpolate(x, f.function_space());
> >> >> >>   x.disp();
> >> >> >>
> >> >> >>   message("Interpolating to the function vector");
> >> >> >>   f.interpolate(f.vector(), f.function_space());
> >> >> >>   f.vector().disp();
> >> >> >>
> >> >> >>   message("Interpolating using initialising by an external
> >> >> >> function"); MyFunction f_(f);
> >> >> >>   f.interpolate(f_.vector(), f.function_space());
> >> >> >>   f.vector().disp();
> >> >> >>
> >> >> >>   message("Returning Function by reference");
> >> >> >>   FunctionContainer fc(V);
> >> >> >>   Function f2 = fc.get_function();
> >> >> >> }
> >> >> >>
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> 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)
> >> >> >
> >> >> > iEYEARECAAYFAkmcOUcACgkQTuwUCDsYZdE/tACghYR+pHvXwurxKi2rKdcAPrtr
> >> >> > XaEAnihNPT9ar+ZLx07ltK+uZM03Ntlc
> >> >> > =8wBa
> >> >> > -----END PGP SIGNATURE-----
> >> >> >
> >> >> > _______________________________________________
> >> >> > DOLFIN-dev mailing list
> >> >> > DOLFIN-dev@xxxxxxxxxx
> >> >> > http://www.fenics.org/mailman/listinfo/dolfin-dev