1 static char help[] = "The main goal of this code is to retrieve the original element numbers as found in the " 2 "initial partitions (sInitialPartition)... but after the call to DMPlexDistribute"; 3 4 #include <petsc.h> 5 6 PetscReal sCoords2x5Mesh[18][2] = { 7 {0.00000000000000000e+00, 0.00000000000000000e+00}, 8 {2.00000000000000000e+00, 0.00000000000000000e+00}, 9 {0.00000000000000000e+00, 1.00000000000000000e+00}, 10 {2.00000000000000000e+00, 1.00000000000000000e+00}, 11 {9.99999999997387978e-01, 0.00000000000000000e+00}, 12 {9.99999999997387978e-01, 1.00000000000000000e+00}, 13 {0.00000000000000000e+00, 2.00000000000000011e-01}, 14 {0.00000000000000000e+00, 4.00000000000000022e-01}, 15 {0.00000000000000000e+00, 5.99999999999999978e-01}, 16 {0.00000000000000000e+00, 8.00000000000000044e-01}, 17 {2.00000000000000000e+00, 2.00000000000000011e-01}, 18 {2.00000000000000000e+00, 4.00000000000000022e-01}, 19 {2.00000000000000000e+00, 5.99999999999999978e-01}, 20 {2.00000000000000000e+00, 8.00000000000000044e-01}, 21 {9.99999999997387756e-01, 2.00000000000000011e-01}, 22 {9.99999999997387978e-01, 4.00000000000000022e-01}, 23 {9.99999999997387978e-01, 6.00000000000000089e-01}, 24 {9.99999999997388089e-01, 8.00000000000000044e-01} 25 }; 26 27 //Connectivity of a 2x5 rectangular mesh of quads : 28 const PetscInt sConnectivity2x5Mesh[10][4] = { 29 {0, 4, 14, 6 }, 30 {6, 14, 15, 7 }, 31 {7, 15, 16, 8 }, 32 {8, 16, 17, 9 }, 33 {9, 17, 5, 2 }, 34 {4, 1, 10, 14}, 35 {14, 10, 11, 15}, 36 {15, 11, 12, 16}, 37 {16, 12, 13, 17}, 38 {17, 13, 3, 5 } 39 }; 40 41 const PetscInt sInitialPartition2x5Mesh[2][5] = { 42 {0, 2, 4, 6, 8}, 43 {1, 3, 5, 7, 9} 44 }; 45 46 const PetscInt sNLoclCells2x5Mesh = 5; 47 const PetscInt sNGlobVerts2x5Mesh = 18; 48 49 int main(int argc, char **argv) { 50 const PetscInt Nc = sNLoclCells2x5Mesh; //Same on each rank for this example... 51 const PetscInt Nv = sNGlobVerts2x5Mesh; 52 const PetscInt *InitPartForRank[2] = {&sInitialPartition2x5Mesh[0][0], &sInitialPartition2x5Mesh[1][0]}; 53 const PetscInt(*Conn)[4] = sConnectivity2x5Mesh; 54 55 const PetscInt Ncor = 4; 56 const PetscInt dim = 2; 57 DM dm, idm, ddm; 58 PetscSF sfVert, sfMig, sfPart; 59 PetscPartitioner part; 60 PetscSection s; 61 PetscInt *cells, c; 62 PetscMPIInt size, rank; 63 PetscBool box = PETSC_FALSE, field = PETSC_FALSE; 64 65 PetscFunctionBeginUser; 66 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 67 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 68 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 69 PetscCheck(size == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "This is a 2 processors example only"); 70 PetscCall(PetscOptionsGetBool(NULL, NULL, "-box", &box, NULL)); 71 PetscCall(PetscOptionsGetBool(NULL, NULL, "-field", &field, NULL)); 72 73 PetscCall(DMPlexCreate(PETSC_COMM_WORLD, &dm)); 74 if (box) { 75 PetscCall(DMSetType(dm, DMPLEX)); 76 PetscCall(DMSetFromOptions(dm)); 77 } else { 78 PetscCall(PetscMalloc1(Nc * Ncor, &cells)); 79 for (c = 0; c < Nc; ++c) { 80 PetscInt cell = (InitPartForRank[rank])[c], cor; 81 82 for (cor = 0; cor < Ncor; ++cor) cells[c * Ncor + cor] = Conn[cell][cor]; 83 } 84 PetscCall(DMSetDimension(dm, dim)); 85 PetscCall(DMPlexBuildFromCellListParallel(dm, Nc, PETSC_DECIDE, Nv, Ncor, cells, &sfVert, NULL)); 86 PetscCall(PetscSFDestroy(&sfVert)); 87 PetscCall(PetscFree(cells)); 88 PetscCall(DMPlexInterpolate(dm, &idm)); 89 PetscCall(DMDestroy(&dm)); 90 dm = idm; 91 } 92 PetscCall(DMSetUseNatural(dm, PETSC_TRUE)); 93 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 94 95 if (field) { 96 const PetscInt Nf = 1; 97 const PetscInt numComp[1] = {1}; 98 const PetscInt numDof[3] = {0, 0, 1}; 99 const PetscInt numBC = 0; 100 101 PetscCall(DMSetNumFields(dm, Nf)); 102 PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &s)); 103 PetscCall(DMSetLocalSection(dm, s)); 104 PetscCall(PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD)); 105 PetscCall(PetscSectionDestroy(&s)); 106 } 107 108 PetscCall(DMPlexGetPartitioner(dm, &part)); 109 PetscCall(PetscPartitionerSetFromOptions(part)); 110 111 PetscCall(DMPlexDistribute(dm, 0, &sfMig, &ddm)); 112 PetscCall(PetscSFView(sfMig, PETSC_VIEWER_STDOUT_WORLD)); 113 PetscCall(PetscSFCreateInverseSF(sfMig, &sfPart)); 114 PetscCall(PetscObjectSetName((PetscObject)sfPart, "Inverse Migration SF")); 115 PetscCall(PetscSFView(sfPart, PETSC_VIEWER_STDOUT_WORLD)); 116 117 Vec lGlobalVec, lNatVec; 118 PetscScalar *lNatVecArray; 119 120 { 121 PetscSection s; 122 123 PetscCall(DMGetGlobalSection(dm, &s)); 124 PetscCall(PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD)); 125 } 126 PetscCall(DMGetGlobalVector(dm, &lNatVec)); 127 PetscCall(PetscObjectSetName((PetscObject)lNatVec, "Natural Vector (initial partition)")); 128 129 //Copying the initial partition into the "natural" vector: 130 PetscCall(VecGetArray(lNatVec, &lNatVecArray)); 131 for (c = 0; c < Nc; ++c) lNatVecArray[c] = (InitPartForRank[rank])[c]; 132 PetscCall(VecRestoreArray(lNatVec, &lNatVecArray)); 133 134 PetscCall(DMGetGlobalVector(ddm, &lGlobalVec)); 135 PetscCall(PetscObjectSetName((PetscObject)lGlobalVec, "Global Vector (reordered element numbers in the petsc distributed order)")); 136 PetscCall(VecZeroEntries(lGlobalVec)); 137 138 // The call to DMPlexNaturalToGlobalBegin/End does not produce our expected result... 139 // In lGlobalVec, we expect to have: 140 /* 141 * Process [0] 142 * 2. 143 * 4. 144 * 8. 145 * 3. 146 * 9. 147 * Process [1] 148 * 1. 149 * 5. 150 * 7. 151 * 0. 152 * 6. 153 * 154 * but we obtained: 155 * 156 * Process [0] 157 * 2. 158 * 4. 159 * 8. 160 * 0. 161 * 0. 162 * Process [1] 163 * 0. 164 * 0. 165 * 0. 166 * 0. 167 * 0. 168 */ 169 170 { 171 PetscSF nsf; 172 173 PetscCall(DMPlexGetGlobalToNaturalSF(ddm, &nsf)); 174 PetscCall(PetscSFView(nsf, NULL)); 175 } 176 PetscCall(DMPlexNaturalToGlobalBegin(ddm, lNatVec, lGlobalVec)); 177 PetscCall(DMPlexNaturalToGlobalEnd(ddm, lNatVec, lGlobalVec)); 178 179 PetscCall(VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD)); 180 PetscCall(VecView(lGlobalVec, PETSC_VIEWER_STDOUT_WORLD)); 181 182 PetscCall(DMRestoreGlobalVector(dm, &lNatVec)); 183 PetscCall(DMRestoreGlobalVector(ddm, &lGlobalVec)); 184 185 PetscCall(PetscSFDestroy(&sfMig)); 186 PetscCall(PetscSFDestroy(&sfPart)); 187 PetscCall(DMDestroy(&dm)); 188 PetscCall(DMDestroy(&ddm)); 189 PetscCall(PetscFinalize()); 190 return 0; 191 } 192 193 /*TEST 194 195 testset: 196 args: -field -petscpartitioner_type simple 197 nsize: 2 198 199 test: 200 suffix: 0 201 args: 202 203 test: 204 suffix: 1 205 args: -box -dm_plex_simplex 0 -dm_plex_box_faces 2,5 -dm_distribute 206 207 TEST*/ 208