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