# [SciPy-dev] FE+sparse utilities

Robert Cimrman cimrman3 at ntc.zcu.cz
Thu Dec 14 03:58:44 CST 2006

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

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