dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Thursday 19 February 2009 20:39:39 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-r
> >> >> >> >> >>>>efe rence/ d6/ da4/c
> >> >> >> >> >>>> lassdolfin_1_1GenericVector.html#bc4ecede55f2846ecadae0242
> >> >> >> >> >>>>6df 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:

vector() calls init() which asserts that a FunctionSpace is attached.

Johan

> 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