xref: /petsc/src/dm/impls/plex/plexnatural.c (revision 0e03b746557e2551025fde0294144c0532d12f68) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 /*@
4   DMPlexSetMigrationSF - Sets the SF for migrating from a parent DM into this DM
5 
6   Input Parameters:
7 + dm        - The DM
8 - naturalSF - The PetscSF
9 
10   Note: It is necessary to call this in order to have DMCreateSubDM() or DMCreateSuperDM() build the Global-To-Natural map
11 
12   Level: intermediate
13 
14 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexGetMigrationSF()
15 @*/
16 PetscErrorCode DMPlexSetMigrationSF(DM dm, PetscSF migrationSF)
17 {
18   PetscErrorCode ierr;
19   PetscFunctionBegin;
20   dm->sfMigration = migrationSF;
21   ierr = PetscObjectReference((PetscObject) migrationSF);CHKERRQ(ierr);
22   PetscFunctionReturn(0);
23 }
24 
25 /*@
26   DMPlexGetMigrationSF - Gets the SF for migrating from a parent DM into this DM
27 
28   Input Parameter:
29 . dm          - The DM
30 
31   Output Parameter:
32 . migrationSF - The PetscSF
33 
34   Level: intermediate
35 
36 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateMigrationSF(), DMPlexSetMigrationSF
37 @*/
38 PetscErrorCode DMPlexGetMigrationSF(DM dm, PetscSF *migrationSF)
39 {
40   PetscFunctionBegin;
41   *migrationSF = dm->sfMigration;
42   PetscFunctionReturn(0);
43 }
44 
45 /*@
46   DMPlexSetGlobalToNaturalSF - Sets the SF for mapping Global Vec to the Natural Vec
47 
48   Input Parameters:
49 + dm          - The DM
50 - naturalSF   - The PetscSF
51 
52   Level: intermediate
53 
54 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexGetGlobaltoNaturalSF()
55 @*/
56 PetscErrorCode DMPlexSetGlobalToNaturalSF(DM dm, PetscSF naturalSF)
57 {
58   PetscErrorCode ierr;
59   PetscFunctionBegin;
60   dm->sfNatural = naturalSF;
61   ierr = PetscObjectReference((PetscObject) naturalSF);CHKERRQ(ierr);
62   dm->useNatural = PETSC_TRUE;
63   PetscFunctionReturn(0);
64 }
65 
66 /*@
67   DMPlexGetGlobalToNaturalSF - Gets the SF for mapping Global Vec to the Natural Vec
68 
69   Input Parameter:
70 . dm          - The DM
71 
72   Output Parameter:
73 . naturalSF   - The PetscSF
74 
75   Level: intermediate
76 
77 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexCreateGlobalToNaturalSF(), DMPlexSetGlobaltoNaturalSF
78 @*/
79 PetscErrorCode DMPlexGetGlobalToNaturalSF(DM dm, PetscSF *naturalSF)
80 {
81   PetscFunctionBegin;
82   *naturalSF = dm->sfNatural;
83   PetscFunctionReturn(0);
84 }
85 
86 /*@
87   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec
88 
89   Input Parameters:
90 + dm          - The DM
91 . section     - The PetscSection describing the Vec before the mesh was distributed
92 - sfMigration - The PetscSF used to distribute the mesh
93 
94   Output Parameter:
95 . sfNatural   - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering
96 
97   Note: This is not typically called by the user.
98 
99   Level: intermediate
100 
101 .seealso: DMPlexDistribute(), DMPlexDistributeField()
102  @*/
103 PetscErrorCode DMPlexCreateGlobalToNaturalSF(DM dm, PetscSection section, PetscSF sfMigration, PetscSF *sfNatural)
104 {
105   MPI_Comm       comm;
106   Vec            gv;
107   PetscSF        sf, sfEmbed, sfSeqToNatural, sfField, sfFieldInv;
108   PetscSection   gSection, sectionDist, gLocSection;
109   PetscInt      *spoints, *remoteOffsets;
110   PetscInt       ssize, pStart, pEnd, p;
111   PetscLayout    map;
112   PetscErrorCode ierr;
113 
114   PetscFunctionBegin;
115   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
116   /* ierr = PetscPrintf(comm, "Point migration SF\n");CHKERRQ(ierr);
117    ierr = PetscSFView(sfMigration, 0);CHKERRQ(ierr); */
118   /* Create a new section from distributing the original section */
119   ierr = PetscSectionCreate(comm, &sectionDist);CHKERRQ(ierr);
120   ierr = PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);CHKERRQ(ierr);
121   /* ierr = PetscPrintf(comm, "Distributed Section\n");CHKERRQ(ierr);
122    ierr = PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
123   ierr = DMSetLocalSection(dm, sectionDist);CHKERRQ(ierr);
124   /* Get a pruned version of migration SF */
125   ierr = DMGetGlobalSection(dm, &gSection);CHKERRQ(ierr);
126   ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr);
127   for (p = pStart, ssize = 0; p < pEnd; ++p) {
128     PetscInt dof, off;
129 
130     ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
131     ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
132     if ((dof > 0) && (off >= 0)) ++ssize;
133   }
134   ierr = PetscMalloc1(ssize, &spoints);CHKERRQ(ierr);
135   for (p = pStart, ssize = 0; p < pEnd; ++p) {
136     PetscInt dof, off;
137 
138     ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
139     ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
140     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
141   }
142   ierr = PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);CHKERRQ(ierr);
143   ierr = PetscFree(spoints);CHKERRQ(ierr);
144   /* ierr = PetscPrintf(comm, "Embedded SF\n");CHKERRQ(ierr);
145    ierr = PetscSFView(sfEmbed, 0);CHKERRQ(ierr); */
146   /* Create the SF for seq to natural */
147   ierr = DMGetGlobalVector(dm, &gv);CHKERRQ(ierr);
148   ierr = VecGetLayout(gv,&map);CHKERRQ(ierr);
149   /* Note that entries of gv are leaves in sfSeqToNatural, entries of the seq vec are roots */
150   ierr = PetscSFCreate(comm, &sfSeqToNatural);CHKERRQ(ierr);
151   ierr = PetscSFSetGraphWithPattern(sfSeqToNatural, map, PETSCSF_PATTERN_GATHER);CHKERRQ(ierr);
152   ierr = DMRestoreGlobalVector(dm, &gv);CHKERRQ(ierr);
153   /* ierr = PetscPrintf(comm, "Seq-to-Natural SF\n");CHKERRQ(ierr);
154    ierr = PetscSFView(sfSeqToNatural, 0);CHKERRQ(ierr); */
155   /* Create the SF associated with this section */
156   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
157   ierr = PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);CHKERRQ(ierr);
158   ierr = PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);CHKERRQ(ierr);
159   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
160   ierr = PetscSFDestroy(&sfEmbed);CHKERRQ(ierr);
161   ierr = PetscSectionDestroy(&gLocSection);CHKERRQ(ierr);
162   ierr = PetscSectionDestroy(&sectionDist);CHKERRQ(ierr);
163   /* ierr = PetscPrintf(comm, "Field SF\n");CHKERRQ(ierr);
164    ierr = PetscSFView(sfField, 0);CHKERRQ(ierr); */
165   /* Invert the field SF so it's now from distributed to sequential */
166   ierr = PetscSFCreateInverseSF(sfField, &sfFieldInv);CHKERRQ(ierr);
167   ierr = PetscSFDestroy(&sfField);CHKERRQ(ierr);
168   /* ierr = PetscPrintf(comm, "Inverse Field SF\n");CHKERRQ(ierr);
169    ierr = PetscSFView(sfFieldInv, 0);CHKERRQ(ierr); */
170   /* Multiply the sfFieldInv with the */
171   ierr = PetscSFComposeInverse(sfFieldInv, sfSeqToNatural, sfNatural);CHKERRQ(ierr);
172   ierr = PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");CHKERRQ(ierr);
173   /* Clean up */
174   ierr = PetscSFDestroy(&sfFieldInv);CHKERRQ(ierr);
175   ierr = PetscSFDestroy(&sfSeqToNatural);CHKERRQ(ierr);
176   PetscFunctionReturn(0);
177 }
178 
179 /*@
180   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
181 
182   Collective on dm
183 
184   Input Parameters:
185 + dm - The distributed DMPlex
186 - gv - The global Vec
187 
188   Output Parameters:
189 . nv - Vec in the canonical ordering distributed over all processors associated with gv
190 
191   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
192 
193   Level: intermediate
194 
195 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalEnd()
196 @*/
197 PetscErrorCode DMPlexGlobalToNaturalBegin(DM dm, Vec gv, Vec nv)
198 {
199   const PetscScalar *inarray;
200   PetscScalar       *outarray;
201   PetscMPIInt        size;
202   PetscErrorCode     ierr;
203 
204   PetscFunctionBegin;
205   ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
206   ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
207   if (dm->sfNatural) {
208     ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);
209     ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
210     ierr = PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
211     ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
212     ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
213   } else if (size == 1) {
214     ierr = VecCopy(nv, gv);CHKERRQ(ierr);
215   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
216   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
217   ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
218   PetscFunctionReturn(0);
219 }
220 
221 /*@
222   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
223 
224   Collective on dm
225 
226   Input Parameters:
227 + dm - The distributed DMPlex
228 - gv - The global Vec
229 
230   Output Parameters:
231 . nv - The natural Vec
232 
233   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
234 
235   Level: intermediate
236 
237  .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
238  @*/
239 PetscErrorCode DMPlexGlobalToNaturalEnd(DM dm, Vec gv, Vec nv)
240 {
241   const PetscScalar *inarray;
242   PetscScalar       *outarray;
243   PetscMPIInt        size;
244   PetscErrorCode     ierr;
245 
246   PetscFunctionBegin;
247   ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
248   ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
249   if (dm->sfNatural) {
250     ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
251     ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);
252     ierr = PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
253     ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
254     ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
255   } else if (size == 1) {
256   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
257   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
258   ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
259   PetscFunctionReturn(0);
260 }
261 
262 /*@
263   DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
264 
265   Collective on dm
266 
267   Input Parameters:
268 + dm - The distributed DMPlex
269 - nv - The natural Vec
270 
271   Output Parameters:
272 . gv - The global Vec
273 
274   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
275 
276   Level: intermediate
277 
278 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(),DMPlexGlobalToNaturalEnd()
279 @*/
280 PetscErrorCode DMPlexNaturalToGlobalBegin(DM dm, Vec nv, Vec gv)
281 {
282   const PetscScalar *inarray;
283   PetscScalar       *outarray;
284   PetscMPIInt        size;
285   PetscErrorCode     ierr;
286 
287   PetscFunctionBegin;
288   ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
289   ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
290   if (dm->sfNatural) {
291     /* We only have acces to the SF that goes from Global to Natural.
292        Instead of inverting dm->sfNatural, we can call PetscSFReduceBegin/End with MPI_Op MPI_SUM.
293        Here the SUM really does nothing since sfNatural is one to one, as long as gV is set to zero first. */
294     ierr = VecZeroEntries(gv);CHKERRQ(ierr);
295     ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);
296     ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
297     ierr = PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
298     ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
299     ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
300   } else if (size == 1) {
301     ierr = VecCopy(nv, gv);CHKERRQ(ierr);
302   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
303   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
304   ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
305   PetscFunctionReturn(0);
306 }
307 
308 /*@
309   DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
310 
311   Collective on dm
312 
313   Input Parameters:
314 + dm - The distributed DMPlex
315 - nv - The natural Vec
316 
317   Output Parameters:
318 . gv - The global Vec
319 
320   Note: The user must call DMSetUseNatural(dm, PETSC_TRUE) before DMPlexDistribute().
321 
322   Level: intermediate
323 
324 .seealso: DMPlexDistribute(), DMPlexDistributeField(), DMPlexNaturalToGlobalBegin(), DMPlexGlobalToNaturalBegin()
325  @*/
326 PetscErrorCode DMPlexNaturalToGlobalEnd(DM dm, Vec nv, Vec gv)
327 {
328   const PetscScalar *inarray;
329   PetscScalar       *outarray;
330   PetscMPIInt        size;
331   PetscErrorCode     ierr;
332 
333   PetscFunctionBegin;
334   ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
335   ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
336   if (dm->sfNatural) {
337     ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
338     ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);
339     ierr = PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
340     ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
341     ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
342   } else if (size == 1) {
343   } else if (dm->useNatural) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "DM global to natural SF not present.\nIf DMPlexDistribute() was called, report to petsc-maint@mcs.anl.gov.\n");
344   else SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM global to natural SF was not created.\nYou must call DMSetUseNatural() before DMPlexDistribute().\n");
345   ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
346   PetscFunctionReturn(0);
347 }
348