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 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."); 242 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 243 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0)); 244 PetscFunctionReturn(0); 245 } 246 247 /*@ 248 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order. 249 250 Collective on dm 251 252 Input Parameters: 253 + dm - The distributed DMPlex 254 - gv - The global Vec 255 256 Output Parameter: 257 . nv - The natural Vec 258 259 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 260 261 Level: intermediate 262 263 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 264 @*/ 265 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv) 266 { 267 const PetscScalar *inarray; 268 PetscScalar *outarray; 269 PetscMPIInt size; 270 271 PetscFunctionBegin; 272 PetscCall(PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 273 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 274 if (dm->sfNatural) { 275 PetscCall(VecGetArrayRead(gv, &inarray)); 276 PetscCall(VecGetArray(nv, &outarray)); 277 PetscCall(PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray,MPI_REPLACE)); 278 PetscCall(VecRestoreArrayRead(gv, &inarray)); 279 PetscCall(VecRestoreArray(nv, &outarray)); 280 } else if (size == 1) { 281 } 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."); 282 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 283 PetscCall(PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0)); 284 PetscFunctionReturn(0); 285 } 286 287 /*@ 288 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order. 289 290 Collective on dm 291 292 Input Parameters: 293 + dm - The distributed DMPlex 294 - nv - The natural Vec 295 296 Output Parameters: 297 . gv - The global Vec 298 299 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 300 301 Level: intermediate 302 303 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalEnd()` 304 @*/ 305 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv) 306 { 307 const PetscScalar *inarray; 308 PetscScalar *outarray; 309 PetscMPIInt size; 310 311 PetscFunctionBegin; 312 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 313 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 314 if (dm->sfNatural) { 315 /* We only have access to the SF that goes from Global to Natural. 316 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM. 317 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */ 318 PetscCall(VecZeroEntries(gv)); 319 PetscCall(VecGetArray(gv, &outarray)); 320 PetscCall(VecGetArrayRead(nv, &inarray)); 321 PetscCall(PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 322 PetscCall(VecRestoreArrayRead(nv, &inarray)); 323 PetscCall(VecRestoreArray(gv, &outarray)); 324 } else if (size == 1) { 325 PetscCall(VecCopy(nv, gv)); 326 } 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."); 327 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 328 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0)); 329 PetscFunctionReturn(0); 330 } 331 332 /*@ 333 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order. 334 335 Collective on dm 336 337 Input Parameters: 338 + dm - The distributed DMPlex 339 - nv - The natural Vec 340 341 Output Parameters: 342 . gv - The global Vec 343 344 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 345 346 Level: intermediate 347 348 .seealso: `DMPlexDistribute()`, `DMPlexDistributeField()`, `DMPlexNaturalToGlobalBegin()`, `DMPlexGlobalToNaturalBegin()` 349 @*/ 350 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv) 351 { 352 const PetscScalar *inarray; 353 PetscScalar *outarray; 354 PetscMPIInt size; 355 356 PetscFunctionBegin; 357 PetscCall(PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 358 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size)); 359 if (dm->sfNatural) { 360 PetscCall(VecGetArrayRead(nv, &inarray)); 361 PetscCall(VecGetArray(gv, &outarray)); 362 PetscCall(PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM)); 363 PetscCall(VecRestoreArrayRead(nv, &inarray)); 364 PetscCall(VecRestoreArray(gv, &outarray)); 365 } else if (size == 1) { 366 } 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."); 367 else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute()."); 368 PetscCall(PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0)); 369 PetscFunctionReturn(0); 370 } 371