dolfin team mailing list archive

[A Navaei <axnavaei@xxxxxxxxxxxxxx>]

2009/2/19 Anders Logg <logg@xxxxxxxxx>:
> On Thu, Feb 19, 2009 at 08:37:14PM +0000, A Navaei wrote:
>> 2009/2/19 Anders Logg <logg@xxxxxxxxx>:
>> > On Thu, Feb 19, 2009 at 07:39:39PM +0000, A Navaei wrote:
>> >> 2009/2/19 Johan Hake <hake@xxxxxxxxx>:
>> >> > On Thursday 19 February 2009 17:51:12 A Navaei wrote:
>> >> >> 2009/2/19 Anders Logg <logg@xxxxxxxxx>:
>> >> >> > On Thu, Feb 19, 2009 at 02:53:52PM +0000, A Navaei wrote:
>> >> >> >> ---------- Forwarded message ----------
>> >> >> >> From: A Navaei <axnavaei@xxxxxxxxxxxxxx>
>> >> >> >> Date: 2009/2/19
>> >> >> >> Subject: Re: [DOLFIN-dev] A minimal c++ Function test and some bugs
>> >> >> >> To: Johan Hake <hake@xxxxxxxxx>
>> >> >> >>
>> >> >> >> 2009/2/19 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> > On Thursday 19 February 2009 14:54:57 A Navaei wrote:
>> >> >> >> >> 2009/2/19 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> >> > On Thursday 19 February 2009 14:46:18 A Navaei wrote:
>> >> >> >> >> >> 2009/2/19 Garth N. Wells <gnw20@xxxxxxxxx>:
>> >> >> >> >> >> > A Navaei wrote:
>> >> >> >> >> >> >> 2009/2/19 Johan Hake <hake@xxxxxxxxx>:
>> >> >> >> >> >> >>> The previous email was sent a bit premature...
>> >> >> >> >> >> >>>
>> >> >> >> >> >> >>> [snip]
>> >> >> >> >> >> >>>
>> >> >> >> >> >> >>>>> 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.
>> >> >> >> >> >> >>>
>> >> >> >> >> >> >>> Yes but you do not know how the mesh is ordered.
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> Consider a 2D image data mapped on a structured grid, with size
>> >> >> >> >> >> >> (number of nodes in each direction) of (sx, sy), the value of
>> >> >> >> >> >> >> image brightness at coordinate (i, j) is given by I(k), where:
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> k = i + sx * j
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> Is this not enough information to create the mesh?
>> >> >> >> >> >> >>
>> >> >> >> >> >> >>>> It seems _vector.set()
>> >> >> >> >> >> >>>>
>> >> >> >> >> >> >>>> (http://www.fenics.org/pub/documents/dolfin/dolfin-progr-refe
>> >> >> >> >> >> >>>>rence/ d6/ da4/c
>> >> >> >> >> >> >>>> lassdolfin_1_1GenericVector.html#bc4ecede55f2846ecadae02426df
>> >> >> >> >> >> >>>>8f63) accepts
>> >> >> >> >> >> >>>> the pointer to each row, is that right?
>> >> >> >> >> >> >>>
>> >> >> >> >> >> >>> There is no such thing as a row in a vector. It stores the
>> >> >> >> >> >> >>> data in a contigouos way (if not distributed).
>> >> >> >> >> >> >>
>> >> >> >> >> >> >> No idea why it says rows: virtual void
>> >> >> >> >> >> >> dolfin::GenericVector::set (const double *block, const uint
>> >> >> >> >> >> >> *num_rows, const uint *const * rows).
>> >> >> >> >> >> >
>> >> >> >> >> >> > You can set as few or as many entries as you wish at once.
>> >> >> >> >> >> > 'block' contains the values you want to set, 'num_rows' are the
>> >> >> >> >> >> > number of entries you wish to set, and 'rows' are the positions
>> >> >> >> >> >> > into which the entries should be inserted.
>> >> >> >> >> >>
>> >> >> >> >> >> I hate to say this, but _vector is private!
>> >> >> >> >> >
>> >> >> >> >> > ... and perfectly accesable through vector()
>> >> >> >> >>
>> >> >> >> >> ... and perfectly read-only through vector()
>> >> >> >> >
>> >> >> >> > I past from Function.h
>> >> >> >> >
>> >> >> >> >    /// Return the vector of expansion coefficients, automatically
>> >> >> >> >    /// initialized to zero if coefficients have not been computed
>> >> >> >> > (non-const version)
>> >> >> >> >    GenericVector& vector();
>> >> >> >> >
>> >> >> >> >    /// Return the vector of expansion coefficients, automatically
>> >> >> >> >    /// initialized to zero if coefficients have not been computed
>> >> >> >> > (const version)
>> >> >> >> >    const GenericVector& vector() const;
>> >> >> >> >
>> >> >> >> > So one const version and one that is not.
>> >> >> >>
>> >> >> >> Would this be the right way of initialising _vector in the subclass? :
>> >> >> >>
>> >> >> >>      DefaultFactory factory;
>> >> >> >>      vector() =
>> >> >> >> *(boost::shared_ptr<GenericVector>(factory.create_vector()));
>> >> >> >>
>> >> >> >> Try attached in your sandbox.
>> >> >> >
>> >> >> > No, just do
>> >> >> >
>> >> >> >  vector();
>> >> >> >
>> >> >> > That will initialize a zero vector of correct size.
>> >> >>
>> >> >> These two work:
>> >> >>
>> >> >>   Function1(boost::shared_ptr<const Image> image) : Function()
>> >> >>   {
>> >> >>   };
>> >> >>
>> >> >>   Function1() : Function()
>> >> >>   {
>> >> >>     vector();
>> >> >>   };
>> >> >>
>> >> >> But the combination generates a segmentation fault:
>> >> >>
>> >> >> class Function2 : public Function
>> >> >> {
>> >> >> public:
>> >> >>   Function2(boost::shared_ptr<const Image> image) : Function()
>> >> >>   {
>> >> >>     vector();
>> >> >>   };
>> >> >> };
>> >> >>
>> >> >> The trace includes /libdolfin.so.0(_ZN6dolfin8Function6vectorEv+0x38).
>> >> >> Example attached.
>> >> >
>> >> > I get:
>> >> >
>> >> > terminate called after throwing an instance of 'std::runtime_error'
>> >> >  what():  *** Assertion (_function_space) [at
>> >> > dolfin/function/Function.cpp:303 in init()]
>> >> >
>> >> > You need to attach a FunctionSpace before you call vector() otherwise it
>> >> > cannot know how long the vector should be.
>> >>
>> >> The vector length is supposed to be initialised directly using the
>> >> image size, not FunctionSpace:
>> >
>> > No, the correct way is the initialize the vector size from the
>> > dimension of the FunctionSpace.
>>
>> How do you fit a FunctionSpace in the below Function subclass? If we
>> want to create a FunctionSpace, it will at least need a mesh and we'll
>> be back to square one again.
>
> As I've said: create the Mesh and the FunctionSpace inside your
> Function constructor.
>
>> As we discussed, if we want to involve a FunctionSpace, the instance
>> should be of FFC-generated FunctionSpace subclass type. Choosing this
>> approach means that, as you mentioned, we have to generate many finite
>> element types which is basically repeating the creating of the already
>> existing pre-compiled finite elements in dolfin/elements/. Also you
>> mentioned that this pre-compiled approach is going to be obsolete
>> soon. So, why should I (1) repeat creating something which is going to
>> be (2) obsolete?
>
> As I've said: create your own FunctionSpace subclasses and instantiate
> the proper one inside your Function constructor. This will not be
> removed (but the stuff in dolfin/elements/ relying on signatures might
> be removed).

I don't think creating a database of pre-compiled FunctionSpaces is a
good idea. Here is a quote from you on this:

'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 do not intend to maintain and support a database of different
FunctionSpaces, the user should be able to create the one required
dynamicallt.

>
> Or if you want to be able to use an ImageFunction on any type of
> FunctionSpace, let the FunctionSpace be an argument to the
> constructor.

I definitely want the feature of ImageFunction to be created with any
type of available FunctionSpace.

This loop is never ending... . 'As I've said': FunctionSpace cannot be
an argument to the constructor of ImageFunction as FunctionSpace
constructor itself requires Mesh, but we want the Mesh to be created
inside ImageFunction constructor -- refer to Problem (2). In adittion,
you wrote below:

'Yes, and the FFC-generated subclasses is what you should use. No
alternatives are needed.'

... and the FFC-generated subclasses do require Mesh in their constructor.

... [FunctionSpace] --(requires in its ctor)--> (Mesh) --(requires to
be created in the ctor of)--> [ImageFunction] --(requires in its
ctor)--> [FunctionSpace] --...

I hope we will not go through this loop again.

>
>> I think the problem here is that dolfin does not still have a robust,
>> general and error-proof framework for creation and assigning different
>> types of finite elements.
>
> What we have now is probably neither robust nor error-proof
> but I don't think that is the problem here.

That's just criticising the library -- I hope you don't take it personally.


-Ali

>
>> Once this exists, creating the image
>> function will be possible and easy.
>
> It is already both possible and easy.
>
>> Just to summarise, here are the
>> problems for designing an image function class as a result of our long
>> discussions:
>>
>> (1) Only FFC-generated subclasses of FunctionSpace are useful in c++.
>> Direct subclassing FunctionSpace will require the use of signatures
>> which are error-prone and going to be obsolete anyway. Currently, no
>> other alternatives are available.
>
> Yes, and the FFC-generated subclasses is what you should use. No
> alternatives are needed.
>
>> (2) An image function requires the mesh to be created internally using
>> the image structure data, The only available recommended way of
>> creating a FunctionSpace subclass (FFC-generated) does needs a mesh in
>> its constructor.
>
> You should not create a FunctionSpace subclass.
>
>> (3) Following (2), we will end up with a data base of pre-compiled
>> finite element types, which is not desirable.
>
> You will end up with a set of precompiled FunctionSpaces, not
> FiniteElements.
>
>> So far over 160 emails have been communicated on this subject and we
>> didn't really get anywhere. I'd suggest to try to target the
>> fundamental problem (1) which would consequently solve other issues.
>> Otherwise, this loop of problems is simply a time waster.
>
> (1) is not a fundamental problem. If you try the hints I've suggested
> you might see that it actually works.
>
> --
> Anders
>
>
>>
>> -Ali
>>
>> >
>> >
>> >> class LinearDiscreteImageFunction : public Function
>> >> {
>> >> public:
>> >>   LinearDiscreteImageFunction(boost::shared_ptr<Image> image) : Function()
>> >>   {
>> >>   int N = image->get_size()[0] * image->get_size()[1];
>> >>   vector(); // segmentation fault
>> >>   vector().resize(N);
>> >>   vector().set(image->get_data());
>> >>   };
>> >> };
>> >>
>> >>
>> >> -Ali
>> >>
>> >> >
>> >> > The first constructor is not called (for what reason I do not know) and wont
>> >> > trigger the assertion.
>> >> >
>> >> > Johan
>> >> >
>> >> _______________________________________________
>> >> 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)
>> >
>> > iEYEARECAAYFAkmdui4ACgkQTuwUCDsYZdGXJwCePo99ya91OUbHu5rDMXrH/L7+
>> > yUEAnR8jdFM4LNwuR0YI8ZuYurxZMlek
>> > =NRhU
>> > -----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
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.9 (GNU/Linux)
>
> iEYEARECAAYFAkmdx1oACgkQTuwUCDsYZdFUfgCaA5BIPHpa7lP8Rx3VT2mPFtPf
> csUAoIeB2rkFy5WRNcnqP245sMsgp1th
> =hoOC
> -----END PGP SIGNATURE-----
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>
>