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