1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  5: /*@
  6:   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec

  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

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

 16:   Level: intermediate

 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;

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

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

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

 94: /*@
 95:   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.

 97:   Collective on dm

 99:   Input Parameters:
100: + dm - The distributed DMPlex
101: - gv - The global Vec

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

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

108:   Level: intermediate

110: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
111: @*/
112: PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
113: {
114:   const PetscScalar *inarray;
115:   PetscScalar       *outarray;
116:   PetscErrorCode     ierr;

119:   PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
120:   if (dm->sfNatural) {
121:     VecGetArray(nv, &outarray);
122:     VecGetArrayRead(gv, &inarray);
123:     PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
124:     VecRestoreArrayRead(gv, &inarray);
125:     VecRestoreArray(nv, &outarray);
126:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
127:   PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);
128:   return(0);
129: }

133: /*@
134:   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.

136:   Collective on dm

138:   Input Parameters:
139: + dm - The distributed DMPlex
140: - gv - The global Vec

142:   Output Parameters:
143: . nv - The natural Vec

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

147:   Level: intermediate

149:  .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
150:  @*/
151: PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
152: {
153:   const PetscScalar *inarray;
154:   PetscScalar       *outarray;
155:   PetscErrorCode     ierr;

158:   PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
159:   if (dm->sfNatural) {
160:     VecGetArrayRead(gv, &inarray);
161:     VecGetArray(nv, &outarray);
162:     PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);
163:     VecRestoreArrayRead(gv, &inarray);
164:     VecRestoreArray(nv, &outarray);
165:   }
166:   PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);
167:   return(0);
168: }

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

175:   Collective on dm

177:   Input Parameters:
178: + dm - The distributed DMPlex
179: - nv - The natural Vec

181:   Output Parameters:
182: . gv - The global Vec

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

186:   Level: intermediate

188: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
189: @*/
190: PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
191: {
192:   const PetscScalar *inarray;
193:   PetscScalar       *outarray;
194:   PetscErrorCode     ierr;

197:   PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
198:   if (dm->sfNatural) {
199:     VecGetArray(gv, &outarray);
200:     VecGetArrayRead(nv, &inarray);
201:     PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
202:     VecRestoreArrayRead(nv, &inarray);
203:     VecRestoreArray(gv, &outarray);
204:   } else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMPlexSetUseNaturalSF() before DMPlexDistribute().\n");
205:   PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);
206:   return(0);
207: }

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

214:   Collective on dm

216:   Input Parameters:
217: + dm - The distributed DMPlex
218: - nv - The natural Vec

220:   Output Parameters:
221: . gv - The global Vec

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

225:   Level: intermediate

227: .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
228:  @*/
229: PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
230: {
231:   const PetscScalar *inarray;
232:   PetscScalar       *outarray;
233:   PetscErrorCode     ierr;

236:   PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
237:   if (dm->sfNatural) {
238:     VecGetArrayRead(nv, &inarray);
239:     VecGetArray(gv, &outarray);
240:     PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);
241:     VecRestoreArrayRead(nv, &inarray);
242:     VecRestoreArray(gv, &outarray);
243:   }
244:   PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);
245:   return(0);
246: }