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 PetscErrorCode ierr; 19 PetscFunctionBegin; 20 dm->sfMigration = migrationSF; 21 ierr = PetscObjectReference((PetscObject) migrationSF);CHKERRQ(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 /*@ 26 DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM 27 28 Input Parameter: 29 . dm - The DM 30 31 Output Parameter: 32 . migrationSF - The PetscSF 33 34 Level: intermediate 35 36 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF 37 @*/ 38 PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF) 39 { 40 PetscFunctionBegin; 41 *migrationSF = dm->sfMigration; 42 PetscFunctionReturn(0); 43 } 44 45 /*@ 46 DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec 47 48 Input Parameters: 49 + dm - The DM 50 - naturalSF - The PetscSF 51 52 Level: intermediate 53 54 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF() 55 @*/ 56 PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF) 57 { 58 PetscErrorCode ierr; 59 PetscFunctionBegin; 60 dm->sfNatural = naturalSF; 61 ierr = PetscObjectReference((PetscObject) naturalSF);CHKERRQ(ierr); 62 dm->useNatural = PETSC_TRUE; 63 PetscFunctionReturn(0); 64 } 65 66 /*@ 67 DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec 68 69 Input Parameter: 70 . dm - The DM 71 72 Output Parameter: 73 . naturalSF - The PetscSF 74 75 Level: intermediate 76 77 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF 78 @*/ 79 PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF) 80 { 81 PetscFunctionBegin; 82 *naturalSF = dm->sfNatural; 83 PetscFunctionReturn(0); 84 } 85 86 /*@ 87 DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec 88 89 Input Parameters: 90 + dm - The DM 91 . section - The PetscSection describing the Vec before the mesh was distributed 92 - sfMigration - The PetscSF used to distribute the mesh 93 94 Output Parameter: 95 . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering 96 97 Note: This is not typically called by the user. 98 99 Level: intermediate 100 101 .seealso: DMPlexDistribute(), DMPlexDistributeField() 102 @*/ 103 PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural) 104 { 105 MPI_Comm comm; 106 Vec gv; 107 PetscSF sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv; 108 PetscSection gSection, sectionDist, gLocSection; 109 PetscInt *spoints, *remoteOffsets; 110 PetscInt ssize, pStart, pEnd, p; 111 PetscLayout map; 112 PetscErrorCode ierr; 113 114 PetscFunctionBegin; 115 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 116 /* ierr = PetscPrintf(comm, "Point migration SF\n");CHKERRQ(ierr); 117 ierr = PetscSFView(sfMigration, 0);CHKERRQ(ierr); */ 118 /* Create a new section from distributing the original section */ 119 ierr = PetscSectionCreate(comm, §ionDist);CHKERRQ(ierr); 120 ierr = PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);CHKERRQ(ierr); 121 /* ierr = PetscPrintf(comm, "Distributed Section\n");CHKERRQ(ierr); 122 ierr = PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 123 ierr = DMSetLocalSection(dm, sectionDist);CHKERRQ(ierr); 124 /* Get a pruned version of migration SF */ 125 ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr); 126 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr); 127 for (p = pStart, ssize = 0; p < pEnd; ++p) { 128 PetscInt dof, off; 129 130 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 131 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr); 132 if ((dof > 0) && (off >= 0)) ++ssize; 133 } 134 ierr = PetscMalloc1(ssize, &spoints);CHKERRQ(ierr); 135 for (p = pStart, ssize = 0; p < pEnd; ++p) { 136 PetscInt dof, off; 137 138 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 139 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr); 140 if ((dof > 0) && (off >= 0)) spoints[ssize++] = p; 141 } 142 ierr = PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);CHKERRQ(ierr); 143 ierr = PetscFree(spoints);CHKERRQ(ierr); 144 /* ierr = PetscPrintf(comm, "Embedded SF\n");CHKERRQ(ierr); 145 ierr = PetscSFView(sfEmbed, 0);CHKERRQ(ierr); */ 146 /* Create the SF for seq to natural */ 147 ierr = DMGetGlobalVector(dm, &gv);CHKERRQ(ierr); 148 ierr = VecGetLayout(gv,&map);CHKERRQ(ierr); 149 /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */ 150 ierr = PetscSFCreate(comm, &sfSeqToNatural);CHKERRQ(ierr); 151 ierr = PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);CHKERRQ(ierr); 152 ierr = DMRestoreGlobalVector(dm, &gv);CHKERRQ(ierr); 153 /* ierr = PetscPrintf(comm, "Seq-to-Natural SF\n");CHKERRQ(ierr); 154 ierr = PetscSFView(sfSeqToNatural, 0);CHKERRQ(ierr); */ 155 /* Create the SF associated with this section */ 156 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 157 ierr = PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);CHKERRQ(ierr); 158 ierr = PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);CHKERRQ(ierr); 159 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 160 ierr = PetscSFDestroy(&sfEmbed);CHKERRQ(ierr); 161 ierr = PetscSectionDestroy(&gLocSection);CHKERRQ(ierr); 162 ierr = PetscSectionDestroy(§ionDist);CHKERRQ(ierr); 163 /* ierr = PetscPrintf(comm, "Field SF\n");CHKERRQ(ierr); 164 ierr = PetscSFView(sfField, 0);CHKERRQ(ierr); */ 165 /* Invert the field SF so it's now from distributed to sequential */ 166 ierr = PetscSFCreateInverseSF(sfField, &sfFieldInv);CHKERRQ(ierr); 167 ierr = PetscSFDestroy(&sfField);CHKERRQ(ierr); 168 /* ierr = PetscPrintf(comm, "Inverse Field SF\n");CHKERRQ(ierr); 169 ierr = PetscSFView(sfFieldInv, 0);CHKERRQ(ierr); */ 170 /* Multiply the sfFieldInv with the */ 171 ierr = PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);CHKERRQ(ierr); 172 ierr = PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");CHKERRQ(ierr); 173 /* Clean up */ 174 ierr = PetscSFDestroy(&sfFieldInv);CHKERRQ(ierr); 175 ierr = PetscSFDestroy(&sfSeqToNatural);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 /*@ 180 DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order. 181 182 Collective on dm 183 184 Input Parameters: 185 + dm - The distributed DMPlex 186 - gv - The global Vec 187 188 Output Parameters: 189 . nv - Vec in the canonical ordering distributed over all processors associated with gv 190 191 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 192 193 Level: intermediate 194 195 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd() 196 @*/ 197 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv) 198 { 199 const PetscScalar *inarray; 200 PetscScalar *outarray; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr); 205 if (dm->sfNatural) { 206 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr); 207 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr); 208 ierr = PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr); 209 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr); 210 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr); 211 } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n"); 212 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 /*@ 217 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order. 218 219 Collective on dm 220 221 Input Parameters: 222 + dm - The distributed DMPlex 223 - gv - The global Vec 224 225 Output Parameters: 226 . nv - The natural Vec 227 228 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 229 230 Level: intermediate 231 232 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin() 233 @*/ 234 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv) 235 { 236 const PetscScalar *inarray; 237 PetscScalar *outarray; 238 PetscErrorCode ierr; 239 240 PetscFunctionBegin; 241 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr); 242 if (dm->sfNatural) { 243 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr); 244 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr); 245 ierr = PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr); 246 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr); 247 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr); 248 } 249 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /*@ 254 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order. 255 256 Collective on dm 257 258 Input Parameters: 259 + dm - The distributed DMPlex 260 - nv - The natural Vec 261 262 Output Parameters: 263 . gv - The global Vec 264 265 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 266 267 Level: intermediate 268 269 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd() 270 @*/ 271 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv) 272 { 273 const PetscScalar *inarray; 274 PetscScalar *outarray; 275 PetscErrorCode ierr; 276 277 PetscFunctionBegin; 278 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr); 279 if (dm->sfNatural) { 280 /* We only have acces to the SF that goes from Global to Natural. 281 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM. 282 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */ 283 ierr = VecZeroEntries(gv);CHKERRQ(ierr); 284 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr); 285 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr); 286 ierr = PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr); 287 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr); 288 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr); 289 } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n"); 290 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 /*@ 295 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order. 296 297 Collective on dm 298 299 Input Parameters: 300 + dm - The distributed DMPlex 301 - nv - The natural Vec 302 303 Output Parameters: 304 . gv - The global Vec 305 306 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 307 308 Level: intermediate 309 310 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin() 311 @*/ 312 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv) 313 { 314 const PetscScalar *inarray; 315 PetscScalar *outarray; 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr); 320 if (dm->sfNatural) { 321 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr); 322 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr); 323 ierr = PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr); 324 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr); 325 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr); 326 } 327 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330