Actual source code: plexnatural.c
petsc-3.7.3 2016-08-01
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, §ionDist);
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(§ionDist);
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: }