1:  #include <petsc/private/dmpleximpl.h>

  3: /*@
  4:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

  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

 11:   Output Parameters:
 12: . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering

 14:   Level: intermediate

 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;

 29:   PetscObjectGetComm((PetscObject) dm, &comm);
 30:   /* PetscPrintf(comm, "Point migration SF\n");
 31:    PetscSFView(sfMigration, 0); */
 32:   /* Create a new section from distributing the original section */
 33:   PetscSectionCreate(comm, &sectionDist);
 34:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
 35:   /* PetscPrintf(comm, "Distributed Section\n");
 36:    PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
 37:   DMSetDefaultSection(dm, sectionDist);
 38:   /* Get a pruned version of migration SF */
 39:   DMGetDefaultGlobalSection(dm, &gSection);
 40:   PetscSectionGetChart(gSection, &pStart, &pEnd);
 41:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
 42:     PetscInt dof, off;

 44:     PetscSectionGetDof(gSection, p, &dof);
 45:     PetscSectionGetOffset(gSection, p, &off);
 46:     if ((dof > 0) && (off >= 0)) ++ssize;
 47:   }
 48:   PetscMalloc1(ssize, &spoints);
 49:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
 50:     PetscInt dof, off;

 52:     PetscSectionGetDof(gSection, p, &dof);
 53:     PetscSectionGetOffset(gSection, p, &off);
 54:     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
 55:   }
 56:   PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
 57:   PetscFree(spoints);
 58:   /* PetscPrintf(comm, "Embedded SF\n");
 59:    PetscSFView(sfEmbed, 0); */
 60:   /* Create the SF for seq to natural */
 61:   DMGetGlobalVector(dm, &gv);
 62:   PetscSFCreateFromZero(comm, gv, &sfSeqToNatural);
 63:   DMRestoreGlobalVector(dm, &gv);
 64:   /* PetscPrintf(comm, "Seq-to-Natural SF\n");
 65:    PetscSFView(sfSeqToNatural, 0); */
 66:   /* Create the SF associated with this section */
 67:   DMGetPointSF(dm, &sf);
 68:   PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
 69:   PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
 70:   PetscFree(remoteOffsets);
 71:   PetscSFDestroy(&sfEmbed);
 72:   PetscSectionDestroy(&gLocSection);
 73:   PetscSectionDestroy(&sectionDist);
 74:   /* PetscPrintf(comm, "Field SF\n");
 75:    PetscSFView(sfField, 0); */
 76:   /* Invert the field SF so it's now from distributed to sequential */
 77:   PetscSFCreateInverseSF(sfField, &sfFieldInv);
 78:   PetscSFDestroy(&sfField);
 79:   /* PetscPrintf(comm, "Inverse Field SF\n");
 80:    PetscSFView(sfFieldInv, 0); */
 81:   /* Multiply the sfFieldInv with the */
 82:   PetscSFCompose(sfFieldInv, sfSeqToNatural, sfNatural);
 83:   PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
 84:   /* Clean up */
 85:   PetscSFDestroy(&sfFieldInv);
 86:   PetscSFDestroy(&sfSeqToNatural);
 87:   return(0);
 88: }

 90: /*@
 91:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

 93:   Collective on dm

 95:   Input Parameters:
 96: + dm - The distributed DMPlex
 97: - gv - The global Vec

 99:   Output Parameters:
100: . nv - Vec in the canonical ordering distributed over all processors associated with gv

102:   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().

104:   Level: intermediate

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;

115:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
116:   if (dm->sfNatural) {
117:     VecGetArray(nv, &outarray);
118:     VecGetArrayRead(gv, &inarray);
119:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
120:     VecRestoreArrayRead(gv, &inarray);
121:     VecRestoreArray(nv, &outarray);
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:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
124:   return(0);
125: }

127: /*@
128:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

130:   Collective on dm

132:   Input Parameters:
133: + dm - The distributed DMPlex
134: - gv - The global Vec

136:   Output Parameters:
137: . nv - The natural Vec

139:   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().

141:   Level: intermediate

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;

152:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
153:   if (dm->sfNatural) {
154:     VecGetArrayRead(gv, &inarray);
155:     VecGetArray(nv, &outarray);
156:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
157:     VecRestoreArrayRead(gv, &inarray);
158:     VecRestoreArray(nv, &outarray);
159:   }
160:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
161:   return(0);
162: }

164: /*@
165:   DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.

167:   Collective on dm

169:   Input Parameters:
170: + dm - The distributed DMPlex
171: - nv - The natural Vec

173:   Output Parameters:
174: . gv - The global Vec

176:   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().

178:   Level: intermediate

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;

189:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
190:   if (dm->sfNatural) {
191:     VecGetArray(gv, &outarray);
192:     VecGetArrayRead(nv, &inarray);
193:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
194:     VecRestoreArrayRead(nv, &inarray);
195:     VecRestoreArray(gv, &outarray);
196:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMPlexSetUseNaturalSF() before DMPlexDistribute().\n");
197:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
198:   return(0);
199: }

201: /*@
202:   DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.

204:   Collective on dm

206:   Input Parameters:
207: + dm - The distributed DMPlex
208: - nv - The natural Vec

210:   Output Parameters:
211: . gv - The global Vec

213:   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().

215:   Level: intermediate

217: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
218:  @*/
219: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
220: {
221:   const PetscScalar *inarray;
222:   PetscScalar       *outarray;
223:   PetscErrorCode     ierr;

226:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
227:   if (dm->sfNatural) {
228:     VecGetArrayRead(nv, &inarray);
229:     VecGetArray(gv, &outarray);
230:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
231:     VecRestoreArrayRead(nv, &inarray);
232:     VecRestoreArray(gv, &outarray);
233:   }
234:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
235:   return(0);
236: }