dolfin team mailing list archive

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

2007/12/12, Anders Logg <logg@xxxxxxxxx>:
> On Wed, Dec 12, 2007 at 11:45:27AM +0100, Martin Sandve Alnæs wrote:
> > 2007/12/12, Anders Logg <logg@xxxxxxxxx>:
> > > On Wed, Dec 12, 2007 at 09:03:20AM +0100, Martin Sandve Alnæs wrote:
> > > > 2007/12/11, Anders Logg <logg@xxxxxxxxx>:
> > > > > On Tue, Dec 11, 2007 at 02:58:58PM +0000, Garth N. Wells wrote:
> > > > > >
> > > > > >
> > > > > > Ola Skavhaug wrote:
> > > > > > > Hi,
> > > > > > >
> > > > > > > I'm trying to assemble dolfin/FFC forms into a user given matrix format
> > > > > > > CRSNew, which inherits GenericTensor as it should. However, the assemble.py in
> > > > > > > src/python in dolfin does not allow this.
> > > > > > >
> > > > > >
> > > > > > assembly.py has rather limited functionality compared to
> > > > > > dolfin::assemble. Does assembly into CRSNew using c++ work?.
> > > > >
> > > > > Yes, assemble() in PyDOLFIN currently assumes that you give it an FFC
> > > > > form and a Mesh. It then calls the JIT compiler on the form and then
> > > > > calls the wrapped C++ assembler on the generated UFC code.
> > > > >
> > > > > Does it work if you use cpp_assemble() instead of assemble()?
> > > > >
> > > > > We should try and extend assemble() in Python so it can recognize
> > > > > various combinations of input arguments (an FFC form or a UFC etc).
> > > > >
> > > > > The assemble() function is just a few lines of code and it's defined
> > > > > in assemble.py in src/pydolfin/.
> > > > >
> > > >
> > > > Have you tried using the Assembler and not the function assemble(..)?
> > > >
> > > > Earlier this autumn this worked for me with a ufc form from syfi (but
> > > > of course not with the new CRSNew matrix):
> > > >
> > > >
> > > > from dolfin import *
> > > >
> > > > mesh = UnitSquare(3, 3)
> > > > asm  = Assembler(mesh)
> > > >
> > > > coefficients = ArrayFunctionPtr()
> > > > cell_domains           = None
> > > > exterior_facet_domains = None
> > > > interior_facet_domains = None
> > > > reset_tensor = True
> > > >
> > > > A = Matrix() #A = CRSNew()
> > > > print type(A)
> > > > print isinstance(A, dolfin.GenericTensor)
> > > > asm.assemble(A, form, coefficients, cell_domains,
> > > > exterior_facet_domains, interior_facet_domains, reset_tensor)
> > > > A.disp()
> > > >
> > > >
> > > > Note that the default arguments for domains etc didn't work then, this
> > > > should be possible to achieve. Also, a typemap to handle a python
> > > > sequence of coefficients instead of the ArrayFunctionPtr would make it
> > > > much better. I have some typemaps in my pycc Assembler which checks
> > > > all items in a list for type and allows f.ex.
> > > >
> > > > A = asm.assemble_matrix(form, [1.23, (1.0, 2.0, 3.0), my_function,
> > > > my_numpy_array])
> > > >
> > > > It wraps constants and constant tuples into my equivalent of Function
> > > > (for those who know pycc, not the regular pycc::ScalarFunction etc).
> > > > For a GenericVector (again my equivalent of the dolfin one) it just
> > > > checks the length and assumes everything is fine.
> > > >
> > > > Anders, I can show you my typemap code, I think you can reuse much of
> > > > it for this.
> > >
> > > Yes, that would be nice. But I haven't checked Ola's code yet, maybe
> > > he has already fixed this? Sending Python lists instead of
> > > ArrayFunctionPtr (which should not really be exposed in the Python
> > > interface) would be very nice.
> >
> > If he has fixed it, it is probably only in assemble(...) and not in
> > the Assembler. But it might be a better idea to do the type
> > conversions in a python layer on top of Assembler than in typemaps
> > anyway. The advantage of typemaps is that they distribute the
> > list->ArrayFunctionPtr mapping to all C++ signatures, and we don't
> > need to add a python layer, but on the other side they are harder to
> > understand and maintain.
>
> I think the best is to add a small Python layer on top, something like
> what we have now in the assemble() function. It's easier than typemaps
> (which I don't understand very well) and it's easy to do special
> things like calling the JIT compiler before assembling, looking at
> arguments etc.

It's easier than typemaps, but the backside is that typemaps are
automatically applied to all C++ functions that take a
ArrayFunctionPtr argument, while this must be applied manually for
each case in the python layer. As long as it's a limited number of
cases, I guess it's ok. If it needs to scale with a lot of functions,
it means some additional maintenance.

-- 
Martin