[SciPy-dev] FE+sparse utilities

[Lisandro Dalcin]

On 12/13/06, Robert Cimrman <cimrman3 at ntc.zcu.cz> wrote:
> I have always wanted to have PETSc accessible in Python so that I could
> try it easily, great!

You can find it at PyPI. I routinely use it to solve NS equations with
monolitic FEM formulation using millons of elements (in a cluster, of
course). The element matrices computation and global matrix assemble
are both in the C side. But the Python side is great to manage the
linear/nonlinear (KSP/SNES) solution.

An let me say that PETSc is really a great, great library, specially
if you want to go to multiprocessor architectures. For sigle-proc
machines, you don't even need MPI.

> Feutils is a very simple package good for:
> 1. given finite element connectivity of degrees of freedom (= mesh
> connectivity, if 1 DOF per node), prepare a global CSR matrix with
> nonzero pattern exactly corresponding to that connectivity. Thus, when
> assembling, no reallocations, data shifts etc. are needed;
> 2. given a bunch of element contributions, assemble them into the global
> matrix.

So I understand you mean this bunch should be a contiguous n-dim
array. What (3, 1, 4, 4) stands for in he line below?

 mtxInEls = nm.ones( (3, 1, 4, 4), dtype = nm.float64 )

> The actual matrix allocation as well as element contribution assembling
> are done in C (via swig), too. But I use standard numpy (array) and
> scipy (CSR sparse matrix) data types to store the data. The key point is
> to compute in one C call many (at least 1000) local element
> contributions, and them assemble all of them also in one function call
> that way the time spent in Python is minimized.

Ok! So you work with scipy CSR matrices. I fully understand your
approach now. I think I will follow your approach in my own project
for a next release.

But I still think many users will feel more confortable of being able
to do something like

A[rows, cols] = values

where rows, cols, values should be interpreted as stacked arrays with
the right sizes. But then it there should be a way to specify if you
want to just put or instead add values in the global matrix structure.

Regards,



-- 
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594