On Thu, Feb 19, 2009 at 02:53:52PM +0000, A Navaei wrote: > ---------- Forwarded message ---------- > From: A Navaei <[email protected]> > Date: 2009/2/19 > Subject: Re: [DOLFIN-dev] A minimal c++ Function test and some bugs > To: Johan Hake <[email protected]> > > > 2009/2/19 Johan Hake <[email protected]>: > > On Thursday 19 February 2009 14:54:57 A Navaei wrote: > >> 2009/2/19 Johan Hake <[email protected]>: > >> > On Thursday 19 February 2009 14:46:18 A Navaei wrote: > >> >> 2009/2/19 Garth N. Wells <[email protected]>: > >> >> > A Navaei wrote: > >> >> >> 2009/2/19 Johan Hake <[email protected]>: > >> >> >>> 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
> >> >> >>>>>>>> [email protected]
> >> >> >>>>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >> >> >>>>>>
> >> >> >>>>>> _______________________________________________
> >> >> >>>>>> DOLFIN-dev mailing list
> >> >> >>>>>> [email protected]
> >> >> >>>>>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >> >> >>
> >> >> >> _______________________________________________
> >> >> >> DOLFIN-dev mailing list
> >> >> >> [email protected]
> >> >> >> 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
> [email protected]
> http://www.fenics.org/mailman/listinfo/dolfin-dev
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
