[petsc-users] IS from DA by coordinates
Ethan Coon
ecoon at lanl.gov
Wed Dec 15 15:04:29 CST 2010
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/
-------------------------------------
More information about the petsc-users
mailing list