xref: /petsc/src/dm/impls/plex/tests/ex64.c (revision 061e922f3926be00487707c73b78dd3d40309129) 
1 static char help[] = "Test FEM layout and GlobalToNaturalSF\n\n";
2 
3 /*
4   In order to see the vectors which are being tested, use
5 
6      -ua_vec_view -s_vec_view
7 */
8 
9 #include <petsc.h>
10 #include <exodusII.h>
11 
12 #include <petsc/private/dmpleximpl.h>
13 
main(int argc,char ** argv)14 int main(int argc, char **argv)
15 {
16   DM              dm, pdm, dmU, dmA, dmS, dmUA, dmUA2, *dmList;
17   DM              seqdmU, seqdmA, seqdmS, seqdmUA, seqdmUA2, *seqdmList;
18   Vec             X, U, A, S, UA, UA2;
19   Vec             seqX, seqU, seqA, seqS, seqUA, seqUA2;
20   IS              isU, isA, isS, isUA;
21   IS              seqisU, seqisA, seqisS, seqisUA;
22   PetscSection    section;
23   const PetscInt  fieldU     = 0;
24   const PetscInt  fieldA     = 2;
25   const PetscInt  fieldS     = 1;
26   const PetscInt  fieldUA[2] = {0, 2};
27   char            ifilename[PETSC_MAX_PATH_LEN];
28   IS              csIS;
29   const PetscInt *csID;
30   PetscInt       *pStartDepth, *pEndDepth;
31   PetscInt        order = 1;
32   PetscInt        sdim, d, pStart, pEnd, p, numCS, set;
33   PetscMPIInt     rank, size;
34   PetscSF         migrationSF;
35 
36   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
38   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
39   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex64");
40   PetscCall(PetscOptionsString("-i", "Filename to read", "ex64", ifilename, ifilename, sizeof(ifilename), NULL));
41   PetscOptionsEnd();
42 
43   /* Read the mesh from a file in any supported format */
44   PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, "ex64_plex", PETSC_TRUE, &dm));
45   PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
46   PetscCall(DMSetFromOptions(dm));
47   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
48   PetscCall(DMGetDimension(dm, &sdim));
49 
50   /* Create the main section containing all fields */
51   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &section));
52   PetscCall(PetscSectionSetNumFields(section, 3));
53   PetscCall(PetscSectionSetFieldName(section, fieldU, "U"));
54   PetscCall(PetscSectionSetFieldName(section, fieldA, "Alpha"));
55   PetscCall(PetscSectionSetFieldName(section, fieldS, "Sigma"));
56   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
57   PetscCall(PetscSectionSetChart(section, pStart, pEnd));
58   PetscCall(PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth));
59   for (d = 0; d <= sdim; ++d) PetscCall(DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d]));
60   /* Vector field U, Scalar field Alpha, Tensor field Sigma */
61   PetscCall(PetscSectionSetFieldComponents(section, fieldU, sdim));
62   PetscCall(PetscSectionSetFieldComponents(section, fieldA, 1));
63   PetscCall(PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2));
64 
65   /* Going through cell sets then cells, and setting up storage for the sections */
66   PetscCall(DMGetLabelSize(dm, "Cell Sets", &numCS));
67   PetscCall(DMGetLabelIdIS(dm, "Cell Sets", &csIS));
68   if (csIS) PetscCall(ISGetIndices(csIS, &csID));
69   for (set = 0; set < numCS; set++) {
70     IS              cellIS;
71     const PetscInt *cellID;
72     PetscInt        numCells, cell, closureSize, *closureA = NULL;
73 
74     PetscCall(DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells));
75     PetscCall(DMGetStratumIS(dm, "Cell Sets", csID[set], &cellIS));
76     if (numCells > 0) {
77       /* dof layout ordered by increasing height in the DAG: cell, face, edge, vertex */
78       PetscInt  dofUP1Tri[]  = {2, 0, 0};
79       PetscInt  dofAP1Tri[]  = {1, 0, 0};
80       PetscInt  dofUP2Tri[]  = {2, 2, 0};
81       PetscInt  dofAP2Tri[]  = {1, 1, 0};
82       PetscInt  dofUP1Quad[] = {2, 0, 0};
83       PetscInt  dofAP1Quad[] = {1, 0, 0};
84       PetscInt  dofUP2Quad[] = {2, 2, 2};
85       PetscInt  dofAP2Quad[] = {1, 1, 1};
86       PetscInt  dofS2D[]     = {0, 0, 3};
87       PetscInt  dofUP1Tet[]  = {3, 0, 0, 0};
88       PetscInt  dofAP1Tet[]  = {1, 0, 0, 0};
89       PetscInt  dofUP2Tet[]  = {3, 3, 0, 0};
90       PetscInt  dofAP2Tet[]  = {1, 1, 0, 0};
91       PetscInt  dofUP1Hex[]  = {3, 0, 0, 0};
92       PetscInt  dofAP1Hex[]  = {1, 0, 0, 0};
93       PetscInt  dofUP2Hex[]  = {3, 3, 3, 3};
94       PetscInt  dofAP2Hex[]  = {1, 1, 1, 1};
95       PetscInt  dofS3D[]     = {0, 0, 0, 6};
96       PetscInt *dofU, *dofA, *dofS;
97 
98       switch (sdim) {
99       case 2:
100         dofS = dofS2D;
101         break;
102       case 3:
103         dofS = dofS3D;
104         break;
105       default:
106         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No layout for dimension %" PetscInt_FMT, sdim);
107       }
108 
109       /* Identify cell type based on closure size only. This works for Tri/Tet/Quad/Hex meshes
110          It will not be enough to identify more exotic elements like pyramid or prisms...  */
111       PetscCall(ISGetIndices(cellIS, &cellID));
112       PetscCall(DMPlexGetTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));
113       switch (closureSize) {
114       case 7: /* Tri */
115         if (order == 1) {
116           dofU = dofUP1Tri;
117           dofA = dofAP1Tri;
118         } else {
119           dofU = dofUP2Tri;
120           dofA = dofAP2Tri;
121         }
122         break;
123       case 9: /* Quad */
124         if (order == 1) {
125           dofU = dofUP1Quad;
126           dofA = dofAP1Quad;
127         } else {
128           dofU = dofUP2Quad;
129           dofA = dofAP2Quad;
130         }
131         break;
132       case 15: /* Tet */
133         if (order == 1) {
134           dofU = dofUP1Tet;
135           dofA = dofAP1Tet;
136         } else {
137           dofU = dofUP2Tet;
138           dofA = dofAP2Tet;
139         }
140         break;
141       case 27: /* Hex */
142         if (order == 1) {
143           dofU = dofUP1Hex;
144           dofA = dofAP1Hex;
145         } else {
146           dofU = dofUP2Hex;
147           dofA = dofAP2Hex;
148         }
149         break;
150       default:
151         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Unknown element with closure size %" PetscInt_FMT, closureSize);
152       }
153       PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA));
154 
155       for (cell = 0; cell < numCells; cell++) {
156         PetscInt *closure = NULL;
157 
158         PetscCall(DMPlexGetTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
159         for (p = 0; p < closureSize; ++p) {
160           /* Find depth of p */
161           for (d = 0; d <= sdim; ++d) {
162             if ((closure[2 * p] >= pStartDepth[d]) && (closure[2 * p] < pEndDepth[d])) {
163               PetscCall(PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d]));
164               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d]));
165               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d]));
166               PetscCall(PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d]));
167             }
168           }
169         }
170         PetscCall(DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure));
171       }
172       PetscCall(ISRestoreIndices(cellIS, &cellID));
173       PetscCall(ISDestroy(&cellIS));
174     }
175   }
176   if (csIS) PetscCall(ISRestoreIndices(csIS, &csID));
177   PetscCall(ISDestroy(&csIS));
178   PetscCall(PetscSectionSetUp(section));
179   PetscCall(DMSetLocalSection(dm, section));
180   PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view"));
181   PetscCall(PetscSectionDestroy(&section));
182 
183   /* Get DM and IS for each field of dm */
184   PetscCall(DMCreateSubDM(dm, 1, &fieldU, &seqisU, &seqdmU));
185   PetscCall(DMCreateSubDM(dm, 1, &fieldA, &seqisA, &seqdmA));
186   PetscCall(DMCreateSubDM(dm, 1, &fieldS, &seqisS, &seqdmS));
187   PetscCall(DMCreateSubDM(dm, 2, fieldUA, &seqisUA, &seqdmUA));
188 
189   PetscCall(PetscMalloc1(2, &seqdmList));
190   seqdmList[0] = seqdmU;
191   seqdmList[1] = seqdmA;
192 
193   PetscCall(DMCreateSuperDM(seqdmList, 2, NULL, &seqdmUA2));
194   PetscCall(PetscFree(seqdmList));
195 
196   PetscCall(DMGetGlobalVector(dm, &seqX));
197   PetscCall(DMGetGlobalVector(seqdmU, &seqU));
198   PetscCall(DMGetGlobalVector(seqdmA, &seqA));
199   PetscCall(DMGetGlobalVector(seqdmS, &seqS));
200   PetscCall(DMGetGlobalVector(seqdmUA, &seqUA));
201   PetscCall(DMGetGlobalVector(seqdmUA2, &seqUA2));
202 
203   PetscCall(PetscObjectSetName((PetscObject)seqX, "seqX"));
204   PetscCall(PetscObjectSetName((PetscObject)seqU, "seqU"));
205   PetscCall(PetscObjectSetName((PetscObject)seqA, "seqAlpha"));
206   PetscCall(PetscObjectSetName((PetscObject)seqS, "seqSigma"));
207   PetscCall(PetscObjectSetName((PetscObject)seqUA, "seqUAlpha"));
208   PetscCall(PetscObjectSetName((PetscObject)seqUA2, "seqUAlpha2"));
209   PetscCall(VecSet(seqX, -111.));
210 
211   /* Setting u to [x,y,z]  and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */
212   {
213     PetscSection sectionUA;
214     Vec          UALoc;
215     PetscSection coordSection;
216     Vec          coord;
217     PetscScalar *cval, *xyz;
218     PetscInt     clSize, i, j;
219 
220     PetscCall(DMGetLocalSection(seqdmUA, &sectionUA));
221     PetscCall(DMGetLocalVector(seqdmUA, &UALoc));
222     PetscCall(VecGetArray(UALoc, &cval));
223     PetscCall(DMGetCoordinateSection(seqdmUA, &coordSection));
224     PetscCall(DMGetCoordinatesLocal(seqdmUA, &coord));
225     PetscCall(DMPlexGetChart(seqdmUA, &pStart, &pEnd));
226 
227     for (p = pStart; p < pEnd; ++p) {
228       PetscInt dofUA, offUA;
229 
230       PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA));
231       if (dofUA > 0) {
232         xyz = NULL;
233         PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA));
234         PetscCall(DMPlexVecGetClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz));
235         cval[offUA + sdim] = 0.;
236         for (i = 0; i < sdim; ++i) {
237           cval[offUA + i] = 0;
238           for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i];
239           cval[offUA + i] = cval[offUA + i] * sdim / clSize;
240           cval[offUA + sdim] += PetscSqr(cval[offUA + i]);
241         }
242         PetscCall(DMPlexVecRestoreClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz));
243       }
244     }
245     PetscCall(VecRestoreArray(UALoc, &cval));
246     PetscCall(DMLocalToGlobalBegin(seqdmUA, UALoc, INSERT_VALUES, seqUA));
247     PetscCall(DMLocalToGlobalEnd(seqdmUA, UALoc, INSERT_VALUES, seqUA));
248     PetscCall(DMRestoreLocalVector(seqdmUA, &UALoc));
249 
250     /* Update X */
251     PetscCall(VecISCopy(seqX, seqisUA, SCATTER_FORWARD, seqUA));
252     PetscCall(VecViewFromOptions(seqUA, NULL, "-ua_vec_view"));
253 
254     /* Restrict to U and Alpha */
255     PetscCall(VecISCopy(seqX, seqisU, SCATTER_REVERSE, seqU));
256     PetscCall(VecISCopy(seqX, seqisA, SCATTER_REVERSE, seqA));
257 
258     /* restrict to UA2 */
259     PetscCall(VecISCopy(seqX, seqisUA, SCATTER_REVERSE, seqUA2));
260     PetscCall(VecViewFromOptions(seqUA2, NULL, "-ua2_vec_view"));
261   }
262 
263   {
264     PetscInt         ovlp = 0;
265     PetscPartitioner part;
266 
267     PetscCall(DMSetUseNatural(dm, PETSC_TRUE));
268     PetscCall(DMPlexGetPartitioner(dm, &part));
269     PetscCall(PetscPartitionerSetFromOptions(part));
270     PetscCall(DMPlexDistribute(dm, ovlp, &migrationSF, &pdm));
271     if (!pdm) pdm = dm;
272     if (migrationSF) {
273       PetscCall(DMPlexSetMigrationSF(pdm, migrationSF));
274       PetscCall(PetscSFDestroy(&migrationSF));
275     }
276     PetscCall(DMViewFromOptions(pdm, NULL, "-dm_view"));
277   }
278 
279   /* Get DM and IS for each field of dm */
280   PetscCall(DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU));
281   PetscCall(DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA));
282   PetscCall(DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS));
283   PetscCall(DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA));
284 
285   PetscCall(PetscMalloc1(2, &dmList));
286   dmList[0] = dmU;
287   dmList[1] = dmA;
288 
289   PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2));
290   PetscCall(PetscFree(dmList));
291 
292   PetscCall(DMGetGlobalVector(pdm, &X));
293   PetscCall(DMGetGlobalVector(dmU, &U));
294   PetscCall(DMGetGlobalVector(dmA, &A));
295   PetscCall(DMGetGlobalVector(dmS, &S));
296   PetscCall(DMGetGlobalVector(dmUA, &UA));
297   PetscCall(DMGetGlobalVector(dmUA2, &UA2));
298 
299   PetscCall(PetscObjectSetName((PetscObject)X, "X"));
300   PetscCall(PetscObjectSetName((PetscObject)U, "U"));
301   PetscCall(PetscObjectSetName((PetscObject)A, "Alpha"));
302   PetscCall(PetscObjectSetName((PetscObject)S, "Sigma"));
303   PetscCall(PetscObjectSetName((PetscObject)UA, "UAlpha"));
304   PetscCall(PetscObjectSetName((PetscObject)UA2, "UAlpha2"));
305   PetscCall(VecSet(X, -111.));
306 
307   /* Setting u to [x,y,z]  and alpha to x^2+y^2+z^2 by writing in UAlpha then restricting to U and Alpha */
308   {
309     PetscSection sectionUA;
310     Vec          UALoc;
311     PetscSection coordSection;
312     Vec          coord;
313     PetscScalar *cval, *xyz;
314     PetscInt     clSize, i, j;
315 
316     PetscCall(DMGetLocalSection(dmUA, &sectionUA));
317     PetscCall(DMGetLocalVector(dmUA, &UALoc));
318     PetscCall(VecGetArray(UALoc, &cval));
319     PetscCall(DMGetCoordinateSection(dmUA, &coordSection));
320     PetscCall(DMGetCoordinatesLocal(dmUA, &coord));
321     PetscCall(DMPlexGetChart(dmUA, &pStart, &pEnd));
322 
323     for (p = pStart; p < pEnd; ++p) {
324       PetscInt dofUA, offUA;
325 
326       PetscCall(PetscSectionGetDof(sectionUA, p, &dofUA));
327       if (dofUA > 0) {
328         xyz = NULL;
329         PetscCall(PetscSectionGetOffset(sectionUA, p, &offUA));
330         PetscCall(DMPlexVecGetClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
331         cval[offUA + sdim] = 0.;
332         for (i = 0; i < sdim; ++i) {
333           cval[offUA + i] = 0;
334           for (j = 0; j < clSize / sdim; ++j) cval[offUA + i] += xyz[j * sdim + i];
335           cval[offUA + i] = cval[offUA + i] * sdim / clSize;
336           cval[offUA + sdim] += PetscSqr(cval[offUA + i]);
337         }
338         PetscCall(DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz));
339       }
340     }
341     PetscCall(VecRestoreArray(UALoc, &cval));
342     PetscCall(DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA));
343     PetscCall(DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA));
344     PetscCall(DMRestoreLocalVector(dmUA, &UALoc));
345 
346     /* Update X */
347     PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA));
348     PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view"));
349 
350     /* Restrict to U and Alpha */
351     PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U));
352     PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A));
353 
354     /* restrict to UA2 */
355     PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2));
356     PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view"));
357   }
358 
359   /* Verification */
360 
361   Vec       Xnat, Unat, Anat, UAnat, Snat, UA2nat;
362   PetscReal norm;
363   PetscCall(DMGetGlobalVector(dm, &Xnat));
364   PetscCall(DMGetGlobalVector(seqdmU, &Unat));
365   PetscCall(DMGetGlobalVector(seqdmA, &Anat));
366   PetscCall(DMGetGlobalVector(seqdmUA, &UAnat));
367   PetscCall(DMGetGlobalVector(seqdmS, &Snat));
368   PetscCall(DMGetGlobalVector(seqdmUA2, &UA2nat));
369 
370   PetscCall(DMPlexGlobalToNaturalBegin(pdm, X, Xnat));
371   PetscCall(DMPlexGlobalToNaturalEnd(pdm, X, Xnat));
372   PetscCall(VecAXPY(Xnat, -1.0, seqX));
373   PetscCall(VecNorm(Xnat, NORM_INFINITY, &norm));
374   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "X ||Vin - Vout|| = %g", (double)norm);
375 
376   PetscCall(DMPlexGlobalToNaturalBegin(dmU, U, Unat));
377   PetscCall(DMPlexGlobalToNaturalEnd(dmU, U, Unat));
378   PetscCall(VecAXPY(Unat, -1.0, seqU));
379   PetscCall(VecNorm(Unat, NORM_INFINITY, &norm));
380   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "U ||Vin - Vout|| = %g", (double)norm);
381 
382   PetscCall(DMPlexGlobalToNaturalBegin(dmA, A, Anat));
383   PetscCall(DMPlexGlobalToNaturalEnd(dmA, A, Anat));
384   PetscCall(VecAXPY(Anat, -1.0, seqA));
385   PetscCall(VecNorm(Anat, NORM_INFINITY, &norm));
386   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "A ||Vin - Vout|| = %g", (double)norm);
387 
388   PetscCall(DMPlexGlobalToNaturalBegin(dmUA, UA, UAnat));
389   PetscCall(DMPlexGlobalToNaturalEnd(dmUA, UA, UAnat));
390   PetscCall(VecAXPY(UAnat, -1.0, seqUA));
391   PetscCall(VecNorm(UAnat, NORM_INFINITY, &norm));
392   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UA ||Vin - Vout|| = %g", (double)norm);
393 
394   PetscCall(DMPlexGlobalToNaturalBegin(dmS, S, Snat));
395   PetscCall(DMPlexGlobalToNaturalEnd(dmS, S, Snat));
396   PetscCall(VecAXPY(Snat, -1.0, seqS));
397   PetscCall(VecNorm(Snat, NORM_INFINITY, &norm));
398   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "S ||Vin - Vout|| = %g", (double)norm);
399 
400   PetscCall(DMPlexGlobalToNaturalBegin(dmUA2, UA2, UA2nat));
401   PetscCall(DMPlexGlobalToNaturalEnd(dmUA2, UA2, UA2nat));
402   PetscCall(VecAXPY(UA2nat, -1.0, seqUA2));
403   PetscCall(VecNorm(UA2nat, NORM_INFINITY, &norm));
404   PetscCheck(norm <= PETSC_SQRT_MACHINE_EPSILON, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "UAlpha2 ||Vin - Vout|| = %g", (double)norm);
405 
406   PetscCall(DMRestoreGlobalVector(dm, &Xnat));
407   PetscCall(DMRestoreGlobalVector(seqdmU, &Unat));
408   PetscCall(DMRestoreGlobalVector(seqdmA, &Anat));
409   PetscCall(DMRestoreGlobalVector(seqdmUA, &UAnat));
410   PetscCall(DMRestoreGlobalVector(seqdmS, &Snat));
411   PetscCall(DMRestoreGlobalVector(seqdmUA2, &UA2nat));
412 
413   PetscCall(DMRestoreGlobalVector(seqdmUA2, &seqUA2));
414   PetscCall(DMRestoreGlobalVector(seqdmUA, &seqUA));
415   PetscCall(DMRestoreGlobalVector(seqdmS, &seqS));
416   PetscCall(DMRestoreGlobalVector(seqdmA, &seqA));
417   PetscCall(DMRestoreGlobalVector(seqdmU, &seqU));
418   PetscCall(DMRestoreGlobalVector(dm, &seqX));
419   PetscCall(DMDestroy(&seqdmU));
420   PetscCall(ISDestroy(&seqisU));
421   PetscCall(DMDestroy(&seqdmA));
422   PetscCall(ISDestroy(&seqisA));
423   PetscCall(DMDestroy(&seqdmS));
424   PetscCall(ISDestroy(&seqisS));
425   PetscCall(DMDestroy(&seqdmUA));
426   PetscCall(ISDestroy(&seqisUA));
427   PetscCall(DMDestroy(&seqdmUA2));
428 
429   PetscCall(DMRestoreGlobalVector(dmUA2, &UA2));
430   PetscCall(DMRestoreGlobalVector(dmUA, &UA));
431   PetscCall(DMRestoreGlobalVector(dmS, &S));
432   PetscCall(DMRestoreGlobalVector(dmA, &A));
433   PetscCall(DMRestoreGlobalVector(dmU, &U));
434   PetscCall(DMRestoreGlobalVector(pdm, &X));
435   PetscCall(DMDestroy(&dmU));
436   PetscCall(ISDestroy(&isU));
437   PetscCall(DMDestroy(&dmA));
438   PetscCall(ISDestroy(&isA));
439   PetscCall(DMDestroy(&dmS));
440   PetscCall(ISDestroy(&isS));
441   PetscCall(DMDestroy(&dmUA));
442   PetscCall(ISDestroy(&isUA));
443   PetscCall(DMDestroy(&dmUA2));
444   PetscCall(DMDestroy(&pdm));
445   if (size > 1) PetscCall(DMDestroy(&dm));
446   PetscCall(PetscFree2(pStartDepth, pEndDepth));
447   PetscCall(PetscFinalize());
448   return 0;
449 }
450 
451 /*TEST
452 
453   build:
454     requires: !complex parmetis exodusii pnetcdf
455   # 2D seq
456   test:
457     suffix: 0
458     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis
459   test:
460     suffix: 1
461     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis
462   test:
463     suffix: 2
464     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis
465 
466   # 2D par
467   test:
468     suffix: 3
469     nsize: 2
470     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareT-large.exo -dm_view -petscpartitioner_type parmetis
471   test:
472     suffix: 4
473     nsize: 2
474     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareQ-large.exo -dm_view -petscpartitioner_type parmetis
475   test:
476     suffix: 5
477     nsize: 2
478     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourSquareH-large.exo -dm_view -petscpartitioner_type parmetis
479 
480   #3d seq
481   test:
482     suffix: 6
483     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis
484   test:
485     suffix: 7
486     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis
487 
488   #3d par
489   test:
490     suffix: 8
491     nsize: 2
492     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickHex-large.exo -dm_view -petscpartitioner_type parmetis
493   test:
494     suffix: 9
495     nsize: 2
496     args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/FourBrickTet-large.exo -dm_view -petscpartitioner_type parmetis
497 
498 TEST*/
499