dolfin team mailing list archive

[A Navaei <axnavaei@xxxxxxxxxxxxxx>]

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


-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.
>> >> >>>>>>>>
>> >> >>>>>>>> --
>> >> >>>>>>>> Anders
>> >> >>>>>>>>
>> >> >>>>>>>> -----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);
}