Actual source code: ex64.c

  1: static char help[] = "Test FEM layout and GlobalToNaturalSF\n\n";

  3: /*
  4:   In order to see the vectors which are being tested, use

  6:      -ua_vec_view -s_vec_view
  7: */

  9: #include <petsc.h>
 10: #include <exodusII.h>

 12: #include <petsc/private/dmpleximpl.h>

 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;

 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();

 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));

 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));

 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;

 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;

 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:       }

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));

155:       for (cell = 0; cell < numCells; cell++) {
156:         PetscInt *closure = NULL;

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));

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));

189:   PetscCall(PetscMalloc1(2, &seqdmList));
190:   seqdmList[0] = seqdmU;
191:   seqdmList[1] = seqdmA;

193:   PetscCall(DMCreateSuperDM(seqdmList, 2, NULL, &seqdmUA2));
194:   PetscCall(PetscFree(seqdmList));

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));

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.));

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;

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));

227:     for (p = pStart; p < pEnd; ++p) {
228:       PetscInt dofUA, offUA;

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));

250:     /* Update X */
251:     PetscCall(VecISCopy(seqX, seqisUA, SCATTER_FORWARD, seqUA));
252:     PetscCall(VecViewFromOptions(seqUA, NULL, "-ua_vec_view"));

254:     /* Restrict to U and Alpha */
255:     PetscCall(VecISCopy(seqX, seqisU, SCATTER_REVERSE, seqU));
256:     PetscCall(VecISCopy(seqX, seqisA, SCATTER_REVERSE, seqA));

258:     /* restrict to UA2 */
259:     PetscCall(VecISCopy(seqX, seqisUA, SCATTER_REVERSE, seqUA2));
260:     PetscCall(VecViewFromOptions(seqUA2, NULL, "-ua2_vec_view"));
261:   }

263:   {
264:     PetscInt         ovlp = 0;
265:     PetscPartitioner part;

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:   }

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));

285:   PetscCall(PetscMalloc1(2, &dmList));
286:   dmList[0] = dmU;
287:   dmList[1] = dmA;

289:   PetscCall(DMCreateSuperDM(dmList, 2, NULL, &dmUA2));
290:   PetscCall(PetscFree(dmList));

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));

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.));

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;

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));

323:     for (p = pStart; p < pEnd; ++p) {
324:       PetscInt dofUA, offUA;

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));

346:     /* Update X */
347:     PetscCall(VecISCopy(X, isUA, SCATTER_FORWARD, UA));
348:     PetscCall(VecViewFromOptions(UA, NULL, "-ua_vec_view"));

350:     /* Restrict to U and Alpha */
351:     PetscCall(VecISCopy(X, isU, SCATTER_REVERSE, U));
352:     PetscCall(VecISCopy(X, isA, SCATTER_REVERSE, A));

354:     /* restrict to UA2 */
355:     PetscCall(VecISCopy(X, isUA, SCATTER_REVERSE, UA2));
356:     PetscCall(VecViewFromOptions(UA2, NULL, "-ua2_vec_view"));
357:   }

359:   /* Verification */

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));

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);

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);

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);

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);

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);

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);

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));

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));

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: }

451: /*TEST

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

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

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

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

498: TEST*/