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 DM 89 . section - The PetscSection describing the Vec before the mesh was distributed, 90 or NULL if not available 91 - sfMigration - The PetscSF used to distribute the mesh, or NULL if it cannot be computed 92 93 Output Parameter: 94 . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering 95 96 Note: This is not typically called by the user. 97 98 Level: intermediate 99 100 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()` 101 @*/ 102 PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural) 103 { 104 MPI_Comm comm; 105 Vec gv, tmpVec; 106 PetscSF sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv; 107 PetscSection gSection, sectionDist, gLocSection; 108 PetscInt *spoints, *remoteOffsets; 109 PetscInt ssize, pStart, pEnd, p, globalSize; 110 PetscLayout map; 111 PetscBool destroyFlag = PETSC_FALSE; 112 113 PetscFunctionBegin; 114 PetscCall(PetscObjectGetComm((PetscObject) dm, &comm)); 115 if (!sfMigration) { 116 /* If sfMigration is missing, 117 sfNatural cannot be computed and is set to NULL */ 118 *sfNatural = NULL; 119 PetscFunctionReturn(0); 120 } else if (!section) { 121 /* If the sequential section is not provided (NULL), 122 it is reconstructed from the parallel section */ 123 PetscSF sfMigrationInv; 124 PetscSection localSection; 125 126 PetscCall(DMGetLocalSection(dm, &localSection)); 127 PetscCall(PetscSFCreateInverseSF(sfMigration, &sfMigrationInv)); 128 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion)); 129 PetscCall(PetscSFDistributeSection(sfMigrationInv, localSection, NULL, section)); 130 PetscCall(PetscSFDestroy(&sfMigrationInv)); 131 destroyFlag = PETSC_TRUE; 132 } 133 /* PetscCall(PetscPrintf(comm, "Point migration SF\n")); 134 PetscCall(PetscSFView(sfMigration, 0)); */ 135 /* Create a new section from distributing the original section */ 136 PetscCall(PetscSectionCreate(comm, §ionDist)); 137 PetscCall(PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist)); 138 /* PetscCall(PetscPrintf(comm, "Distributed Section\n")); 139 PetscCall(PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD)); */ 140 PetscCall(DMSetLocalSection(dm, sectionDist)); 141 /* If a sequential section is provided but no dof is affected, 142 sfNatural cannot be computed and is set to NULL */ 143 PetscCall(DMCreateGlobalVector(dm, &tmpVec)); 144 PetscCall(VecGetSize(tmpVec, &globalSize)); 145 PetscCall(DMRestoreGlobalVector(dm, &tmpVec)); 146 if (globalSize) { 147 /* Get a pruned version of migration SF */ 148 PetscCall(DMGetGlobalSection(dm, &gSection)); 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(PetscMalloc1(ssize, &spoints)); 158 for (p = pStart, ssize = 0; p < pEnd; ++p) { 159 PetscInt dof, off; 160 161 PetscCall(PetscSectionGetDof(gSection, p, &dof)); 162 PetscCall(PetscSectionGetOffset(gSection, p, &off)); 163 if ((dof > 0) && (off >= 0)) spoints[ssize++] = p; 164 } 165 PetscCall(PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed)); 166 PetscCall(PetscFree(spoints)); 167 /* PetscCall(PetscPrintf(comm, "Embedded SF\n")); 168 PetscCall(PetscSFView(sfEmbed, 0)); */ 169 /* Create the SF for seq to natural */ 170 PetscCall(DMGetGlobalVector(dm, &gv)); 171 PetscCall(VecGetLayout(gv,&map)); 172 /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */ 173 PetscCall(PetscSFCreate(comm, &sfSeqToNatural)); 174 PetscCall(PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER)); 175 PetscCall(DMRestoreGlobalVector(dm, &gv)); 176 /* PetscCall(PetscPrintf(comm, "Seq-to-Natural SF\n")); 177 PetscCall(PetscSFView(sfSeqToNatural, 0)); */ 178 /* Create the SF associated with this section */ 179 PetscCall(DMGetPointSF(dm, &sf)); 180 PetscCall(PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection)); 181 PetscCall(PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField)); 182 PetscCall(PetscSFDestroy(&sfEmbed)); 183 PetscCall(PetscSectionDestroy(&gLocSection)); 184 /* PetscCall(PetscPrintf(comm, "Field SF\n")); 185 PetscCall(PetscSFView(sfField, 0)); */ 186 /* Invert the field SF so it's now from distributed to sequential */ 187 PetscCall(PetscSFCreateInverseSF(sfField, &sfFieldInv)); 188 PetscCall(PetscSFDestroy(&sfField)); 189 /* PetscCall(PetscPrintf(comm, "Inverse Field SF\n")); 190 PetscCall(PetscSFView(sfFieldInv, 0)); */ 191 /* Multiply the sfFieldInv with the */ 192 PetscCall(PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural)); 193 PetscCall(PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view")); 194 /* Clean up */ 195 PetscCall(PetscSFDestroy(&sfFieldInv)); 196 PetscCall(PetscSFDestroy(&sfSeqToNatural)); 197 } else { 198 *sfNatural = NULL; 199 } 200 PetscCall(PetscSectionDestroy(§ionDist)); 201 PetscCall(PetscFree(remoteOffsets)); 202 if (destroyFlag) PetscCall(PetscSectionDestroy(§ion)); 203 PetscFunctionReturn(0); 204 } 205 206 /*@ 207 DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order. 208 209 Collective on dm 210 211 Input Parameters: 212 + dm - The distributed DMPlex 213 - gv - The global Vec 214 215 Output Parameters: 216 . nv - Vec in the canonical ordering distributed over all processors associated with gv 217 218 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 219 220 Level: intermediate 221 222 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()` 223 @*/ 224 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv) 225 { 226 const PetscScalar *inarray; 227 PetscScalar *outarray; 228 PetscMPIInt size; 229 230 PetscFunctionBegin; 231 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0)); 232 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 233 if (dm->sfNatural) { 234 PetscCall(VecGetArray(nv, &outarray)); 235 PetscCall(VecGetArrayRead(gv, &inarray)); 236 PetscCall(PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE)); 237 PetscCall(VecRestoreArrayRead(gv, &inarray)); 238 PetscCall(VecRestoreArray(nv, &outarray)); 239 } else if (size == 1) { 240 PetscCall(VecCopy(gv, nv)); 241 } else { 242 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."); 243 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 244 } 245 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0)); 246 PetscFunctionReturn(0); 247 } 248 249 /*@ 250 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order. 251 252 Collective on dm 253 254 Input Parameters: 255 + dm - The distributed DMPlex 256 - gv - The global Vec 257 258 Output Parameter: 259 . nv - The natural Vec 260 261 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 262 263 Level: intermediate 264 265 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 266 @*/ 267 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv) 268 { 269 const PetscScalar *inarray; 270 PetscScalar *outarray; 271 PetscMPIInt size; 272 273 PetscFunctionBegin; 274 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 275 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 276 if (dm->sfNatural) { 277 PetscCall(VecGetArrayRead(gv, &inarray)); 278 PetscCall(VecGetArray(nv, &outarray)); 279 PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE)); 280 PetscCall(VecRestoreArrayRead(gv, &inarray)); 281 PetscCall(VecRestoreArray(nv, &outarray)); 282 } else if (size == 1) { 283 } else { 284 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."); 285 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 286 } 287 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 288 PetscFunctionReturn(0); 289 } 290 291 /*@ 292 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order. 293 294 Collective on dm 295 296 Input Parameters: 297 + dm - The distributed DMPlex 298 - nv - The natural Vec 299 300 Output Parameters: 301 . gv - The global Vec 302 303 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 304 305 Level: intermediate 306 307 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()` 308 @*/ 309 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv) 310 { 311 const PetscScalar *inarray; 312 PetscScalar *outarray; 313 PetscMPIInt size; 314 315 PetscFunctionBegin; 316 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 317 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 318 if (dm->sfNatural) { 319 /* We only have access to the SF that goes from Global to Natural. 320 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM. 321 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */ 322 PetscCall(VecZeroEntries(gv)); 323 PetscCall(VecGetArray(gv, &outarray)); 324 PetscCall(VecGetArrayRead(nv, &inarray)); 325 PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 326 PetscCall(VecRestoreArrayRead(nv, &inarray)); 327 PetscCall(VecRestoreArray(gv, &outarray)); 328 } else if (size == 1) { 329 PetscCall(VecCopy(nv, gv)); 330 } else { 331 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."); 332 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 333 } 334 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 335 PetscFunctionReturn(0); 336 } 337 338 /*@ 339 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order. 340 341 Collective on dm 342 343 Input Parameters: 344 + dm - The distributed DMPlex 345 - nv - The natural Vec 346 347 Output Parameters: 348 . gv - The global Vec 349 350 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 351 352 Level: intermediate 353 354 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 355 @*/ 356 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv) 357 { 358 const PetscScalar *inarray; 359 PetscScalar *outarray; 360 PetscMPIInt size; 361 362 PetscFunctionBegin; 363 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 364 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 365 if (dm->sfNatural) { 366 PetscCall(VecGetArrayRead(nv, &inarray)); 367 PetscCall(VecGetArray(gv, &outarray)); 368 PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 369 PetscCall(VecRestoreArrayRead(nv, &inarray)); 370 PetscCall(VecRestoreArray(gv, &outarray)); 371 } else if (size == 1) { 372 } else { 373 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."); 374 SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 375 } 376 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 377 PetscFunctionReturn(0); 378 } 379