hi
I'm trying to use ISGlobalToLocalMappingApplyBlock() and am a bit
puzzled about the results it's giving.
I've attached a small test to illustrate. It just sets up a
local-to-global mapping with 10 elements. Running on two processes the
first has global indices 0 - 4 and the the second has 5 - 9. I then try
to find the local index corresponding to global index 8.
If I set the blocksize parameter to 1, it correctly gives the results -1
on rank 0 and 3 on rank 1.
But if I set the blocksize to 2 (or more), the results are -253701943 on
rank 0 and -1 on rank 1. Neither of these are what I expected- I thought
they should be the same as in the blocksize 1 case.
I'm presuming the global indices I pass in to
ISGlobalToLocalMappingApplyBlock() should be global block indices (i.e.
not scaled up by blocksize). If I do scale them up it doesn't give the
answers I expect either.
Or am I wrong to expect this to give the same results regardless of
blocksize?
Cheers, Adrian
--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: a.crouc...@auckland.ac.nz
tel: +64 (0)9 923 4611
program testl2g
! Test PETSc local-to-global mapping
#include <petsc/finclude/petsc.h>
use petsc
implicit none
PetscInt, parameter :: dof = 1, overlap = 1, n = 30
PetscInt :: dim, i, f, start_face, end_face, num_cells
PetscReal, parameter :: eps = 1.e-3
PetscReal :: area
PetscReal, pointer :: centroid(:), normal(:)
PetscMPIInt :: rank
DM :: dm, dist_dm, ghost_dm
PetscErrorCode :: ierr
character(40) :: filename = "col30.msh"
PetscInt, parameter :: bdy_face_index = 94
character(len = 16) :: boundary_label_name = "boundary"
PetscDS :: ds
PetscFV :: fvm
PetscInt :: global_indices(0: n-1), local_indices(0: n-1), nout
ISLocalToGlobalMapping :: l2g
call PetscInitialize(PETSC_NULL_CHARACTER, ierr); CHKERRQ(ierr)
call MPI_comm_rank(PETSC_COMM_WORLD, rank, ierr)
call DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, &
PETSC_TRUE, dm, ierr); CHKERRQ(ierr)
call DMPlexSetAdjacencyUseCone(dm, PETSC_TRUE, ierr); CHKERRQ(ierr)
call DMPlexSetAdjacencyUseClosure(dm, PETSC_FALSE, ierr); CHKERRQ(ierr)
call DMGetDimension(dm, dim, ierr); CHKERRQ(ierr)
call DMPlexDistribute(dm, overlap, PETSC_NULL_SF, dist_dm, ierr)
CHKERRQ(ierr)
if (dist_dm .ne. PETSC_NULL_DM) then
call DMDestroy(dm, ierr); CHKERRQ(ierr)
dm = dist_dm
end if
call DMCreateLabel(dm, boundary_label_name, ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(dm, 1, start_face, end_face, ierr); CHKERRQ(ierr)
allocate(centroid(3), normal(3))
do f = start_face, end_face - 1
call DMPlexGetSupportSize(dm, f, num_cells, ierr); CHKERRQ(ierr)
if (num_cells == 1) then
call DMPlexComputeCellGeometryFVM(dm, f, area, centroid, normal, ierr)
if (normal(dim) > eps) then
call DMSetLabelValue(dm, boundary_label_name, f, 1, ierr)
end if
end if
end do
deallocate(centroid, normal)
call DMPlexConstructGhostCells(dm, boundary_label_name, &
PETSC_NULL_INTEGER, ghost_dm, ierr); CHKERRQ(ierr)
if (ghost_dm .ne. PETSC_NULL_DM) then
call DMDestroy(dm, ierr); CHKERRQ(ierr);
dm = ghost_dm
end if
call DMView(dm, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)
call PetscFVCreate(PETSC_COMM_WORLD, fvm, ierr); CHKERRQ(ierr)
call PetscFVSetNumComponents(fvm, dof, ierr); CHKERRQ(ierr)
call PetscFVSetSpatialDimension(fvm, dim, ierr); CHKERRQ(ierr)
call DMGetDS(dm, ds, ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(ds, fvm, ierr); CHKERRQ(ierr)
call DMGetLocalToGlobalMapping(dm, l2g, ierr); CHKERRQ(ierr)
call ISLocalToGlobalMappingView(l2g, PETSC_VIEWER_STDOUT_WORLD, ierr); CHKERRQ(ierr)
global_indices = [(i, i = 0, n-1)]
call ISGlobalToLocalMappingApplyBlock(l2g, IS_GTOLM_MASK, n, global_indices, &
nout, local_indices, ierr); CHKERRQ(ierr)
write(*,*) 'rank:', rank, '(global, local):', &
[(global_indices(i), local_indices(i), i = 0, n - 1)]
call DMDestroy(dm, ierr); CHKERRQ(ierr)
call PetscFinalize(ierr); CHKERRQ(ierr)
end program testl2g
