[petsc-users] use of VecGetArray() with dof>1 (block size > 1)

[Dave May]

On 27 May 2016 at 20:34, Ed Bueler <elbueler at alaska.edu> wrote:

> Dear PETSc --
>
> This is an "am I using it correctly" question.  Probably the API has the
> current design because of something I am missing.
>
> First, a quote from the PETSc manual which I fully understand; it works
> great and gives literate code (to the extent possible...):
>
> """
> The recommended approach for multi-component PDEs is to declare a struct
> representing the fields defined
> at each node of the grid, e.g.
>
> typedef struct {
> PetscScalar u,v,omega,temperature;
> } Node;
>
> and write residual evaluation using
>
> Node **f,**u;
> DMDAVecGetArray(DM da,Vec local,&u);
> DMDAVecGetArray(DM da,Vec global,&f);
> ...
> f[i][j].omega = ...
>

Note that here the indexing should be
  f[ *j* ][ *i* ].omega



> ...
> DMDAVecRestoreArray(DM da,Vec local,&u);
> DMDAVecRestoreArray(DM da,Vec global,&f);
> """
>
> Now the three questions:
>
> 1)  The third argument to DMDAVec{Get,Restore}Array() is of type "void
> *".  It makes the above convenient.  But the third argument of the
> unstructured version Vec{Get,Restore}Array() is of type "PetscScalar **",
> which means that in an unstructured case, with the same Node struct, I
> would write
> "VecGetArray(DM da,Vec local,(PetscScalar **)&u);"
> to get the same functionality.  Why is it this way?  More specifically,
> why not have the argument to VecGetArray() be of type "void *"?
>

I would say the reason why the last arg to VecGetArray() is not void* is
because it is intended to give you direct access to the pointer associated
with the entries within the vector - these are also PetscScalar's


>
> 2) Given that the "recommended approach"
>

I don't believe it is ever recommended anywhere to do the following:
  VecGetArray(DM da,Vec local,(PetscScalar**)&u)

Trying to trick the compile with such a cast is just begging for memory
corruption to occur.

above works just fine, why do DMDAVec{Get,Restore}ArrayDOF() exist?  (I.e.
> is there something I am missing about C indexing?)
>

As an additional point, DMDAVec{Get,Restore}ArrayDOF() return

void *array

so that the same API will work for 1D, 2D and 3D DMDA's which would require
PetscScalar **data, PetscScalar ***data, PetscScalar ****data respectively.


Cheers,
  Dave


> 3) There are parts of the PETSc API that refer to "dof" and parts that
> refer to "block size".  Is this a systematic distinction with an underlying
> reason?  It seems "block size" is more generic, but also it seems that it
> could replace "dof" everywhere.
>
> Thanks for addressing silly questions.
>
> Ed
>
>
> --
> Ed Bueler
> Dept of Math and Stat and Geophysical Institute
> University of Alaska Fairbanks
> Fairbanks, AK 99775-6660
> 301C Chapman and 410D Elvey
> 907 474-7693 and 907 474-7199  (fax 907 474-5394)
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160527/04c50e27/attachment.html>