← Back to team overview

dolfin team mailing list archive

Re: A minimal c++ Function test and some bugs

 

2009/2/19 Anders Logg <logg@xxxxxxxxx>:
> On Thu, Feb 19, 2009 at 08:45:38AM +0000, A Navaei wrote:
>> 2009/2/19 Anders Logg <logg@xxxxxxxxx>:
>> > On Wed, Feb 18, 2009 at 07:56:59PM +0000, 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;
>> >>
>> >> 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
>> >
>> > I've already said that we can make the _function_space protected in
>> > the Function class. I've done so now and also added a new
>> > interpolation function. See if that helps.
>>
>> I should have missed that, I didn't notice the message from hg to
>> dolfin-dev about this (still cannot find it). I'm afraid, while this
>> would help, it's not still enough. Let's revise the design of an image
>> function again:
>>
>> My guess is that the following constructor for ImageFunction is
>> necessary and sufficient:
>>
>> Image image;
>> FiniteElement element;
>> ImageFunction v(image, element);
>
> Why not simply
>
>  Image image;
>  ImageFunction v(image);
>
> ?
>
> The class ImageFunction could take some optional arguments for how to
> represent the function, like degree, continuous or discontinuous etc.

Yes, ImageFunction can come with a default finite element and the user
can override. Again, this will need the use of FiniteElement and
DofMap internally, which you suggested to avoid.

>
>> The user should be able to specify the input image, and the type of
>> finite element of interest. The mesh completely depends on image size
>> and dimension and the user does not need to define it externally.
>> DofMap is also a function of Mesh, and involving it in the constructor
>> may be redundant. What do you think?
>
> I suggest avoiding FiniteElement and DofMap. Use generated subclasses of
> FunctionSpace instead.

Like I explained, the generated subclasses have a ctor which requires mesh:

PoissonBilinearFormFunctionSpace0(const dolfin::Mesh& mesh)

We want a mesh-free ctor for ImageFunction.

>
>> Another way of doing this could be by the use of an existing FunctionSpace:
>>
>>   UnitSquare dummy_mesh(1, 1);
>>   PoissonFunctionSpace V(dummy_mesh);
>>   ImageFunction v(image, V);
>>
>> Then, in the constructor of ImageFunction, V.element and V.dofmap can
>> be used to create another FunctionSpace which has a mesh created using
>> image:
>>
>>   ImageToFunction(Image image, const FunctionSpace& V)
>>   {
>>     // Create the function space
>>     UnitSquare mesh(image.get_size()[0] - 1, image.get_size()[1] - 1);
>>     FunctionSpace IV(mesh, V.element(), V.dofmap());
>>
>>     // ...
>>   };
>>
>> The problem with this approach is that it involves the use of a dummy mesh.
>>
>> A mesh-independent constructor added to FunctionSpace could help.
>> Alternatively, if a (protected) default (empty) constructor is added
>> to FunctionSpace, ImageFunctionSpace:FunctionSpace can have a
>> mesh-independent constructor. However, the FFC-generated function
>> spaces, eg PoissonFunctionSpace, still need a mesh.
>>
>> Hope this makes the problem more clear now.
>
> Create the mesh and the FunctionSpace subclass inside the
> ImageFunction constructor. Neither the mesh nor the function space
> need to be visible outside.

Again, there is no mesh-free ctor for FunctionSpace, and it doesn't
come with a default ctor so that the subclass can implement a
mesh-free ctor. I quote the above again:

' A mesh-independent constructor added to FunctionSpace could help.
 Alternatively, if a (protected) default (empty) constructor is added
 to FunctionSpace, ImageFunctionSpace:FunctionSpace can have a
 mesh-independent constructor.'


-Ali

>
> --
> Anders
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkmdI/IACgkQTuwUCDsYZdFDiwCfXL2JM2AkimAcuxut86eu9wTq
> +kEAn37EQK7B7dmehsUHvwiwjvhz6ojp
> =qP3/
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>
>


Follow ups

References