[petsc-users] petsc4py DA with dof>1 confusion

[Lisandro Dalcin]

On 4 July 2011 16:19, Matthew Emmett <memmett at unc.edu> wrote:
> On Mon, Jul 4, 2011 at 2:52 PM, Lisandro Dalcin <dalcinl at gmail.com> wrote:
>> No, sorry... It has nothing to do with periodicity... It is actually a
>> C/Fortran ordering issue I need to fix...
>
> Ah, I see.  Thanks for looking into this so quickly.  Let me know if
> there is anything I can do on my end to help.  I will poke around the
> code for petsc4py, but you obviously know it better than I do.
>
>> In Fortran 90, it seems you index a DA Vec array as A[dof,x,y,z]... ,
>> However, I think that for Python we should follow a more C-ish
>> indexing A[x,y,z,dof]. Or we could do it like PETSc in C, that is
>> A[z,y,x,dof] (wich is the transpose of the Fortran way) but it is
>> counter-intuitive to C (and likely Python) programmers ...
>>
>> What do others think about this?
>
> I think A[x,y,z,dof] is probably the most intuitive for Python.
>

OK. I pushed a fix solving the issue. As a consequence, the underlying
NumPy array do have a mixed ordering (no C nor Fortran) where the
latest dim varies fastest, then first, second, third. So, you have to
be very careful about accessing the underlying array (I mean, doing
va.array[...] in your previous code), I would recomend you against
that, and code like this instead:


  with da.getVecArray(lvec) as va:
    if dof == 1:
      va[:,:] = rank+1
    else:
      va[:,:,0] = rank+1
  da.localToGlobal(lvec, gvec)

  da.globalToLocal(gvec, lvec)
  with da.getVecArray(lvec) as va:
    if dof == 1:
      print va[:,:]
    else:
      print va[:,:,0]




-- 
Lisandro Dalcin
---------------
CIMEC (INTEC/CONICET-UNL)
Predio CONICET-Santa Fe
Colectora RN 168 Km 472, Paraje El Pozo
3000 Santa Fe, Argentina
Tel: +54-342-4511594 (ext 1011)
Tel/Fax: +54-342-4511169