dolfin team mailing list archive

[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-reference/
> >> >> >>>>d6/ da4/c
> >> >> >>>> lassdolfin_1_1GenericVector.html#bc4ecede55f2846ecadae02426df8f63)
> >> >> >>>> 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.

-- 
Anders


> 
> -Ali
> 
> >
> > Johan
> >
> >>
> >> -Ali
> >>
> >> > Johan
> >> >
> >> >> -Ali
> >> >>
> >> >> > Garth
> >> >> >
> >> >> >>>> Is there a better way of assigning the whole data at once?
> >> >> >>>
> >> >> >>> Maybe but you then have to be sure that your data is arrange in
> >> >> >>> exact same
> >> >> >>> manner as the mesh.
> >> >> >>
> >> >> >> As long as we work with complete rectangular images, which is the
> >> >> >> frequently used case, it should be fine. Given this, how is it
> >> >> >> possible to directly assign the image data pointer to _vector?
> >> >> >>
> >> >> >>>>> 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?
> >> >> >>>
> >> >> >>> If your raw image data is represented by a linear FunctionSpace any
> >> >> >>> interpolation using this space would be linear.
> >> >> >>
> >> >> >> Right, that's yet another problem.
> >> >> >>
> >> >> >>
> >> >> >> -Ali
> >> >> >>
> >> >> >>> Johan
> >> >> >>>
> >> >> >>>> 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.
> >> >> >>>>>>>>
> >> >> >>>>>>>>
> >> >> >>>>>>>> -----BEGIN PGP SIGNATURE-----
> >> >> >>>>>>>> Version: GnuPG v1.4.9 (GNU/Linux)
> >> >> >>>>>>>>
> >> >> >>>>>>>> iEYEARECAAYFAkmdMN8ACgkQTuwUCDsYZdGnOQCdFBJSD4FymLnVPbheRt63aJJ
> >> >> >>>>>>>>a 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
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> DOLFIN-dev mailing list
> >> >> >> DOLFIN-dev@xxxxxxxxxx
> >> >> >> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >
> >
> >

> // Place for random tests
> 
> #include <dolfin.h>
> #include "Poisson.h"
> 
> using namespace dolfin;
> 
> class Image
> {
> public:
> 
>   Image()
>   {
>     _size[0] = 10;
>     _size[1] = 10;
>     
>     for(int i=0; i<_size[0]; i++)
>       for(int j=0; j<_size[1]; j++)
>       {
>         int ij = i + _size[0] * j;
>         _data[ij] = 500.0*exp(-(i*i + j*j) / 0.02);
>       }
>   };
>   
>   int *get_size()
>   {
>     return _size;
>   }
>   
>   double *get_data()
>   {
>     return _data;
>   }
> 
> protected:
>   int _size[2];
>   double _data[100];
> };
> 
> class LinearDiscreteImageFunction : public Function
> {
> public:
>   LinearDiscreteImageFunction(boost::shared_ptr<Image> image) : Function()
>   {
>     if(!has_vector())
>     {
>       int N = image->get_size()[0] * image->get_size()[1];
>       
>       DefaultFactory factory;
>       message("[debug -1]");
>       vector() = *(boost::shared_ptr<GenericVector>(factory.create_vector()));
>       message("[debug  0]");
>       vector().resize(N);
>       message("[debug  1]");
>       vector().set(image->get_data());
>       message("[debug  2]");
>     }
>     else
>     {
>       message("This should never happen.");
>     }
>   };
> };
> 
> 
> int main()
> {  
>   message("Creating linear discrete image function");
>   boost::shared_ptr<Image> image(new Image());
>   LinearDiscreteImageFunction v(image);
> }
> 

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature