dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Thursday 19 February 2009 22:22:04 A Navaei wrote:
> 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-p
> >> >> >> >> >> >> >>>>rogr-refe rence/ d6/ da4/c
> >> >> >> >> >> >> >>>> lassdolfin_1_1GenericVector.html#bc4ecede55f2846ecad
> >> >> >> >> >> >> >>>>ae02426df 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.

If you want it dynamically and in c++ you need a precompiled 
Element/FunctionSpace library. Or you have to do the FunctionSpace part in 
PyDOLFIN.

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

How many different elements are of importance? Excuse my ignorance on image 
theory :) Isn't piecwise Lagrange the only element that is of importance when 
you _create_ your ImageFunction, as the pixels are defined at the vertices.

If this is the case you can hardcode this into the constructor and return to 
python. Here you can interpolate your image function to what ever 
FunctionSpace you want, created dynamically using the JIT compilation 
feature.

Johan

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