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