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

  3: /*@
  4:   DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM

  6: + dm          - The DM
  7: . naturalSF   - The PetscSF
  8:   Level: intermediate

 10: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexGetMigrationSF()
 11: @*/
 12: PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
 13: {
 16:   dm->sfMigration = migrationSF;
 17:   PetscObjectReference((PetscObject) migrationSF);
 18:   return(0);
 19: }

 21: /*@
 22:   DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM

 24: + dm          - The DM
 25: . *migrationSF   - The PetscSF
 26:   Level: intermediate

 28: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
 29: @*/
 30: PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
 31: {
 33:   *migrationSF = dm->sfMigration;
 34:   return(0);
 35: }

 37: /*@
 38:   DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec

 40: + dm          - The DM
 41: . naturalSF   - The PetscSF
 42:   Level: intermediate

 44: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
 45: @*/
 46: PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
 47: {
 50:   dm->sfNatural = naturalSF;
 51:   PetscObjectReference((PetscObject) naturalSF);
 52:   dm->useNatural = PETSC_TRUE;
 53:   return(0);
 54: }

 56: /*@
 57:   DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec

 59: + dm          - The DM
 60: . *naturalSF   - The PetscSF
 61:   Level: intermediate

 63: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
 64: @*/
 65: PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
 66: {
 68:   *naturalSF = dm->sfNatural;
 69:   return(0);
 70: }

 72: /*@
 73:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

 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

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

 83:   Level: intermediate

 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;

 98:   PetscObjectGetComm((PetscObject) dm, &comm);
 99:   /* PetscPrintf(comm, "Point migration SF\n");
100:    PetscSFView(sfMigration, 0); */
101:   /* Create a new section from distributing the original section */
102:   PetscSectionCreate(comm, &sectionDist);
103:   PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);
104:   /* PetscPrintf(comm, "Distributed Section\n");
105:    PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD); */
106:   DMSetDefaultSection(dm, sectionDist);
107:   /* Get a pruned version of migration SF */
108:   DMGetDefaultGlobalSection(dm, &gSection);
109:   PetscSectionGetChart(gSection, &pStart, &pEnd);
110:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
111:     PetscInt dof, off;

113:     PetscSectionGetDof(gSection, p, &dof);
114:     PetscSectionGetOffset(gSection, p, &off);
115:     if ((dof > 0) && (off >= 0)) ++ssize;
116:   }
117:   PetscMalloc1(ssize, &spoints);
118:   for (p = pStart, ssize = 0; p < pEnd; ++p) {
119:     PetscInt dof, off;

121:     PetscSectionGetDof(gSection, p, &dof);
122:     PetscSectionGetOffset(gSection, p, &off);
123:     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
124:   }
125:   PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);
126:   PetscFree(spoints);
127:   /* PetscPrintf(comm, "Embedded SF\n");
128:    PetscSFView(sfEmbed, 0); */
129:   /* Create the SF for seq to natural */
130:   DMGetGlobalVector(dm, &gv);
131:   PetscSFCreateFromZero(comm, gv, &sfSeqToNatural);
132:   DMRestoreGlobalVector(dm, &gv);
133:   /* PetscPrintf(comm, "Seq-to-Natural SF\n");
134:    PetscSFView(sfSeqToNatural, 0); */
135:   /* Create the SF associated with this section */
136:   DMGetPointSF(dm, &sf);
137:   PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);
138:   PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);
139:   PetscFree(remoteOffsets);
140:   PetscSFDestroy(&sfEmbed);
141:   PetscSectionDestroy(&gLocSection);
142:   PetscSectionDestroy(&sectionDist);
143:   /* PetscPrintf(comm, "Field SF\n");
144:    PetscSFView(sfField, 0); */
145:   /* Invert the field SF so it's now from distributed to sequential */
146:   PetscSFCreateInverseSF(sfField, &sfFieldInv);
147:   PetscSFDestroy(&sfField);
148:   /* PetscPrintf(comm, "Inverse Field SF\n");
149:    PetscSFView(sfFieldInv, 0); */
150:   /* Multiply the sfFieldInv with the */
151:   PetscSFCompose(sfFieldInv, sfSeqToNatural, sfNatural);
152:   PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");
153:   /* Clean up */
154:   PetscSFDestroy(&sfFieldInv);
155:   PetscSFDestroy(&sfSeqToNatural);
156:   return(0);
157: }

159: /*@
160:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

162:   Collective on dm

164:   Input Parameters:
165: + dm - The distributed DMPlex
166: - gv - The global Vec

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

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

173:   Level: intermediate

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;

184:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
185:   if (dm->sfNatural) {
186:     VecGetArray(nv, &outarray);
187:     VecGetArrayRead(gv, &inarray);
188:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
189:     VecRestoreArrayRead(gv, &inarray);
190:     VecRestoreArray(nv, &outarray);
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:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
193:   return(0);
194: }

196: /*@
197:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

199:   Collective on dm

201:   Input Parameters:
202: + dm - The distributed DMPlex
203: - gv - The global Vec

205:   Output Parameters:
206: . nv - The natural Vec

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

210:   Level: intermediate

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;

221:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
222:   if (dm->sfNatural) {
223:     VecGetArrayRead(gv, &inarray);
224:     VecGetArray(nv, &outarray);
225:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
226:     VecRestoreArrayRead(gv, &inarray);
227:     VecRestoreArray(nv, &outarray);
228:   }
229:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
230:   return(0);
231: }

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

236:   Collective on dm

238:   Input Parameters:
239: + dm - The distributed DMPlex
240: - nv - The natural Vec

242:   Output Parameters:
243: . gv - The global Vec

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

247:   Level: intermediate

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;

258:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
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:     VecZeroEntries(gv);
264:     VecGetArray(gv, &outarray);
265:     VecGetArrayRead(nv, &inarray);
266:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
267:     VecRestoreArrayRead(nv, &inarray);
268:     VecRestoreArray(gv, &outarray);
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:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
271:   return(0);
272: }

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

277:   Collective on dm

279:   Input Parameters:
280: + dm - The distributed DMPlex
281: - nv - The natural Vec

283:   Output Parameters:
284: . gv - The global Vec

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

288:   Level: intermediate

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;

299:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
300:   if (dm->sfNatural) {
301:     VecGetArrayRead(nv, &inarray);
302:     VecGetArray(gv, &outarray);
303:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
304:     VecRestoreArrayRead(nv, &inarray);
305:     VecRestoreArray(gv, &outarray);
306:   }
307:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
308:   return(0);
309: }