[petsc-users] IS from DA by coordinates

[Barry Smith]

  Ethan,

  I don't think there is any reason to create a new DA or call DAGetGlobalIndices().. Just call DAGetLocalToGlobalMappingBlock() on the original DA  Then 

for i in range(xs, xs+xl):
 for j in range(ys, ys+yl):
   for k in range(zs, zs+zl):
     if vecc_a[k,j,i,:] == whatever_coordinate:
      local_indices.append( convert i,j,k, to local numbering with something like (k-zs)*mx*my + (j-ys)*mx + ..
...

Now apply ISLocalToGlobalMappingApply to local_indices and you have a list of global indices depending on what you want you do with this beast you may need to scale by bs or 1/bs

   Barry


On Dec 15, 2010, at 1:09 PM, Ethan Coon wrote:

> Hi all,
> 
> Is there a cleaner way to create an IS to a global vector on a DA for a
> subset of nodes using a coordinate value than the following?  (Pardon my
> pseudo-code combo of python and c and imprecise arguments)  
> 
> // get the coordinates of the nodes
> DAGetCoordinateDA(da, dac)
> DAGetCoordinates(da, vecc)
> DAVecGetArray(dac, vecc, vecc_a)
> 
> // generate a one-dof da with no ghosts and the same parallel 
> // structure as the da
> DAGetOwnershipRanges(da, lx[], ly[], lz[])
> DACreate3D(comm, M,N,P,len(lx),len(ly),len(lz),1,0, lx[], ly[], \
>  lz[], da_one)
> 
> // get the global indices of the one-dof da, noting that because
> // we set the stencil size to zero, the local array of global 
> // indices is the same size as the local portion of the global array 
> // of coordinates
> DAGetCorners(xs, ys, zs, xl, yl, zl)
> DAGetGlobalIndices(da_one, xl*yl*zl, indices[])
> 
> // loop over the local array returned by 
> // DAGetGlobalIndices and the coordinates, 
> // comparing to the test
> global_block_indices = []
> for i in range(xs, xs+xl):
>  for j in range(ys, ys+yl):
>    for k in range(zs, zs+zl):
>      if vecc_a[k,j,i,:] == whatever_coordinate:
>        global_block_indices.append(indices[k,j,i])
> 
> 
> // make the IS
> ISCreateBlock(comm, ndofs, global_block_indices, coord_is)
> 
> // restore and destroy etc.
> 
> This just seems quite complicated with the construction of the one-dof
> da to get global indices of the block.  There might be a better way with
> ISLocalToGlobalMapping, but I wasn't sure how.  Any suggestions?
> 
> Thanks,
> 
> Ethan
> 
> 
> 
> 
> -- 
> -------------------------------------
> Ethan Coon
> Post-Doctoral Researcher
> Mathematical Modeling and Analysis
> Los Alamos National Laboratory
> 505-665-8289
> 
> http://www.ldeo.columbia.edu/~ecoon/
> -------------------------------------
>