static char help[] = "The main goal of this code is to retrieve the original element numbers as found in the " "initial partitions (sInitialPartition)... but after the call to DMPlexDistribute"; #include PetscReal sCoords2x5Mesh[18][2] = { {0.00000000000000000e+00, 0.00000000000000000e+00}, {2.00000000000000000e+00, 0.00000000000000000e+00}, {0.00000000000000000e+00, 1.00000000000000000e+00}, {2.00000000000000000e+00, 1.00000000000000000e+00}, {9.99999999997387978e-01, 0.00000000000000000e+00}, {9.99999999997387978e-01, 1.00000000000000000e+00}, {0.00000000000000000e+00, 2.00000000000000011e-01}, {0.00000000000000000e+00, 4.00000000000000022e-01}, {0.00000000000000000e+00, 5.99999999999999978e-01}, {0.00000000000000000e+00, 8.00000000000000044e-01}, {2.00000000000000000e+00, 2.00000000000000011e-01}, {2.00000000000000000e+00, 4.00000000000000022e-01}, {2.00000000000000000e+00, 5.99999999999999978e-01}, {2.00000000000000000e+00, 8.00000000000000044e-01}, {9.99999999997387756e-01, 2.00000000000000011e-01}, {9.99999999997387978e-01, 4.00000000000000022e-01}, {9.99999999997387978e-01, 6.00000000000000089e-01}, {9.99999999997388089e-01, 8.00000000000000044e-01}}; //Connectivity of a 2x5 rectangular mesh of quads : const PetscInt sConnectivity2x5Mesh[10][4] = { {0,4,14,6}, {6,14,15,7}, {7,15,16,8}, {8,16,17,9}, {9,17,5,2}, {4,1,10,14}, {14,10,11,15}, {15,11,12,16}, {16,12,13,17}, {17,13,3,5}}; const PetscInt sInitialPartition2x5Mesh[2][5] = { {0,2,4,6,8}, {1,3,5,7,9} }; const PetscInt sNLoclCells2x5Mesh = 5; const PetscInt sNGlobVerts2x5Mesh = 18; int main(int argc, char **argv) { const PetscInt Nc = sNLoclCells2x5Mesh; //Same on each rank for this example... const PetscInt Nv = sNGlobVerts2x5Mesh; const PetscInt* InitPartForRank[2] = {&sInitialPartition2x5Mesh[0][0], &sInitialPartition2x5Mesh[1][0]}; const PetscInt (*Conn)[4] = sConnectivity2x5Mesh; const PetscInt Ncor = 4; const PetscInt dim = 2; DM dm, idm, ddm; PetscSF sfVert, sfMig, sfPart; PetscPartitioner part; PetscSection s; PetscInt *cells, c; PetscMPIInt size, rank; PetscBool box = PETSC_FALSE, field = PETSC_FALSE; PetscErrorCode ierr; ierr = PetscInitialize(&argc, &argv, NULL, help);if (ierr) return ierr; ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRMPI(ierr); ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRMPI(ierr); if (size != 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a 2 processors example only"); ierr = PetscOptionsGetBool(NULL, NULL, "-box", &box, NULL);CHKERRQ(ierr); ierr = PetscOptionsGetBool(NULL, NULL, "-field", &field, NULL);CHKERRQ(ierr); ierr = DMPlexCreate(PETSC_COMM_WORLD, &dm);CHKERRQ(ierr); if (box) { ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr); ierr = DMSetFromOptions(dm);CHKERRQ(ierr); } else { ierr = PetscMalloc1(Nc * Ncor, &cells);CHKERRQ(ierr); for (c = 0; c < Nc; ++c) { PetscInt cell = (InitPartForRank[rank])[c], cor; for (cor = 0; cor < Ncor; ++cor) { cells[c*Ncor + cor] = Conn[cell][cor]; } } ierr = DMSetDimension(dm, dim);CHKERRQ(ierr); ierr = DMPlexBuildFromCellListParallel(dm, Nc, PETSC_DECIDE, Nv, Ncor, cells, &sfVert);CHKERRQ(ierr); //ierr = DMPlexBuildCoordinatesFromCellListParallel(dm, dim, sfVert, coords);CHKERRQ(ierr); ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr); ierr = PetscFree(cells);CHKERRQ(ierr); ierr = DMPlexInterpolate(dm, &idm);CHKERRQ(ierr); ierr = DMDestroy(&dm);CHKERRQ(ierr); dm = idm; } ierr = DMSetUseNatural(dm, PETSC_TRUE);CHKERRQ(ierr); ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); if (field) { const PetscInt Nf = 1; const PetscInt numComp[1] = {1}; const PetscInt numDof[3] = {0, 0, 1}; const PetscInt numBC = 0; ierr = DMSetNumFields(dm, Nf); CHKERRQ(ierr); ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &s); CHKERRQ(ierr); ierr = DMSetLocalSection(dm, s); CHKERRQ(ierr); ierr = PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); } ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); ierr = DMPlexDistribute(dm, 0, &sfMig, &ddm); CHKERRQ(ierr); ierr = PetscSFView(sfMig, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); ierr = PetscSFCreateInverseSF(sfMig, &sfPart); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) sfPart, "Inverse Migration SF");CHKERRQ(ierr); ierr = PetscSFView(sfPart, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); Vec lNatVec ; Vec lGlobalVec; ierr = DMGetGlobalVector(dm,&lNatVec); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) lNatVec, "Natural Vector (initial partition)");CHKERRQ(ierr); PetscScalar *lNatVecArray; ierr = VecGetArray(lNatVec, &lNatVecArray); CHKERRQ(ierr); //Copying the initial partition into the "natural" vector: for (c = 0; c < Nc; ++c) { lNatVecArray[c] = (InitPartForRank[rank])[c]; } ierr = VecRestoreArray(lNatVec, &lNatVecArray); CHKERRQ(ierr); ierr = DMGetGlobalVector(ddm,&lGlobalVec); CHKERRQ(ierr); ierr = PetscObjectSetName((PetscObject) lGlobalVec, "Global Vector (reordered element numbers in the petsc distributed order)");CHKERRQ(ierr); ierr = VecZeroEntries(lGlobalVec); CHKERRQ(ierr); // The call to DMPlexNaturalToGlobalBegin/End does not produce our expected result... // In lGlobalVec, we expect to have: /* * Process [0] * 2. * 4. * 8. * 3. * 9. * Process [1] * 1. * 5. * 7. * 0. * 6. */ ierr = DMPlexNaturalToGlobalBegin(ddm, lNatVec, lGlobalVec); CHKERRQ(ierr); ierr = DMPlexNaturalToGlobalEnd (ddm, lNatVec, lGlobalVec); CHKERRQ(ierr); ierr = VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); ierr = VecView(lGlobalVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); // Returning the Global vector to the natural one should give the initial values... ierr = VecZeroEntries(lNatVec); CHKERRQ(ierr); ierr = DMPlexGlobalToNaturalBegin(ddm, lGlobalVec, lNatVec); CHKERRQ(ierr); ierr = DMPlexGlobalToNaturalEnd (ddm, lGlobalVec, lNatVec); CHKERRQ(ierr); ierr = VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); //Verification ierr = VecGetArray(lNatVec, &lNatVecArray); CHKERRQ(ierr); //Verify that we got the initial partition back into the "natural" vector: PetscBool lResultBad = PETSC_FALSE; for (c = 0; c < Nc; ++c) { lResultBad |= lNatVecArray[c] != (InitPartForRank[rank])[c]; } ierr = VecRestoreArray(lNatVec, &lNatVecArray); CHKERRQ(ierr); // Now create an Overlap: //ierr = DMSetBasicAdjacency(ddm, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This works: it give only 2 cells for the overlap //ierr = DMSetAdjacency(ddm, 0, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This works too, but why do I have to use field #0 ??? Is it a convention? ierr = DMSetAdjacency(ddm, PETSC_DEFAULT, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr); //This does not works: it give 3 cells for the overlap instead of 2. DM ddmo; PetscSF sfo; ierr = DMPlexDistributeOverlap(ddm, 1, &sfo, &ddmo); CHKERRQ(ierr); IS lCellsWithOvl = 0; ierr = DMPlexGetCellNumbering(ddmo, &lCellsWithOvl); CHKERRQ(ierr); PetscScalar *lGlobVecArray = 0; ierr = VecGetArray(lGlobalVec, &lGlobVecArray); CHKERRQ(ierr); PetscInt n = 0; const PetscInt *lCellsPetscNum = 0; ierr = ISGetLocalSize(lCellsWithOvl, &n); CHKERRQ(ierr); ierr = ISGetIndices(lCellsWithOvl, &lCellsPetscNum); CHKERRQ(ierr); int i; for (i = 0; i < n; i++) { const int lGlobCellNum = lCellsPetscNum[i]; const int lLocalCellNum = lGlobCellNum - rank*5; //hard coded 5 elements by process if (lGlobCellNum >= 0) { const int lInitialCellNum = lGlobVecArray[lLocalCellNum]; printf("[%d] petsc_cell_# %d is initial (as in sInitialPartition) cell num %d\n", rank, lCellsPetscNum[i], lInitialCellNum); } else { printf("[%d] overlap petsc_cell_# %d \n", rank, -lCellsPetscNum[i]-1); } } ierr = ISRestoreIndices(lCellsWithOvl, &lCellsPetscNum); ierr = ISView(lCellsWithOvl,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr); ierr = VecRestoreArray(lGlobalVec, &lGlobVecArray); CHKERRQ(ierr); ierr = DMRestoreGlobalVector(dm,&lNatVec); CHKERRQ(ierr); ierr = DMRestoreGlobalVector(ddm,&lGlobalVec); CHKERRQ(ierr); ierr = PetscSFDestroy(&sfMig);CHKERRQ(ierr); ierr = PetscSFDestroy(&sfPart);CHKERRQ(ierr); ierr = DMDestroy(&dm);CHKERRQ(ierr); ierr = DMDestroy(&ddm);CHKERRQ(ierr); ierr = PetscFinalize(); return ierr || lResultBad; } /*TEST testset: args: -field nsize: 2 test: suffix: 0 args: test: suffix: 1 args: -box -dm_plex_simplex 0 -dm_plex_box_faces 7,10 -dm_distribute TEST*/