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: PetscInitialize(&argc, &argv, NULL, help);
37: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
38: MPI_Comm_size(PETSC_COMM_WORLD, &size);
39: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex64");
40: 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: DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, "ex64_plex", PETSC_TRUE, &dm);
45: DMPlexDistributeSetDefault(dm, PETSC_FALSE);
46: DMSetFromOptions(dm);
47: DMViewFromOptions(dm, NULL, "-dm_view");
48: DMGetDimension(dm, &sdim);
50: /* Create the main section containing all fields */
51: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
52: PetscSectionSetNumFields(section, 3);
53: PetscSectionSetFieldName(section, fieldU, "U");
54: PetscSectionSetFieldName(section, fieldA, "Alpha");
55: PetscSectionSetFieldName(section, fieldS, "Sigma");
56: DMPlexGetChart(dm, &pStart, &pEnd);
57: PetscSectionSetChart(section, pStart, pEnd);
58: PetscMalloc2(sdim + 1, &pStartDepth, sdim + 1, &pEndDepth);
59: for (d = 0; d <= sdim; ++d) DMPlexGetDepthStratum(dm, d, &pStartDepth[d], &pEndDepth[d]);
60: /* Vector field U, Scalar field Alpha, Tensor field Sigma */
61: PetscSectionSetFieldComponents(section, fieldU, sdim);
62: PetscSectionSetFieldComponents(section, fieldA, 1);
63: PetscSectionSetFieldComponents(section, fieldS, sdim * (sdim + 1) / 2);
65: /* Going through cell sets then cells, and setting up storage for the sections */
66: DMGetLabelSize(dm, "Cell Sets", &numCS);
67: DMGetLabelIdIS(dm, "Cell Sets", &csIS);
68: if (csIS) 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: DMGetStratumSize(dm, "Cell Sets", csID[set], &numCells);
75: 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: ISGetIndices(cellIS, &cellID);
112: 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: DMPlexRestoreTransitiveClosure(dm, cellID[0], PETSC_TRUE, &closureSize, &closureA);
155: for (cell = 0; cell < numCells; cell++) {
156: PetscInt *closure = NULL;
158: 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: PetscSectionSetDof(section, closure[2 * p], dofU[d] + dofA[d] + dofS[d]);
164: PetscSectionSetFieldDof(section, closure[2 * p], fieldU, dofU[d]);
165: PetscSectionSetFieldDof(section, closure[2 * p], fieldA, dofA[d]);
166: PetscSectionSetFieldDof(section, closure[2 * p], fieldS, dofS[d]);
167: }
168: }
169: }
170: DMPlexRestoreTransitiveClosure(dm, cellID[cell], PETSC_TRUE, &closureSize, &closure);
171: }
172: ISRestoreIndices(cellIS, &cellID);
173: ISDestroy(&cellIS);
174: }
175: }
176: if (csIS) ISRestoreIndices(csIS, &csID);
177: ISDestroy(&csIS);
178: PetscSectionSetUp(section);
179: DMSetLocalSection(dm, section);
180: PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view");
181: PetscSectionDestroy(§ion);
183: /* Get DM and IS for each field of dm */
184: DMCreateSubDM(dm, 1, &fieldU, &seqisU, &seqdmU);
185: DMCreateSubDM(dm, 1, &fieldA, &seqisA, &seqdmA);
186: DMCreateSubDM(dm, 1, &fieldS, &seqisS, &seqdmS);
187: DMCreateSubDM(dm, 2, fieldUA, &seqisUA, &seqdmUA);
189: PetscMalloc1(2, &seqdmList);
190: seqdmList[0] = seqdmU;
191: seqdmList[1] = seqdmA;
193: DMCreateSuperDM(seqdmList, 2, NULL, &seqdmUA2);
194: PetscFree(seqdmList);
196: DMGetGlobalVector(dm, &seqX);
197: DMGetGlobalVector(seqdmU, &seqU);
198: DMGetGlobalVector(seqdmA, &seqA);
199: DMGetGlobalVector(seqdmS, &seqS);
200: DMGetGlobalVector(seqdmUA, &seqUA);
201: DMGetGlobalVector(seqdmUA2, &seqUA2);
203: PetscObjectSetName((PetscObject)seqX, "seqX");
204: PetscObjectSetName((PetscObject)seqU, "seqU");
205: PetscObjectSetName((PetscObject)seqA, "seqAlpha");
206: PetscObjectSetName((PetscObject)seqS, "seqSigma");
207: PetscObjectSetName((PetscObject)seqUA, "seqUAlpha");
208: PetscObjectSetName((PetscObject)seqUA2, "seqUAlpha2");
209: 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: DMGetLocalSection(seqdmUA, §ionUA);
221: DMGetLocalVector(seqdmUA, &UALoc);
222: VecGetArray(UALoc, &cval);
223: DMGetCoordinateSection(seqdmUA, &coordSection);
224: DMGetCoordinatesLocal(seqdmUA, &coord);
225: DMPlexGetChart(seqdmUA, &pStart, &pEnd);
227: for (p = pStart; p < pEnd; ++p) {
228: PetscInt dofUA, offUA;
230: PetscSectionGetDof(sectionUA, p, &dofUA);
231: if (dofUA > 0) {
232: xyz = NULL;
233: PetscSectionGetOffset(sectionUA, p, &offUA);
234: 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: DMPlexVecRestoreClosure(seqdmUA, coordSection, coord, p, &clSize, &xyz);
243: }
244: }
245: VecRestoreArray(UALoc, &cval);
246: DMLocalToGlobalBegin(seqdmUA, UALoc, INSERT_VALUES, seqUA);
247: DMLocalToGlobalEnd(seqdmUA, UALoc, INSERT_VALUES, seqUA);
248: DMRestoreLocalVector(seqdmUA, &UALoc);
250: /* Update X */
251: VecISCopy(seqX, seqisUA, SCATTER_FORWARD, seqUA);
252: VecViewFromOptions(seqUA, NULL, "-ua_vec_view");
254: /* Restrict to U and Alpha */
255: VecISCopy(seqX, seqisU, SCATTER_REVERSE, seqU);
256: VecISCopy(seqX, seqisA, SCATTER_REVERSE, seqA);
258: /* restrict to UA2 */
259: VecISCopy(seqX, seqisUA, SCATTER_REVERSE, seqUA2);
260: VecViewFromOptions(seqUA2, NULL, "-ua2_vec_view");
261: }
263: {
264: PetscInt ovlp = 0;
265: PetscPartitioner part;
267: DMSetUseNatural(dm, PETSC_TRUE);
268: DMPlexGetPartitioner(dm, &part);
269: PetscPartitionerSetFromOptions(part);
270: DMPlexDistribute(dm, ovlp, &migrationSF, &pdm);
271: if (!pdm) pdm = dm;
272: if (migrationSF) {
273: DMPlexSetMigrationSF(pdm, migrationSF);
274: PetscSFDestroy(&migrationSF);
275: }
276: DMViewFromOptions(pdm, NULL, "-dm_view");
277: }
279: /* Get DM and IS for each field of dm */
280: DMCreateSubDM(pdm, 1, &fieldU, &isU, &dmU);
281: DMCreateSubDM(pdm, 1, &fieldA, &isA, &dmA);
282: DMCreateSubDM(pdm, 1, &fieldS, &isS, &dmS);
283: DMCreateSubDM(pdm, 2, fieldUA, &isUA, &dmUA);
285: PetscMalloc1(2, &dmList);
286: dmList[0] = dmU;
287: dmList[1] = dmA;
289: DMCreateSuperDM(dmList, 2, NULL, &dmUA2);
290: PetscFree(dmList);
292: DMGetGlobalVector(pdm, &X);
293: DMGetGlobalVector(dmU, &U);
294: DMGetGlobalVector(dmA, &A);
295: DMGetGlobalVector(dmS, &S);
296: DMGetGlobalVector(dmUA, &UA);
297: DMGetGlobalVector(dmUA2, &UA2);
299: PetscObjectSetName((PetscObject)X, "X");
300: PetscObjectSetName((PetscObject)U, "U");
301: PetscObjectSetName((PetscObject)A, "Alpha");
302: PetscObjectSetName((PetscObject)S, "Sigma");
303: PetscObjectSetName((PetscObject)UA, "UAlpha");
304: PetscObjectSetName((PetscObject)UA2, "UAlpha2");
305: 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: DMGetLocalSection(dmUA, §ionUA);
317: DMGetLocalVector(dmUA, &UALoc);
318: VecGetArray(UALoc, &cval);
319: DMGetCoordinateSection(dmUA, &coordSection);
320: DMGetCoordinatesLocal(dmUA, &coord);
321: DMPlexGetChart(dmUA, &pStart, &pEnd);
323: for (p = pStart; p < pEnd; ++p) {
324: PetscInt dofUA, offUA;
326: PetscSectionGetDof(sectionUA, p, &dofUA);
327: if (dofUA > 0) {
328: xyz = NULL;
329: PetscSectionGetOffset(sectionUA, p, &offUA);
330: 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: DMPlexVecRestoreClosure(dmUA, coordSection, coord, p, &clSize, &xyz);
339: }
340: }
341: VecRestoreArray(UALoc, &cval);
342: DMLocalToGlobalBegin(dmUA, UALoc, INSERT_VALUES, UA);
343: DMLocalToGlobalEnd(dmUA, UALoc, INSERT_VALUES, UA);
344: DMRestoreLocalVector(dmUA, &UALoc);
346: /* Update X */
347: VecISCopy(X, isUA, SCATTER_FORWARD, UA);
348: VecViewFromOptions(UA, NULL, "-ua_vec_view");
350: /* Restrict to U and Alpha */
351: VecISCopy(X, isU, SCATTER_REVERSE, U);
352: VecISCopy(X, isA, SCATTER_REVERSE, A);
354: /* restrict to UA2 */
355: VecISCopy(X, isUA, SCATTER_REVERSE, UA2);
356: VecViewFromOptions(UA2, NULL, "-ua2_vec_view");
357: }
359: /* Verification */
361: Vec Xnat, Unat, Anat, UAnat, Snat, UA2nat;
362: PetscReal norm;
363: DMGetGlobalVector(dm, &Xnat);
364: DMGetGlobalVector(seqdmU, &Unat);
365: DMGetGlobalVector(seqdmA, &Anat);
366: DMGetGlobalVector(seqdmUA, &UAnat);
367: DMGetGlobalVector(seqdmS, &Snat);
368: DMGetGlobalVector(seqdmUA2, &UA2nat);
370: DMPlexGlobalToNaturalBegin(pdm, X, Xnat);
371: DMPlexGlobalToNaturalEnd(pdm, X, Xnat);
372: VecAXPY(Xnat, -1.0, seqX);
373: VecNorm(Xnat, NORM_INFINITY, &norm);
376: DMPlexGlobalToNaturalBegin(dmU, U, Unat);
377: DMPlexGlobalToNaturalEnd(dmU, U, Unat);
378: VecAXPY(Unat, -1.0, seqU);
379: VecNorm(Unat, NORM_INFINITY, &norm);
382: DMPlexGlobalToNaturalBegin(dmA, A, Anat);
383: DMPlexGlobalToNaturalEnd(dmA, A, Anat);
384: VecAXPY(Anat, -1.0, seqA);
385: VecNorm(Anat, NORM_INFINITY, &norm);
388: DMPlexGlobalToNaturalBegin(dmUA, UA, UAnat);
389: DMPlexGlobalToNaturalEnd(dmUA, UA, UAnat);
390: VecAXPY(UAnat, -1.0, seqUA);
391: VecNorm(UAnat, NORM_INFINITY, &norm);
394: DMPlexGlobalToNaturalBegin(dmS, S, Snat);
395: DMPlexGlobalToNaturalEnd(dmS, S, Snat);
396: VecAXPY(Snat, -1.0, seqS);
397: VecNorm(Snat, NORM_INFINITY, &norm);
400: DMPlexGlobalToNaturalBegin(dmUA2, UA2, UA2nat);
401: DMPlexGlobalToNaturalEnd(dmUA2, UA2, UA2nat);
402: VecAXPY(UA2nat, -1.0, seqUA2);
403: VecNorm(UA2nat, NORM_INFINITY, &norm);
406: DMRestoreGlobalVector(dm, &Xnat);
407: DMRestoreGlobalVector(seqdmU, &Unat);
408: DMRestoreGlobalVector(seqdmA, &Anat);
409: DMRestoreGlobalVector(seqdmUA, &UAnat);
410: DMRestoreGlobalVector(seqdmS, &Snat);
411: DMRestoreGlobalVector(seqdmUA2, &UA2nat);
413: DMRestoreGlobalVector(seqdmUA2, &seqUA2);
414: DMRestoreGlobalVector(seqdmUA, &seqUA);
415: DMRestoreGlobalVector(seqdmS, &seqS);
416: DMRestoreGlobalVector(seqdmA, &seqA);
417: DMRestoreGlobalVector(seqdmU, &seqU);
418: DMRestoreGlobalVector(dm, &seqX);
419: DMDestroy(&seqdmU);
420: ISDestroy(&seqisU);
421: DMDestroy(&seqdmA);
422: ISDestroy(&seqisA);
423: DMDestroy(&seqdmS);
424: ISDestroy(&seqisS);
425: DMDestroy(&seqdmUA);
426: ISDestroy(&seqisUA);
427: DMDestroy(&seqdmUA2);
429: DMRestoreGlobalVector(dmUA2, &UA2);
430: DMRestoreGlobalVector(dmUA, &UA);
431: DMRestoreGlobalVector(dmS, &S);
432: DMRestoreGlobalVector(dmA, &A);
433: DMRestoreGlobalVector(dmU, &U);
434: DMRestoreGlobalVector(pdm, &X);
435: DMDestroy(&dmU);
436: ISDestroy(&isU);
437: DMDestroy(&dmA);
438: ISDestroy(&isA);
439: DMDestroy(&dmS);
440: ISDestroy(&isS);
441: DMDestroy(&dmUA);
442: ISDestroy(&isUA);
443: DMDestroy(&dmUA2);
444: DMDestroy(&pdm);
445: if (size > 1) DMDestroy(&dm);
446: PetscFree2(pStartDepth, pEndDepth);
447: 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*/