dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2007/7/27, Anders Logg <logg@xxxxxxxxx>:
> Martin Sandve Alnæs wrote:
> > 2007/7/26, Anders Logg <logg@xxxxxxxxx>:
> >> Martin Sandve Alnæs wrote:
> >> > 2007/7/25, Anders Logg <logg@xxxxxxxxx>:
> >> >> Martin Sandve Alnæs wrote:
> >> >> > 2007/7/25, Anders Logg <logg@xxxxxxxxx>:
> >> >> >> Martin Sandve Alnæs wrote:
> >> >> >> > 2007/7/25, Anders Logg <logg@xxxxxxxxx>:
> >> >> >> >> Garth N. Wells wrote:
> >> >> >> >>> Anders Logg wrote:
> >> >> >> >>>> I think there is a very simple solution. We just change the
> >> >> >> >>>> constructor of Function to
> >> >> >> >>>>
> >> >> >> >>>>     Function(Mesh& mesh, GenericVector& x, const Form&
> >> form, uint
> >> >> >> i = 1)
> >> >> >> >>>>
> >> >> >> >>>> and the vector member of DiscreteFunction to
> >> >> >> >>>>
> >> >> >> >>>>     GenericVector* x;
> >> >> >> >>>>
> >> >> >> >>>> Then if one needs a Function that uses a particular linear
> >> >> algebra
> >> >> >> >>>> backend, then just create the Function from a PETScVector,
> >> >> >> uBlasVector
> >> >> >> >>>> or EpetraVector. Otherwise, create it from a Vector (tied
> >> to the
> >> >> >> >>>> default backend) or let the Function itself decide (in which
> >> >> case it
> >> >> >> >>>> will be a Vector).
> >> >> >> >>>>
> >> >> >> >>> I thought that this would be a simple solution, but there a
> >> "new
> >> >> >> Vector"
> >> >> >> >>> in DiscreteFunction which requires that the vector type be
> >> >> known. If
> >> >> >> >>> this is somehow re-written, the above approach could be used.
> >> >> >> >>>
> >> >> >> >>> Garth
> >> >> >> >> This is used to create copies of the vector (when creating a sub
> >> >> >> >> function or in the copy constructor). How about adding a factory
> >> >> >> >> function to GenericVector? (And also to GenericFoo in general.)
> >> >> >> >>
> >> >> >> >>      virtual GenericVector* create() = 0;
> >> >> >> >>
> >> >> >> >> Then we can call this instead of new Vector() in
> >> >> DiscreteFunction.cpp.
> >> >> >> >>
> >> >> >> >> /Anders
> >> >> >> >
> >> >> >> > I already have this in pycc::GenericVector/Foo, but I call it
> >> >> "copy",
> >> >> >> > which I think is a better name.
> >> >> >>
> >> >> >> I also thought of naming it copy(), but this is not alwats
> >> >> appropriate.
> >> >> >> In DiscreteFunction, one needs to create a new vector which is
> >> not a
> >> >> >> copy (but which is a copy of a slice of the vector). Having a
> >> copy()
> >> >> >> function might be useful, but one can create a new object and then
> >> >> copy
> >> >> >>    the values by get/set.
> >> >> >
> >> >> > So you mean that A.create() returns a non-initialized matrix with
> >> >> > nothing in common with A? I find that a bit confusing, and see no
> >> >> > reason why A should have this ability.
> >> >>
> >> >> Yes, A.create() returns a new object that uses the same backend. Since
> >> >> the Function class can do things like creating a subfunction and it
> >> has
> >> >> a copy constructor, it needs to be able to look at a given Function
> >> and
> >> >> create a new Function which uses the same linear algebra backend.
> >> (There
> >> >> seems to be a need for this.
> >> >
> >> > Ok.
> >> >
> >> >> Otherwise, we could just let the Function
> >> >> class use the default backend.)
> >> >
> >> > No good.
> >> >
> >> >> > And values aren't the only thing that has to be copied, copy() would
> >> >> > copy or reuse the sparsity pattern and distribution pattern as well.
> >> >>
> >> >> We could have a special factory function for this and it could be
> >> named
> >> >> copy(). So I suggest one create() and one copy().
> >> >
> >> > Ok.
> >> >
> >> > We can do it like this:
> >> >
> >> > class GenericTensor
> >> > {
> >> >  GenericTensor * create() = 0;
> >> > }
> >> >
> >> > class GenericVector: public GenericTensor
> >> > {
> >> >  GenericVector * create() = 0;
> >> > }
> >> >
> >> > This C++-feature is perhaps not so widely known, but it's useful in
> >> > this context.
> >> >
> >> > Martin
> >>
> >> I haven't seen this before. Do you mean that the function returning a
> >> GenericVector will implement the function in GenericTensor returning
> >> GenericTensor?
> >
> > Yes. And if you have a GenericTensor* gt, gt->create() will return a
> > GenericTensor, while a GenericVector* gv; gv->create() returns a
> > GenericVector. It only works when using subclasses as return types, of
> > course. So if f.ex. using a smart pointer class, it wouldn't work. You
> > can even continue:
> >
> > class MyVec: public GenericVector
> > {
> >  MyVec * create();
> > }
> >
> > This way, we won't have to cast the pointers from create() and copy()
> > in code that uses f.ex. PetSCVector directly.
> >
> >> Do you have a reference/link explaining this feature?
> >
> > It's probably in Bjarne Stroustrups book.
> >
> > Martin
>
> ok, very good. And just to confirm, the important point seems to be not
> to write "virtual"?

No, that's an error, it's important to write "virtual" :-)
The important point is that the return types must be pointers
to subclasses of the return type used in the parent signature.

> I think this looks good. Any more comments/objections to adding create()
> and copy() to GenericFoo in la?
>
> Anything else we need to add to make this work for PyCC?

Nothing to do with PyCC, but there's the Epetra initialization.
Reusing Epetra_CrsGraphs for initializing new matrices are important.

-- 
---
Martin Sandve Alnæs
PhD student, Cardiac Computations
Simula Research Laboratory, Norway