On Wed, 2010-12-15 at 13:26 -0600, Barry Smith wrote: > 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
Great, this is what I was missing. This should do the trick. Thanks Barry. Ethan > 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/ > > ------------------------------------- > > > -- ------------------------------------- Ethan Coon Post-Doctoral Researcher Mathematical Modeling and Analysis Los Alamos National Laboratory 505-665-8289 http://www.ldeo.columbia.edu/~ecoon/ -------------------------------------