[petsc-users] accessing individual elements of a global vector created from dmda

[Barry Smith]

  Bishesh,

    Do you need off process values? That is, might one process ask for a value that is not in the local part of the vector and belongs on a different process? If you know that the values accessed are always local then you can just use sol[k][j][i].vx   the k,j,i are the global indices over the entire mesh. 

    If any process might ask for any value on the entire mesh then you need to use DMDAGlobalToNaturalBegin/End() followed by VecScatterCreateToAll() and scatter all the values to all processes then you access them locally.  Note this is not scalable.


   Barry

Note that though with sol[k][j][i].vx it looks like you can access any value for k,j,i if you use a k, j, i that do not correspond to a value that belongs to this process you will get a memory error.




On Sep 17, 2013, at 11:06 AM, Bishesh Khanal <bisheshkh at gmail.com> wrote:

> Dear all,
> I have a vector say sol which is a vector created from a dmda with dof = 3.
> Now I need to provide an interface to the part of my code which does not use petsc but needs access to the values of sol. 
> 
> I was thinking of providing an interface such as:
> double solver::getSolutionAtPosition(int x, int y, int z, int field);
> which should return a value of the vector sol at the grid position (x,y,z) and dof = field.
> 
> I know how to get the array and loop through all the values such as:
> DMDALocalInfo info;
> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
> 
> Field ***sol;
> ierr = DMDAVecGetArray(da,b,&rhs);CHKERRQ(ierr);
> 
>     for (PetscInt k = info.zs; k<info.zs+info.zm; ++k) {
>         for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {
>             for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
>                 //access sol[k][j][i].vx   or sol[k][j][i].vy or sol[k][j][i].vz depending on field
> 
>             }
>         }
>     }
> 
> But if I have a particular (x,y,z) based on global co-ordinate system, how do I get the value for the given field ? 
> 
> Thanks,
> Bishesh