[petsc-users] IS from DA by coordinates

[Ethan Coon]

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