[SciPy-dev] FE+sparse utilities

[Robert Cimrman]

Lisandro Dalcin wrote:
> 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.

I do the same thing - logic in Python, real computations (not covered by 
numpy/scipy already) in C.

> 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.

Yeah, it is great. I have looked at it several times already in the 
past, but never really coded anything with it, though :(

>> 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 )

If you look at the next line
'assembleMatrix( mtx, mtxInEls, [0, 10, 50], conn )'
you can notice that the third argument is a list of three numbers, 
saying that mtxInEls contain element contribution of three elements 
conn[[0,10,50],:]. -> this is the '3' in the shape tuple.
next is the number of quadrature points, but you assemble local matrix 
already integrated over the element, so it equals '1'.
The last two numbers are the actual dimensions of the local matrix (4x4).

>> 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.

A[rows, cols] = values should put values, while
A[rows, cols] += values should add them?

I admit that my approach is a bit more low-level, but in the full code, 
this is not visible to the user, who just defines the terms of the 
equations and the rest is automatic (look at the simple example in 
poisson.pdf to get a feel).

And one could wrap it so that it is used under covers in expression like 
yours. :)

Cheers,
r.