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 //Connectivity of a 2x5 rectangular mesh of quads : 27 const PetscInt sConnectivity2x5Mesh[10][4] = { 28 {0,4,14,6}, 29 {6,14,15,7}, 30 {7,15,16,8}, 31 {8,16,17,9}, 32 {9,17,5,2}, 33 {4,1,10,14}, 34 {14,10,11,15}, 35 {15,11,12,16}, 36 {16,12,13,17}, 37 {17,13,3,5}}; 38 39 const PetscInt sInitialPartition2x5Mesh[2][5] = { 40 {0,2,4,6,8}, 41 {1,3,5,7,9} 42 }; 43 44 const PetscInt sNLoclCells2x5Mesh = 5; 45 const PetscInt sNGlobVerts2x5Mesh = 18; 46 47 int main(int argc, char **argv) 48 { 49 const PetscInt Nc = sNLoclCells2x5Mesh; //Same on each rank for this example... 50 const PetscInt Nv = sNGlobVerts2x5Mesh; 51 const PetscInt* InitPartForRank[2] = {&sInitialPartition2x5Mesh[0][0], 52 &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) { 83 cells[c*Ncor + cor] = Conn[cell][cor]; 84 } 85 } 86 PetscCall(DMSetDimension(dm, dim)); 87 PetscCall(DMPlexBuildFromCellListParallel(dm, Nc, PETSC_DECIDE, Nv, Ncor, cells, &sfVert, NULL)); 88 PetscCall(PetscSFDestroy(&sfVert)); 89 PetscCall(PetscFree(cells)); 90 PetscCall(DMPlexInterpolate(dm, &idm)); 91 PetscCall(DMDestroy(&dm)); 92 dm = idm; 93 } 94 PetscCall(DMSetUseNatural(dm, PETSC_TRUE)); 95 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 96 97 if (field) { 98 const PetscInt Nf = 1; 99 const PetscInt numComp[1] = {1}; 100 const PetscInt numDof[3] = {0, 0, 1}; 101 const PetscInt numBC = 0; 102 103 PetscCall(DMSetNumFields(dm, Nf)); 104 PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &s)); 105 PetscCall(DMSetLocalSection(dm, s)); 106 PetscCall(PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD)); 107 PetscCall(PetscSectionDestroy(&s)); 108 } 109 110 PetscCall(DMPlexGetPartitioner(dm, &part)); 111 PetscCall(PetscPartitionerSetFromOptions(part)); 112 113 PetscCall(DMPlexDistribute(dm, 0, &sfMig, &ddm)); 114 PetscCall(PetscSFView(sfMig, PETSC_VIEWER_STDOUT_WORLD)); 115 PetscCall(PetscSFCreateInverseSF(sfMig, &sfPart)); 116 PetscCall(PetscObjectSetName((PetscObject) sfPart, "Inverse Migration SF")); 117 PetscCall(PetscSFView(sfPart, PETSC_VIEWER_STDOUT_WORLD)); 118 119 Vec lGlobalVec, lNatVec; 120 PetscScalar *lNatVecArray; 121 122 { 123 PetscSection s; 124 125 PetscCall(DMGetGlobalSection(dm, &s)); 126 PetscCall(PetscSectionView(s, PETSC_VIEWER_STDOUT_WORLD)); 127 } 128 PetscCall(DMGetGlobalVector(dm, &lNatVec)); 129 PetscCall(PetscObjectSetName((PetscObject) lNatVec, "Natural Vector (initial partition)")); 130 131 //Copying the initial partition into the "natural" vector: 132 PetscCall(VecGetArray(lNatVec, &lNatVecArray)); 133 for (c = 0; c < Nc; ++c) lNatVecArray[c] = (InitPartForRank[rank])[c]; 134 PetscCall(VecRestoreArray(lNatVec, &lNatVecArray)); 135 136 PetscCall(DMGetGlobalVector(ddm,&lGlobalVec)); 137 PetscCall(PetscObjectSetName((PetscObject) lGlobalVec, "Global Vector (reordered element numbers in the petsc distributed order)")); 138 PetscCall(VecZeroEntries(lGlobalVec)); 139 140 // The call to DMPlexNaturalToGlobalBegin/End does not produce our expected result... 141 // In lGlobalVec, we expect to have: 142 /* 143 * Process [0] 144 * 2. 145 * 4. 146 * 8. 147 * 3. 148 * 9. 149 * Process [1] 150 * 1. 151 * 5. 152 * 7. 153 * 0. 154 * 6. 155 * 156 * but we obtained: 157 * 158 * Process [0] 159 * 2. 160 * 4. 161 * 8. 162 * 0. 163 * 0. 164 * Process [1] 165 * 0. 166 * 0. 167 * 0. 168 * 0. 169 * 0. 170 */ 171 172 { 173 PetscSF nsf; 174 175 PetscCall(DMPlexGetGlobalToNaturalSF(ddm, &nsf)); 176 PetscCall(PetscSFView(nsf, NULL)); 177 } 178 PetscCall(DMPlexNaturalToGlobalBegin(ddm, lNatVec, lGlobalVec)); 179 PetscCall(DMPlexNaturalToGlobalEnd (ddm, lNatVec, lGlobalVec)); 180 181 PetscCall(VecView(lNatVec, PETSC_VIEWER_STDOUT_WORLD)); 182 PetscCall(VecView(lGlobalVec, PETSC_VIEWER_STDOUT_WORLD)); 183 184 PetscCall(DMRestoreGlobalVector(dm,&lNatVec)); 185 PetscCall(DMRestoreGlobalVector(ddm,&lGlobalVec)); 186 187 PetscCall(PetscSFDestroy(&sfMig)); 188 PetscCall(PetscSFDestroy(&sfPart)); 189 PetscCall(DMDestroy(&dm)); 190 PetscCall(DMDestroy(&ddm)); 191 PetscCall(PetscFinalize()); 192 return 0; 193 } 194 195 /*TEST 196 197 testset: 198 args: -field -petscpartitioner_type simple 199 nsize: 2 200 201 test: 202 suffix: 0 203 args: 204 205 test: 206 suffix: 1 207 args: -box -dm_plex_simplex 0 -dm_plex_box_faces 2,5 -dm_distribute 208 209 TEST*/ 210