dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Feb 19, 2009 at 10:34:25AM +0000, A Navaei wrote:
> 2009/2/19 Anders Logg <logg@xxxxxxxxx>:
> > On Thu, Feb 19, 2009 at 09:36:01AM +0000, A Navaei wrote:
> >> 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.
> >
> > Yes, it should be avoided. It only leads to trouble. There's no reason
> > to use the FiniteElement class in C++. You should use the
> > FunctionSpace class.
> >
> >> >> 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.
> >
> > Yes, and that's why you should instantiate the FunctionSpace subclass
> > *inside* the ImageFunction constructor.
> 
> Is attached example what you mean by the above?

Yes, this looks good, I'd suggest making _image a reference and also
changing the constructor as follows:

 ImageFunction(const Image& image) : Function(), _image(image)
 {
   boost::shared_ptr<Mesh> mesh(new UnitSquare(image.get_size()[0] - 1, image.get_size()[1] - 1));
   boost::shared_ptr<FunctionSpace> V(new PoissonFunctionSpace(mesh));

   _function_space = V;
 };

> PoissonFunctionSpace
> simply stands for a general auto-generated function space subclass. In
> that case, what can be replaced with PoissonFunctionSpace which is
> general enough for solving an arbitrary PDE? If what you mean is
> generating many of these FunctionSpace subclasses, we will end up with
> something like dolfin/dolfin/elements/ffc_*.h which is not a good
> idea.

Yes, you will need one FunctionSpace subclass for each of the cases
you want to support. Then a switch in the Function constructor can
initialize the correct one depending on some flag.

> I simply want the image data to be directly incorporated in a PDE with
> the least effort and without losing generality.

For that to work, you will need to have a matching FunctionSpace in
the PDE definition. If that's important, it would be better to require
a FunctionSpace argument to the ImageFunction constructor.

-- 
Anders


> 
> -Ali
> 
> >
> >> >> 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:
> >
> > You don't need a mesh-free constructor for FunctionSpace, see above.
> >
> >> ' 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.'
> >
> > It's not needed if you do as I suggest.
> >
> >
> > -----BEGIN PGP SIGNATURE-----
> > Version: GnuPG v1.4.9 (GNU/Linux)
> >
> > iEYEARECAAYFAkmdMN8ACgkQTuwUCDsYZdGnOQCdFBJSD4FymLnVPbheRt63aJJa
> > yyoAn3KDuOmwd8ZX5YR1KucbafvieNBc
> > =lpyI
> > -----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 Image
> {
> public:
> 
>   Image()
>   {
>     _size[0] = 10;
>     _size[1] = 10;
>   };
>   
>   int *get_size()
>   {
>     return _size;
>   }
> 
> protected:
>   int _size[2];
> };
> 
> class ImageFunction : public Function
> {
> public:
>   ImageFunction() :Function() {};
>   ImageFunction(const Image& image) : Function(V)
>   {
>     _image = image;
>     
>     UnitSquare mesh(image.get_size()[0] - 1, image.get_size()[1] - 1);
>     _function_space = PoissonFunctionSpace(mesh);
>   };
>   
>   void eval(double* values, const double* x) const
>   {
>     message("[ImageFunction] Calling eval");
>     // set values from image data
>     // these are dummy
>     double dx = x[0] - 0.5;
>     double dy = x[1] - 0.5;
>     values[0] = 500.0*exp(-(dx*dx + dy*dy) / 0.02);
>   };
>   
> protected:
>   Image _image;
> };
> 
> 
> int main()
> {  
>   message("Creating imge function");  
>   Image image;
>   ImageFunction v(image);
>   v.interpolate();
> }
> 

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature