1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /*@ 4 DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM 5 6 Input Parameters: 7 + dm - The DM 8 - naturalSF - The PetscSF 9 10 Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map 11 12 Level: intermediate 13 14 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexGetMigrationSF()` 15 @*/ 16 PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF) 17 { 18 PetscFunctionBegin; 19 dm->sfMigration = migrationSF; 20 PetscCall(PetscObjectReference((PetscObject) migrationSF)); 21 PetscFunctionReturn(0); 22 } 23 24 /*@ 25 DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM 26 27 Input Parameter: 28 . dm - The DM 29 30 Output Parameter: 31 . migrationSF - The PetscSF 32 33 Level: intermediate 34 35 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateMigrationSF()`, `DMPlexSetMigrationSF` 36 @*/ 37 PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF) 38 { 39 PetscFunctionBegin; 40 *migrationSF = dm->sfMigration; 41 PetscFunctionReturn(0); 42 } 43 44 /*@ 45 DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec 46 47 Input Parameters: 48 + dm - The DM 49 - naturalSF - The PetscSF 50 51 Level: intermediate 52 53 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexGetGlobaltoNaturalSF()` 54 @*/ 55 PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF) 56 { 57 PetscFunctionBegin; 58 dm->sfNatural = naturalSF; 59 PetscCall(PetscObjectReference((PetscObject) naturalSF)); 60 dm->useNatural = PETSC_TRUE; 61 PetscFunctionReturn(0); 62 } 63 64 /*@ 65 DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec 66 67 Input Parameter: 68 . dm - The DM 69 70 Output Parameter: 71 . naturalSF - The PetscSF 72 73 Level: intermediate 74 75 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexCreateGlobalToNaturalSF()`, `DMPlexSetGlobaltoNaturalSF` 76 @*/ 77 PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF) 78 { 79 PetscFunctionBegin; 80 *naturalSF = dm->sfNatural; 81 PetscFunctionReturn(0); 82 } 83 84 /*@ 85 DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec 86 87 Input Parameters: 88 + dm - The redistributed DM 89 . section - The local PetscSection describing the Vec before the mesh was distributed, or NULL if not available 90 - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed 91 92 Output Parameter: 93 . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering 94 95 Note: This is not typically called by the user. 96 97 Level: intermediate 98 99 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()` 100 @*/ 101 PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural) 102 { 103 MPI_Comm comm; 104 Vec tmpVec; 105 PetscSF sf, sfEmbed, sfField; 106 PetscSection gSection, sectionDist, gLocSection; 107 PetscInt *spoints, *remoteOffsets; 108 PetscInt ssize, pStart, pEnd, p, globalSize; 109 PetscBool destroyFlag = PETSC_FALSE, debug = PETSC_FALSE; 110 111 PetscFunctionBegin; 112 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 113 if (!sfMigration) { 114 /* If sfMigration is missing, sfNatural cannot be computed and is set to NULL */ 115 *sfNatural = NULL; 116 PetscFunctionReturn(0); 117 } else if (!section) { 118 /* If the sequential section is not provided (NULL), it is reconstructed from the parallel section */ 119 PetscSF sfMigrationInv; 120 PetscSection localSection; 121 122 PetscCall(DMGetLocalSection(dm, &localSection)); 123 PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv)); 124 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion)); 125 PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section)); 126 PetscCall(PetscSFDestroy(&sfMigrationInv)); 127 destroyFlag = PETSC_TRUE; 128 } 129 if (debug) PetscCall(PetscSFView(sfMigration, NULL)); 130 /* Create a new section from distributing the original section */ 131 PetscCall(PetscSectionCreate(comm, §ionDist)); 132 PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist)); 133 PetscCall(PetscObjectSetName((PetscObject) sectionDist, "Migrated Section")); 134 if (debug) PetscCall(PetscSectionView(sectionDist, NULL)); 135 PetscCall(DMSetLocalSection(dm, sectionDist)); 136 /* If a sequential section is provided but no dof is affected, sfNatural cannot be computed and is set to NULL */ 137 PetscCall(DMCreateGlobalVector(dm, &tmpVec)); 138 PetscCall(VecGetSize(tmpVec, &globalSize)); 139 PetscCall(DMRestoreGlobalVector(dm, &tmpVec)); 140 if (globalSize) { 141 const PetscInt *leaves; 142 PetscInt *sortleaves, *indices; 143 PetscInt Nl; 144 145 /* Get a pruned version of migration SF */ 146 PetscCall(DMGetGlobalSection(dm, &gSection)); 147 if (debug) PetscCall(PetscSectionView(gSection, NULL)); 148 PetscCall(PetscSFGetGraph(sfMigration, NULL, &Nl, &leaves, NULL)); 149 PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd)); 150 for (p = pStart, ssize = 0; p < pEnd; ++p) { 151 PetscInt dof, off; 152 153 PetscCall(PetscSectionGetDof(gSection, p, &dof)); 154 PetscCall(PetscSectionGetOffset(gSection, p, &off)); 155 if ((dof > 0) && (off >= 0)) ++ssize; 156 } 157 PetscCall(PetscMalloc3(ssize, &spoints, Nl, &sortleaves, Nl, &indices)); 158 for (p = 0; p < Nl; ++p) {sortleaves[p] = leaves ? leaves[p] : p; indices[p] = p;} 159 PetscCall(PetscSortIntWithArray(Nl, sortleaves, indices)); 160 for (p = pStart, ssize = 0; p < pEnd; ++p) { 161 PetscInt dof, off, loc; 162 163 PetscCall(PetscSectionGetDof(gSection, p, &dof)); 164 PetscCall(PetscSectionGetOffset(gSection, p, &off)); 165 if ((dof > 0) && (off >= 0)) { 166 PetscCall(PetscFindInt(p, Nl, sortleaves, &loc)); 167 PetscCheck(loc >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %" PetscInt_FMT " with nonzero dof is not a leaf of the migration SF", p); 168 spoints[ssize++] = indices[loc]; 169 } 170 } 171 PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed)); 172 PetscCall(PetscObjectSetName((PetscObject) sfEmbed, "Embedded SF")); 173 PetscCall(PetscFree3(spoints, sortleaves, indices)); 174 if (debug) PetscCall(PetscSFView(sfEmbed, NULL)); 175 /* Create the SF associated with this section 176 Roots are natural dofs, leaves are global dofs */ 177 PetscCall(DMGetPointSF(dm, &sf)); 178 PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection)); 179 PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField)); 180 PetscCall(PetscSFDestroy(&sfEmbed)); 181 PetscCall(PetscSectionDestroy(&gLocSection)); 182 PetscCall(PetscObjectSetName((PetscObject) sfField, "Natural-to-Global SF")); 183 if (debug) PetscCall(PetscSFView(sfField, NULL)); 184 /* Invert the field SF 185 Roots are global dofs, leaves are natural dofs */ 186 PetscCall(PetscSFCreateInverseSF(sfField, sfNatural)); 187 PetscCall(PetscObjectSetName((PetscObject) *sfNatural, "Global-to-Natural SF")); 188 PetscCall(PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view")); 189 /* Clean up */ 190 PetscCall(PetscSFDestroy(&sfField)); 191 } else { 192 *sfNatural = NULL; 193 } 194 PetscCall(PetscSectionDestroy(§ionDist)); 195 PetscCall(PetscFree(remoteOffsets)); 196 if (destroyFlag) PetscCall(PetscSectionDestroy(§ion)); 197 PetscFunctionReturn(0); 198 } 199 200 /*@ 201 DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order. 202 203 Collective on dm 204 205 Input Parameters: 206 + dm - The distributed DMPlex 207 - gv - The global Vec 208 209 Output Parameters: 210 . nv - Vec in the canonical ordering distributed over all processors associated with gv 211 212 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 213 214 Level: intermediate 215 216 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()` 217 @*/ 218 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv) 219 { 220 const PetscScalar *inarray; 221 PetscScalar *outarray; 222 PetscMPIInt size; 223 224 PetscFunctionBegin; 225 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0)); 226 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 227 if (dm->sfNatural) { 228 PetscCall(VecGetArray(nv, &outarray)); 229 PetscCall(VecGetArrayRead(gv, &inarray)); 230 PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE)); 231 PetscCall(VecRestoreArrayRead(gv, &inarray)); 232 PetscCall(VecRestoreArray(nv, &outarray)); 233 } else if (size == 1) { 234 PetscCall(VecCopy(gv, nv)); 235 } else PetscCheck(!dm->useNatural,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov."); 236 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 237 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0)); 238 PetscFunctionReturn(0); 239 } 240 241 /*@ 242 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order. 243 244 Collective on dm 245 246 Input Parameters: 247 + dm - The distributed DMPlex 248 - gv - The global Vec 249 250 Output Parameter: 251 . nv - The natural Vec 252 253 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 254 255 Level: intermediate 256 257 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 258 @*/ 259 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv) 260 { 261 const PetscScalar *inarray; 262 PetscScalar *outarray; 263 PetscMPIInt size; 264 265 PetscFunctionBegin; 266 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 267 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 268 if (dm->sfNatural) { 269 PetscCall(VecGetArrayRead(gv, &inarray)); 270 PetscCall(VecGetArray(nv, &outarray)); 271 PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE)); 272 PetscCall(VecRestoreArrayRead(gv, &inarray)); 273 PetscCall(VecRestoreArray(nv, &outarray)); 274 } else if (size == 1) { 275 } else PetscCheck(!dm->useNatural,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov."); 276 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 277 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 278 PetscFunctionReturn(0); 279 } 280 281 /*@ 282 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order. 283 284 Collective on dm 285 286 Input Parameters: 287 + dm - The distributed DMPlex 288 - nv - The natural Vec 289 290 Output Parameters: 291 . gv - The global Vec 292 293 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 294 295 Level: intermediate 296 297 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()` 298 @*/ 299 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv) 300 { 301 const PetscScalar *inarray; 302 PetscScalar *outarray; 303 PetscMPIInt size; 304 305 PetscFunctionBegin; 306 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 307 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 308 if (dm->sfNatural) { 309 /* We only have access to the SF that goes from Global to Natural. 310 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM. 311 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */ 312 PetscCall(VecZeroEntries(gv)); 313 PetscCall(VecGetArray(gv, &outarray)); 314 PetscCall(VecGetArrayRead(nv, &inarray)); 315 PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 316 PetscCall(VecRestoreArrayRead(nv, &inarray)); 317 PetscCall(VecRestoreArray(gv, &outarray)); 318 } else if (size == 1) { 319 PetscCall(VecCopy(nv, gv)); 320 } else PetscCheck(!dm->useNatural,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov."); 321 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 322 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 323 PetscFunctionReturn(0); 324 } 325 326 /*@ 327 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order. 328 329 Collective on dm 330 331 Input Parameters: 332 + dm - The distributed DMPlex 333 - nv - The natural Vec 334 335 Output Parameters: 336 . gv - The global Vec 337 338 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 339 340 Level: intermediate 341 342 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 343 @*/ 344 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv) 345 { 346 const PetscScalar *inarray; 347 PetscScalar *outarray; 348 PetscMPIInt size; 349 350 PetscFunctionBegin; 351 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 352 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 353 if (dm->sfNatural) { 354 PetscCall(VecGetArrayRead(nv, &inarray)); 355 PetscCall(VecGetArray(gv, &outarray)); 356 PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 357 PetscCall(VecRestoreArrayRead(nv, &inarray)); 358 PetscCall(VecRestoreArray(gv, &outarray)); 359 } else if (size == 1) { 360 } else PetscCheck(!dm->useNatural,PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called and a section was defined, report to petsc-maint@mcs.anl.gov."); 361 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 362 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 363 PetscFunctionReturn(0); 364 } 365