dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Thursday 19 February 2009 12:16:45 A Navaei wrote:
> 2009/2/19 A Navaei <axnavaei@xxxxxxxxxxxxxx>:
> > 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? 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.
> >
> > I simply want the image data to be directly incorporated in a PDE with
> > the least effort and without losing generality.
>
> One more thing, does it make sense to have the image function as an
> already-interpolated function by directly setting _vector using the
> image data? Discretising descrete image data may not be the best
> approach, it will also help us to get rid of all this mesh, element
> and dofmap stuff.

I have also thought of this :)

Then you can take an ordered mesh, iterate over the vertices, and use the 
present implementation of the eval function (which in its present shape is a 
piecewise constant interpolator isn't it?), to fill the new _vector, all done 
in the constructor. Then you have a discrete function in the first place :)

Note that this will only work for a piecewise linear lagrange FunctionSpace. 
But haveing this in place you can always interpolate to other FunctionSpaces. 
However these interpolations will only be linear of course.

Johan


>
>
> -Ali
>
> > -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.
> >>
> >> --
> >> Anders
> >>
> >> -----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
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev