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 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 115 /* ierr = PetscPrintf(comm, "Point migration SF\n");CHKERRQ(ierr); 116 ierr = PetscSFView(sfMigration, 0);CHKERRQ(ierr); */ 117 /* Create a new section from distributing the original section */ 118 ierr = PetscSectionCreate(comm, §ionDist);CHKERRQ(ierr); 119 ierr = PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);CHKERRQ(ierr); 120 /* ierr = PetscPrintf(comm, "Distributed Section\n");CHKERRQ(ierr); 121 ierr = PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 122 ierr = DMSetSection(dm, sectionDist);CHKERRQ(ierr); 123 /* Get a pruned version of migration SF */ 124 ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr); 125 ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr); 126 for (p = pStart, ssize = 0; p < pEnd; ++p) { 127 PetscInt dof, off; 128 129 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 130 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr); 131 if ((dof > 0) && (off >= 0)) ++ssize; 132 } 133 ierr = PetscMalloc1(ssize, &spoints);CHKERRQ(ierr); 134 for (p = pStart, ssize = 0; p < pEnd; ++p) { 135 PetscInt dof, off; 136 137 ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr); 138 ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr); 139 if ((dof > 0) && (off >= 0)) spoints[ssize++] = p; 140 } 141 ierr = PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);CHKERRQ(ierr); 142 ierr = PetscFree(spoints);CHKERRQ(ierr); 143 /* ierr = PetscPrintf(comm, "Embedded SF\n");CHKERRQ(ierr); 144 ierr = PetscSFView(sfEmbed, 0);CHKERRQ(ierr); */ 145 /* Create the SF for seq to natural */ 146 ierr = DMGetGlobalVector(dm, &gv);CHKERRQ(ierr); 147 ierr = PetscSFCreateFromZero(comm, gv, &sfSeqToNatural);CHKERRQ(ierr); 148 ierr = DMRestoreGlobalVector(dm, &gv);CHKERRQ(ierr); 149 /* ierr = PetscPrintf(comm, "Seq-to-Natural SF\n");CHKERRQ(ierr); 150 ierr = PetscSFView(sfSeqToNatural, 0);CHKERRQ(ierr); */ 151 /* Create the SF associated with this section */ 152 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 153 ierr = PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);CHKERRQ(ierr); 154 ierr = PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);CHKERRQ(ierr); 155 ierr = PetscFree(remoteOffsets);CHKERRQ(ierr); 156 ierr = PetscSFDestroy(&sfEmbed);CHKERRQ(ierr); 157 ierr = PetscSectionDestroy(&gLocSection);CHKERRQ(ierr); 158 ierr = PetscSectionDestroy(§ionDist);CHKERRQ(ierr); 159 /* ierr = PetscPrintf(comm, "Field SF\n");CHKERRQ(ierr); 160 ierr = PetscSFView(sfField, 0);CHKERRQ(ierr); */ 161 /* Invert the field SF so it's now from distributed to sequential */ 162 ierr = PetscSFCreateInverseSF(sfField, &sfFieldInv);CHKERRQ(ierr); 163 ierr = PetscSFDestroy(&sfField);CHKERRQ(ierr); 164 /* ierr = PetscPrintf(comm, "Inverse Field SF\n");CHKERRQ(ierr); 165 ierr = PetscSFView(sfFieldInv, 0);CHKERRQ(ierr); */ 166 /* Multiply the sfFieldInv with the */ 167 ierr = PetscSFCompose(sfFieldInv, sfSeqToNatural, sfNatural);CHKERRQ(ierr); 168 ierr = PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");CHKERRQ(ierr); 169 /* Clean up */ 170 ierr = PetscSFDestroy(&sfFieldInv);CHKERRQ(ierr); 171 ierr = PetscSFDestroy(&sfSeqToNatural);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 /*@ 176 DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order. 177 178 Collective on dm 179 180 Input Parameters: 181 + dm - The distributed DMPlex 182 - gv - The global Vec 183 184 Output Parameters: 185 . nv - Vec in the canonical ordering distributed over all processors associated with gv 186 187 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 188 189 Level: intermediate 190 191 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd() 192 @*/ 193 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv) 194 { 195 const PetscScalar *inarray; 196 PetscScalar *outarray; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr); 201 if (dm->sfNatural) { 202 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr); 203 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr); 204 ierr = PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr); 205 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr); 206 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr); 207 } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n"); 208 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 /*@ 213 DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order. 214 215 Collective on dm 216 217 Input Parameters: 218 + dm - The distributed DMPlex 219 - gv - The global Vec 220 221 Output Parameters: 222 . nv - The natural Vec 223 224 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 225 226 Level: intermediate 227 228 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin() 229 @*/ 230 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv) 231 { 232 const PetscScalar *inarray; 233 PetscScalar *outarray; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr); 238 if (dm->sfNatural) { 239 ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr); 240 ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr); 241 ierr = PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr); 242 ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr); 243 ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr); 244 } 245 ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 /*@ 250 DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order. 251 252 Collective on dm 253 254 Input Parameters: 255 + dm - The distributed DMPlex 256 - nv - The natural Vec 257 258 Output Parameters: 259 . gv - The global Vec 260 261 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 262 263 Level: intermediate 264 265 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd() 266 @*/ 267 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv) 268 { 269 const PetscScalar *inarray; 270 PetscScalar *outarray; 271 PetscErrorCode ierr; 272 273 PetscFunctionBegin; 274 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr); 275 if (dm->sfNatural) { 276 /* We only have acces to the SF that goes from Global to Natural. 277 Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM. 278 Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */ 279 ierr = VecZeroEntries(gv);CHKERRQ(ierr); 280 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr); 281 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr); 282 ierr = PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr); 283 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr); 284 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr); 285 } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n"); 286 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 /*@ 291 DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order. 292 293 Collective on dm 294 295 Input Parameters: 296 + dm - The distributed DMPlex 297 - nv - The natural Vec 298 299 Output Parameters: 300 . gv - The global Vec 301 302 Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute(). 303 304 Level: intermediate 305 306 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin() 307 @*/ 308 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv) 309 { 310 const PetscScalar *inarray; 311 PetscScalar *outarray; 312 PetscErrorCode ierr; 313 314 PetscFunctionBegin; 315 ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr); 316 if (dm->sfNatural) { 317 ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr); 318 ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr); 319 ierr = PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr); 320 ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr); 321 ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr); 322 } 323 ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326