dolfin team mailing list archive

[A Navaei <axnavaei@xxxxxxxxxxxxxx>]

2009/2/19 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 :)

Is it possible to set the pointer to _vector data directly where we
can pass the pointer to image data? If this works, no loop will be
involve in data conversion and it will be super fast.

It seems _vector.set()
(http://www.fenics.org/pub/documents/dolfin/dolfin-progr-reference/d6/da4/classdolfin_1_1GenericVector.html#bc4ecede55f2846ecadae02426df8f63)
accepts the pointer to each row, is that right? Is there a better way
of assigning the whole data at once?

>
> Note that this will only work for a piecewise linear lagrange FunctionSpace.

It should be ok, it would be equivalent of having the image function
initialised with Lagrange element, but much more faster.

> But haveing this in place you can always interpolate to other FunctionSpaces.
> However these interpolations will only be linear of course.

It is not possible to interpolate a function to another space using
non-linear interpolators?

To Anders: Do you think if this is more practical and efficient than
the other discussed approach? Would this be general enough to work
with an arbitrary PDE of any finite element type?


-Ali

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