[petsc-users] IS from DA by coordinates

[Ethan Coon]

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