dolfin team mailing list archive

[A Navaei <axnavaei@xxxxxxxxxxxxxx>]

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

>
> 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.

>
>>  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/attachment-0001.cpp


-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
>
>
>