xref: /petsc/src/dm/impls/plex/plexnatural.c (revision 2b8d69ca7ea5fe9190df62c1dce3bbd66fce84dd) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #undef __FUNCT__
4 #define __FUNCT__ "DMPlexCreateGlobalToNaturalSF"
5 /*@
6   DMPlexCreateGlobalToNaturalSF - Creates the SF for mapping Global Vec to the Natural Vec
7 
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
12 
13   Output Parameters:
14 . sfNatural - PetscSF for mapping the Vec in PETSc ordering to the canonical ordering
15 
16   Level: intermediate
17 
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;
28   PetscErrorCode ierr;
29 
30   PetscFunctionBegin;
31   ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr);
32   /* ierr = PetscPrintf(comm, "Point migration SF\n");CHKERRQ(ierr);
33    ierr = PetscSFView(sfMigration, 0);CHKERRQ(ierr); */
34   /* Create a new section from distributing the original section */
35   ierr = PetscSectionCreate(comm, &sectionDist);CHKERRQ(ierr);
36   ierr = PetscSFDistributeSection(sfMigration, section, &remoteOffsets, sectionDist);CHKERRQ(ierr);
37   /* ierr = PetscPrintf(comm, "Distributed Section\n");CHKERRQ(ierr);
38    ierr = PetscSectionView(sectionDist, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
39   ierr = DMSetDefaultSection(dm, sectionDist);CHKERRQ(ierr);
40   /* Get a pruned version of migration SF */
41   ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr);
42   ierr = PetscSectionGetChart(gSection, &pStart, &pEnd);CHKERRQ(ierr);
43   for (p = pStart, ssize = 0; p < pEnd; ++p) {
44     PetscInt dof, off;
45 
46     ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
47     ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
48     if ((dof > 0) && (off >= 0)) ++ssize;
49   }
50   ierr = PetscMalloc1(ssize, &spoints);CHKERRQ(ierr);
51   for (p = pStart, ssize = 0; p < pEnd; ++p) {
52     PetscInt dof, off;
53 
54     ierr = PetscSectionGetDof(gSection, p, &dof);CHKERRQ(ierr);
55     ierr = PetscSectionGetOffset(gSection, p, &off);CHKERRQ(ierr);
56     if ((dof > 0) && (off >= 0)) spoints[ssize++] = p;
57   }
58   ierr = PetscSFCreateEmbeddedLeafSF(sfMigration, ssize, spoints, &sfEmbed);CHKERRQ(ierr);
59   ierr = PetscFree(spoints);CHKERRQ(ierr);
60   /* ierr = PetscPrintf(comm, "Embedded SF\n");CHKERRQ(ierr);
61    ierr = PetscSFView(sfEmbed, 0);CHKERRQ(ierr); */
62   /* Create the SF for seq to natural */
63   ierr = DMGetGlobalVector(dm, &gv);CHKERRQ(ierr);
64   ierr = PetscSFCreateFromZero(comm, gv, &sfSeqToNatural);CHKERRQ(ierr);
65   ierr = DMRestoreGlobalVector(dm, &gv);CHKERRQ(ierr);
66   /* ierr = PetscPrintf(comm, "Seq-to-Natural SF\n");CHKERRQ(ierr);
67    ierr = PetscSFView(sfSeqToNatural, 0);CHKERRQ(ierr); */
68   /* Create the SF associated with this section */
69   ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
70   ierr = PetscSectionCreateGlobalSection(sectionDist, sf, PETSC_FALSE, PETSC_TRUE, &gLocSection);CHKERRQ(ierr);
71   ierr = PetscSFCreateSectionSF(sfEmbed, section, remoteOffsets, gLocSection, &sfField);CHKERRQ(ierr);
72   ierr = PetscFree(remoteOffsets);CHKERRQ(ierr);
73   ierr = PetscSFDestroy(&sfEmbed);CHKERRQ(ierr);
74   ierr = PetscSectionDestroy(&gLocSection);CHKERRQ(ierr);
75   ierr = PetscSectionDestroy(&sectionDist);CHKERRQ(ierr);
76   /* ierr = PetscPrintf(comm, "Field SF\n");CHKERRQ(ierr);
77    ierr = PetscSFView(sfField, 0);CHKERRQ(ierr); */
78   /* Invert the field SF so it's now from distributed to sequential */
79   ierr = PetscSFCreateInverseSF(sfField, &sfFieldInv);CHKERRQ(ierr);
80   ierr = PetscSFDestroy(&sfField);CHKERRQ(ierr);
81   /* ierr = PetscPrintf(comm, "Inverse Field SF\n");CHKERRQ(ierr);
82    ierr = PetscSFView(sfFieldInv, 0);CHKERRQ(ierr); */
83   /* Multiply the sfFieldInv with the */
84   ierr = PetscSFCompose(sfFieldInv, sfSeqToNatural, sfNatural);CHKERRQ(ierr);
85   ierr = PetscObjectViewFromOptions((PetscObject) *sfNatural, NULL, "-globaltonatural_sf_view");CHKERRQ(ierr);
86   /* Clean up */
87   ierr = PetscSFDestroy(&sfFieldInv);CHKERRQ(ierr);
88   ierr = PetscSFDestroy(&sfSeqToNatural);CHKERRQ(ierr);
89   PetscFunctionReturn(0);
90 }
91 
92 #undef __FUNCT__
93 #define __FUNCT__ "DMPlexGlobalToNaturalBegin"
94 /*@
95   DMPlexGlobalToNaturalBegin - Rearranges a global Vector in the natural order.
96 
97   Collective on dm
98 
99   Input Parameters:
100 + dm - The distributed DMPlex
101 - gv - The global Vec
102 
103   Output Parameters:
104 . nv - Vec in the canonical ordering distributed over all processors associated with gv
105 
106   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().
107 
108   Level: intermediate
109 
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;
117 
118   PetscFunctionBegin;
119   ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
120   if (dm->sfNatural) {
121     ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);
122     ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
123     ierr = PetscSFBcastBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
124     ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
125     ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
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   ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalBegin,dm,0,0,0);CHKERRQ(ierr);
128   PetscFunctionReturn(0);
129 }
130 
131 #undef __FUNCT__
132 #define __FUNCT__ "DMPlexGlobalToNaturalEnd"
133 /*@
134   DMPlexGlobalToNaturalEnd - Rearranges a global Vector in the natural order.
135 
136   Collective on dm
137 
138   Input Parameters:
139 + dm - The distributed DMPlex
140 - gv - The global Vec
141 
142   Output Parameters:
143 . nv - The natural Vec
144 
145   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().
146 
147   Level: intermediate
148 
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;
156 
157   PetscFunctionBegin;
158   ierr = PetscLogEventBegin(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
159   if (dm->sfNatural) {
160     ierr = VecGetArrayRead(gv, &inarray);CHKERRQ(ierr);
161     ierr = VecGetArray(nv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr);
162     ierr = PetscSFBcastEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray);CHKERRQ(ierr);
163     ierr = VecRestoreArrayRead(gv, &inarray);CHKERRQ(ierr);
164     ierr = VecRestoreArray(nv, &outarray);CHKERRQ(ierr);
165   }
166   ierr = PetscLogEventEnd(DMPLEX_GlobalToNaturalEnd,dm,0,0,0);CHKERRQ(ierr);
167   PetscFunctionReturn(0);
168 }
169 
170 #undef __FUNCT__
171 #define __FUNCT__ "DMPlexNaturalToGlobalBegin"
172 /*@
173   DMPlexNaturalToGlobalBegin - Rearranges a Vector in the natural order to the Global order.
174 
175   Collective on dm
176 
177   Input Parameters:
178 + dm - The distributed DMPlex
179 - nv - The natural Vec
180 
181   Output Parameters:
182 . gv - The global Vec
183 
184   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().
185 
186   Level: intermediate
187 
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;
195 
196   PetscFunctionBegin;
197   ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
198   if (dm->sfNatural) {
199     ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);
200     ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
201     ierr = PetscSFReduceBegin(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
202     ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
203     ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
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   ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalBegin,dm,0,0,0);CHKERRQ(ierr);
206   PetscFunctionReturn(0);
207 }
208 
209 #undef __FUNCT__
210 #define __FUNCT__ "DMPlexNaturalToGlobalEnd"
211 /*@
212   DMPlexNaturalToGlobalEnd - Rearranges a Vector in the natural order to the Global order.
213 
214   Collective on dm
215 
216   Input Parameters:
217 + dm - The distributed DMPlex
218 - nv - The natural Vec
219 
220   Output Parameters:
221 . gv - The global Vec
222 
223   Note: The user must call DMPlexSetUseNaturalSF(dm, PETSC_TRUE) before DMPlexDistribute().
224 
225   Level: intermediate
226 
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;
234 
235   PetscFunctionBegin;
236   ierr = PetscLogEventBegin(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
237   if (dm->sfNatural) {
238     ierr = VecGetArrayRead(nv, &inarray);CHKERRQ(ierr);
239     ierr = VecGetArray(gv, &outarray);CHKERRQ(ierr);CHKERRQ(ierr);
240     ierr = PetscSFReduceEnd(dm->sfNatural, MPIU_SCALAR, (PetscScalar *) inarray, outarray, MPI_SUM);CHKERRQ(ierr);
241     ierr = VecRestoreArrayRead(nv, &inarray);CHKERRQ(ierr);
242     ierr = VecRestoreArray(gv, &outarray);CHKERRQ(ierr);
243   }
244   ierr = PetscLogEventEnd(DMPLEX_NaturalToGlobalEnd,dm,0,0,0);CHKERRQ(ierr);
245   PetscFunctionReturn(0);
246 }
247