dolfin team mailing list archive

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

2007/11/1, Anders Logg <logg@xxxxxxxxx>:
> On Thu, Nov 01, 2007 at 07:50:15AM +0100, Martin Sandve Alnæs wrote:
> > 2007/11/1, Anders Logg <logg@xxxxxxxxx>:
> > > On Wed, Oct 31, 2007 at 10:10:42AM +0100, Martin Sandve Alnæs wrote:
> > > > 2007/10/31, Martin Sandve Alnæs <martinal@xxxxxxxxx>:
> > > > > 2007/10/30, Ola Skavhaug <skavhaug@xxxxxxxxx>:
> > > > > > Anders Logg skrev den 30/10-2007 følgende:
> > > > > > > The following code gives me a segmentation fault:
> > > > > > >
> > > > > > >   from dolfin import *
> > > > > > >
> > > > > > >   def foo():
> > > > > > >       element = FiniteElement("Lagrange", "triangle", 1)
> > > > > > >       mesh = UnitSquare(2, 2)
> > > > > > >       return Function(element, mesh, 1.0)
> > > > > > >
> > > > > > >   print foo().mesh().numVertices()
> > > > > > >
> > > > > > > If the mesh is declared outside foo() it works fine.
> > > > > > >
> > > > > > > It seems that Python is calling the Mesh destructor when foo()
> > > > > > > returns. It shouldn't since the Function created in foo() keeps a
> > > > > > > reference to the mesh.
> > > > > > >
> > > > > > > Does anyone know how to fix this? Is it possible to tell Python that
> > > > > > > the mesh is referenced by the function?
> > > > > >
> > > > > > Take a look at the swig documentation. It might be possible to add some
> > > > > > directives for this, but it is a major problem. Even the Python bindings for
> > > > > > vtk have this flaw.
> > > > >
> > > > > This is also the same as the memory handling issues we haven't given
> > > > > any serious thought in the GenericTensor hierarchy: A = Matrix();
> > > > > assemble(A); A = A.mat(); who owns the underlying PetSCMatrix or
> > > > > EpetraMatrix now?
> > > > >
> > > > > If you can accept a common base class with a reference count for all
> > > > > dolfin classes, and using str::tr1::shared_ptr in the C++ code, it can
> > > > > be solved (at least for most use cases).
> > > >
> > > > Sorry, that won't work for the linear algebra classes, because it
> > > > would require a reference count located in the Mat object (with
> > > > PetSC).
> > >
> > > The PETSc Mat pointer is created in the constructor of
> > > dolfin::PETScMatrix and destroyed in the destructor so we have full
> > > control over it.
> >
> > I think that's the best solution, since it's easy to understand, but
> > it doesn't solve the problem. Let me just explain what the problem is:
> >
> > # If the user does
> > def foo(...):
> >     A = PETScMatrix()
> >     assemble(A, ...)
> >     return A.mat()
> >
> > A = foo(...)
> > # ... do PETSc stuff only
> >
> > At this point the PETScMatrix is deleted, deleting the Mat, and A is
> > invalid. The solution is for the user to store the PETScMatrix
> > reference somewhere, but that kind of destroys the convenience of
> > reference counting in the first place... This is the same pattern as
> > your mesh example.
>
> This is not a big problem right now since no one (?) is using PETSc
> calls from Python.

Of course. It will become more relevant when using PyTrilinos, but
we'll just have to work around it when we get there. In PyTrilinos,
Epetra.Vector is a subclass of the C++ class Epetra_Vector, so we'll
have to do some trick anyway:

A = EpetraVector()
assemble(A, ...)
A2 = Epetra.Vector(A.vec(), Epetra.View) # or something like this

--
Martin