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