[SciPy-dev] FE+sparse utilities

Lisandro Dalcin dalcinl at gmail.com
Wed Dec 13 13:02:19 CST 2006

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.


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

More information about the Scipy-dev mailing list