Actual source code: ex1.c
2: #include <petscfe.h>
3: #include <petscdmplex.h>
4: #include <petsc/private/hashmap.h>
5: #include <petsc/private/dmpleximpl.h>
6: #include <petsc/private/petscfeimpl.h>
8: const char help[] = "Test PETSCDUALSPACELAGRANGE\n";
10: typedef struct _PetscHashLagKey {
11: PetscInt dim;
12: PetscInt order;
13: PetscInt formDegree;
14: PetscBool trimmed;
15: PetscInt tensor;
16: PetscBool continuous;
17: } PetscHashLagKey;
19: #define PetscHashLagKeyHash(key) \
20: PetscHashCombine(PetscHashCombine(PetscHashCombine(PetscHashInt((key).dim), PetscHashInt((key).order)), PetscHashInt((key).formDegree)), PetscHashCombine(PetscHashCombine(PetscHashInt((key).trimmed), PetscHashInt((key).tensor)), PetscHashInt((key).continuous)))
22: #define PetscHashLagKeyEqual(k1, k2) \
23: (((k1).dim == (k2).dim) ? ((k1).order == (k2).order) ? ((k1).formDegree == (k2).formDegree) ? ((k1).trimmed == (k2).trimmed) ? ((k1).tensor == (k2).tensor) ? ((k1).continuous == (k2).continuous) : 0 : 0 : 0 : 0 : 0)
25: PETSC_HASH_MAP(HashLag, PetscHashLagKey, PetscInt, PetscHashLagKeyHash, PetscHashLagKeyEqual, 0)
27: static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs);
28: static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs);
30: static PetscErrorCode ExpectedNumDofs_Total(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs)
31: {
32: formDegree = PetscAbsInt(formDegree);
33: /* see femtable.org for the source of most of these values */
34: *nDofs = -1;
35: if (tensor == 0) { /* simplex spaces */
36: if (trimmed) {
37: PetscInt rnchooserk;
38: PetscInt rkm1choosek;
40: PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk);
41: PetscDTBinomialInt(order + formDegree - 1, formDegree, &rkm1choosek);
42: *nDofs = rnchooserk * rkm1choosek * nCopies;
43: } else {
44: PetscInt rnchooserk;
45: PetscInt rkchoosek;
47: PetscDTBinomialInt(order + dim, order + formDegree, &rnchooserk);
48: PetscDTBinomialInt(order + formDegree, formDegree, &rkchoosek);
49: *nDofs = rnchooserk * rkchoosek * nCopies;
50: }
51: } else if (tensor == 1) { /* hypercubes */
52: if (trimmed) {
53: PetscInt nchoosek;
54: PetscInt rpowk, rp1pownmk;
56: PetscDTBinomialInt(dim, formDegree, &nchoosek);
57: rpowk = PetscPowInt(order, formDegree);
58: rp1pownmk = PetscPowInt(order + 1, dim - formDegree);
59: *nDofs = nchoosek * rpowk * rp1pownmk * nCopies;
60: } else {
61: PetscInt nchoosek;
62: PetscInt rp1pown;
64: PetscDTBinomialInt(dim, formDegree, &nchoosek);
65: rp1pown = PetscPowInt(order + 1, dim);
66: *nDofs = nchoosek * rp1pown * nCopies;
67: }
68: } else { /* prism */
69: PetscInt tracek = 0;
70: PetscInt tracekm1 = 0;
71: PetscInt fiber0 = 0;
72: PetscInt fiber1 = 0;
74: if (formDegree < dim) {
75: ExpectedNumDofs_Total(dim - 1, order, formDegree, trimmed, 0, 1, &tracek);
76: ExpectedNumDofs_Total(1, order, 0, trimmed, 0, 1, &fiber0);
77: }
78: if (formDegree > 0) {
79: ExpectedNumDofs_Total(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1);
80: ExpectedNumDofs_Total(1, order, 1, trimmed, 0, 1, &fiber1);
81: }
82: *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies;
83: }
84: return 0;
85: }
87: static PetscErrorCode ExpectedNumDofs_Interior(PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscInt nCopies, PetscInt *nDofs)
88: {
89: formDegree = PetscAbsInt(formDegree);
90: /* see femtable.org for the source of most of these values */
91: *nDofs = -1;
92: if (tensor == 0) { /* simplex spaces */
93: if (trimmed) {
94: if (order + formDegree > dim) {
95: PetscInt eorder = order + formDegree - dim - 1;
96: PetscInt eformDegree = dim - formDegree;
97: PetscInt rnchooserk;
98: PetscInt rkchoosek;
100: PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk);
101: PetscDTBinomialInt(eorder + eformDegree, eformDegree, &rkchoosek);
102: *nDofs = rnchooserk * rkchoosek * nCopies;
103: } else {
104: *nDofs = 0;
105: }
107: } else {
108: if (order + formDegree > dim) {
109: PetscInt eorder = order + formDegree - dim;
110: PetscInt eformDegree = dim - formDegree;
111: PetscInt rnchooserk;
112: PetscInt rkm1choosek;
114: PetscDTBinomialInt(eorder + dim, eorder + eformDegree, &rnchooserk);
115: PetscDTBinomialInt(eorder + eformDegree - 1, eformDegree, &rkm1choosek);
116: *nDofs = rnchooserk * rkm1choosek * nCopies;
117: } else {
118: *nDofs = 0;
119: }
120: }
121: } else if (tensor == 1) { /* hypercubes */
122: if (dim < 2) {
123: ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, 0, nCopies, nDofs);
124: } else {
125: PetscInt tracek = 0;
126: PetscInt tracekm1 = 0;
127: PetscInt fiber0 = 0;
128: PetscInt fiber1 = 0;
130: if (formDegree < dim) {
131: ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, dim > 2 ? 1 : 0, 1, &tracek);
132: ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0);
133: }
134: if (formDegree > 0) {
135: ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, dim > 2 ? 1 : 0, 1, &tracekm1);
136: ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1);
137: }
138: *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies;
139: }
140: } else { /* prism */
141: PetscInt tracek = 0;
142: PetscInt tracekm1 = 0;
143: PetscInt fiber0 = 0;
144: PetscInt fiber1 = 0;
146: if (formDegree < dim) {
147: ExpectedNumDofs_Interior(dim - 1, order, formDegree, trimmed, 0, 1, &tracek);
148: ExpectedNumDofs_Interior(1, order, 0, trimmed, 0, 1, &fiber0);
149: }
150: if (formDegree > 0) {
151: ExpectedNumDofs_Interior(dim - 1, order, formDegree - 1, trimmed, 0, 1, &tracekm1);
152: ExpectedNumDofs_Interior(1, order, 1, trimmed, 0, 1, &fiber1);
153: }
154: *nDofs = (tracek * fiber0 + tracekm1 * fiber1) * nCopies;
155: }
156: return 0;
157: }
159: PetscErrorCode testLagrange(PetscHashLag lagTable, DM K, PetscInt dim, PetscInt order, PetscInt formDegree, PetscBool trimmed, PetscInt tensor, PetscBool continuous, PetscInt nCopies)
160: {
161: PetscDualSpace sp;
162: MPI_Comm comm = PETSC_COMM_SELF;
163: PetscInt Nk;
164: PetscHashLagKey key;
165: PetscHashIter iter;
166: PetscBool missing;
167: PetscInt spdim, spintdim, exspdim, exspintdim;
169: PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);
170: PetscDualSpaceCreate(comm, &sp);
171: PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE);
172: PetscDualSpaceSetDM(sp, K);
173: PetscDualSpaceSetOrder(sp, order);
174: PetscDualSpaceSetFormDegree(sp, formDegree);
175: PetscDualSpaceSetNumComponents(sp, nCopies * Nk);
176: PetscDualSpaceLagrangeSetContinuity(sp, continuous);
177: PetscDualSpaceLagrangeSetTensor(sp, (PetscBool)tensor);
178: PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);
179: PetscInfo(NULL, "Input: dim %" PetscInt_FMT ", order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensor %" PetscInt_FMT ", continuous %" PetscInt_FMT ", formDegree %" PetscInt_FMT ", nCopies %" PetscInt_FMT "\n", dim, order, (PetscInt)trimmed, tensor, (PetscInt)continuous, formDegree, nCopies);
180: ExpectedNumDofs_Total(dim, order, formDegree, trimmed, tensor, nCopies, &exspdim);
181: if (continuous && dim > 0 && order > 0) {
182: ExpectedNumDofs_Interior(dim, order, formDegree, trimmed, tensor, nCopies, &exspintdim);
183: } else {
184: exspintdim = exspdim;
185: }
186: PetscDualSpaceSetUp(sp);
187: PetscDualSpaceGetDimension(sp, &spdim);
188: PetscDualSpaceGetInteriorDimension(sp, &spintdim);
191: key.dim = dim;
192: key.formDegree = formDegree;
193: PetscDualSpaceGetOrder(sp, &key.order);
194: PetscDualSpaceLagrangeGetContinuity(sp, &key.continuous);
195: if (tensor == 2) {
196: key.tensor = 2;
197: } else {
198: PetscBool bTensor;
200: PetscDualSpaceLagrangeGetTensor(sp, &bTensor);
201: key.tensor = bTensor;
202: }
203: PetscDualSpaceLagrangeGetTrimmed(sp, &key.trimmed);
204: PetscInfo(NULL, "After setup: order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensor %" PetscInt_FMT ", continuous %" PetscInt_FMT "\n", key.order, (PetscInt)key.trimmed, key.tensor, (PetscInt)key.continuous);
205: PetscHashLagPut(lagTable, key, &iter, &missing);
206: if (missing) {
207: DMPolytopeType type;
209: DMPlexGetCellType(K, 0, &type);
210: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "New space: %s, order %" PetscInt_FMT ", trimmed %" PetscInt_FMT ", tensor %" PetscInt_FMT ", continuous %" PetscInt_FMT ", form degree %" PetscInt_FMT "\n", DMPolytopeTypes[type], order, (PetscInt)trimmed, tensor, (PetscInt)continuous, formDegree);
211: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);
212: {
213: PetscQuadrature intNodes, allNodes;
214: Mat intMat, allMat;
215: MatInfo info;
216: PetscInt i, j, nodeIdxDim, nodeVecDim, nNodes;
217: const PetscInt *nodeIdx;
218: const PetscReal *nodeVec;
220: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data;
222: PetscLagNodeIndicesGetData_Internal(lag->allNodeIndices, &nodeIdxDim, &nodeVecDim, &nNodes, &nodeIdx, &nodeVec);
226: PetscDualSpaceGetAllData(sp, &allNodes, &allMat);
228: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All nodes:\n");
229: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);
230: PetscQuadratureView(allNodes, PETSC_VIEWER_STDOUT_SELF);
231: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);
232: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All node indices:\n");
233: for (i = 0; i < spdim; i++) {
234: PetscPrintf(PETSC_COMM_SELF, "(");
235: for (j = 0; j < nodeIdxDim; j++) PetscPrintf(PETSC_COMM_SELF, " %" PetscInt_FMT ",", nodeIdx[i * nodeIdxDim + j]);
236: PetscPrintf(PETSC_COMM_SELF, "): [");
237: for (j = 0; j < nodeVecDim; j++) PetscPrintf(PETSC_COMM_SELF, " %g,", (double)nodeVec[i * nodeVecDim + j]);
238: PetscPrintf(PETSC_COMM_SELF, "]\n");
239: }
241: MatGetInfo(allMat, MAT_LOCAL, &info);
242: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "All matrix: %" PetscInt_FMT " nonzeros\n", (PetscInt)info.nz_used);
244: PetscDualSpaceGetInteriorData(sp, &intNodes, &intMat);
245: if (intMat && intMat != allMat) {
246: PetscInt intNodeIdxDim, intNodeVecDim, intNnodes;
247: const PetscInt *intNodeIdx;
248: const PetscReal *intNodeVec;
249: PetscBool same;
251: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior nodes:\n");
252: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);
253: PetscQuadratureView(intNodes, PETSC_VIEWER_STDOUT_SELF);
254: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);
256: MatGetInfo(intMat, MAT_LOCAL, &info);
257: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior matrix: %" PetscInt_FMT " nonzeros\n", (PetscInt)info.nz_used);
258: PetscLagNodeIndicesGetData_Internal(lag->intNodeIndices, &intNodeIdxDim, &intNodeVecDim, &intNnodes, &intNodeIdx, &intNodeVec);
261: PetscArraycmp(intNodeIdx, nodeIdx, nodeIdxDim * intNnodes, &same);
263: PetscArraycmp(intNodeVec, nodeVec, nodeVecDim * intNnodes, &same);
265: } else if (intMat) {
266: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior data is the same as all data\n");
269: }
270: }
271: if (dim <= 2 && spintdim) {
272: PetscInt numFaces, o;
274: {
275: DMPolytopeType ct;
276: /* The number of arrangements is no longer based on the number of faces */
277: DMPlexGetCellType(K, 0, &ct);
278: numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2;
279: }
280: for (o = -numFaces; o < numFaces; ++o) {
281: Mat symMat;
283: PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, o, &symMat);
284: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_SELF, "Interior node symmetry matrix for orientation %" PetscInt_FMT ":\n", o);
285: PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_SELF);
286: MatView(symMat, PETSC_VIEWER_STDOUT_SELF);
287: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);
288: MatDestroy(&symMat);
289: }
290: }
291: PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_SELF);
292: }
293: PetscDualSpaceDestroy(&sp);
294: return 0;
295: }
297: int main(int argc, char **argv)
298: {
299: PetscInt dim;
300: PetscHashLag lagTable;
301: PetscInt tensorCell;
302: PetscInt order, ordermin, ordermax;
303: PetscBool continuous;
304: PetscBool trimmed;
305: DM dm;
308: PetscInitialize(&argc, &argv, NULL, help);
309: dim = 3;
310: tensorCell = 0;
311: continuous = PETSC_FALSE;
312: trimmed = PETSC_FALSE;
313: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for PETSCDUALSPACELAGRANGE test", "none");
314: PetscOptionsRangeInt("-dim", "The spatial dimension", "ex1.c", dim, &dim, NULL, 0, 3);
315: PetscOptionsRangeInt("-tensor", "(0) simplex (1) hypercube (2) wedge", "ex1.c", tensorCell, &tensorCell, NULL, 0, 2);
316: PetscOptionsBool("-continuous", "Whether the dual space has continuity", "ex1.c", continuous, &continuous, NULL);
317: PetscOptionsBool("-trimmed", "Whether the dual space matches a trimmed polynomial space", "ex1.c", trimmed, &trimmed, NULL);
318: PetscOptionsEnd();
319: PetscHashLagCreate(&lagTable);
321: if (tensorCell < 2) {
322: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, (PetscBool)!tensorCell), &dm);
323: } else {
324: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DM_POLYTOPE_TRI_PRISM, &dm);
325: }
326: ordermin = trimmed ? 1 : 0;
327: ordermax = tensorCell == 2 ? 4 : tensorCell == 1 ? 3 : dim + 2;
328: for (order = ordermin; order <= ordermax; order++) {
329: PetscInt formDegree;
331: for (formDegree = PetscMin(0, -dim + 1); formDegree <= dim; formDegree++) {
332: PetscInt nCopies;
334: for (nCopies = 1; nCopies <= 3; nCopies++) testLagrange(lagTable, dm, dim, order, formDegree, trimmed, (PetscBool)tensorCell, continuous, nCopies);
335: }
336: }
337: DMDestroy(&dm);
338: PetscHashLagDestroy(&lagTable);
339: PetscFinalize();
340: return 0;
341: }
343: /*TEST
345: test:
346: suffix: 0
347: args: -dim 0
349: test:
350: suffix: 1_discontinuous_full
351: args: -dim 1 -continuous 0 -trimmed 0
353: test:
354: suffix: 1_continuous_full
355: args: -dim 1 -continuous 1 -trimmed 0
357: test:
358: suffix: 2_simplex_discontinuous_full
359: args: -dim 2 -tensor 0 -continuous 0 -trimmed 0
361: test:
362: suffix: 2_simplex_continuous_full
363: args: -dim 2 -tensor 0 -continuous 1 -trimmed 0
365: test:
366: suffix: 2_tensor_discontinuous_full
367: args: -dim 2 -tensor 1 -continuous 0 -trimmed 0
369: test:
370: suffix: 2_tensor_continuous_full
371: args: -dim 2 -tensor 1 -continuous 1 -trimmed 0
373: test:
374: suffix: 3_simplex_discontinuous_full
375: args: -dim 3 -tensor 0 -continuous 0 -trimmed 0
377: test:
378: suffix: 3_simplex_continuous_full
379: args: -dim 3 -tensor 0 -continuous 1 -trimmed 0
381: test:
382: suffix: 3_tensor_discontinuous_full
383: args: -dim 3 -tensor 1 -continuous 0 -trimmed 0
385: test:
386: suffix: 3_tensor_continuous_full
387: args: -dim 3 -tensor 1 -continuous 1 -trimmed 0
389: test:
390: suffix: 3_wedge_discontinuous_full
391: args: -dim 3 -tensor 2 -continuous 0 -trimmed 0
393: test:
394: suffix: 3_wedge_continuous_full
395: args: -dim 3 -tensor 2 -continuous 1 -trimmed 0
397: test:
398: suffix: 1_discontinuous_trimmed
399: args: -dim 1 -continuous 0 -trimmed 1
401: test:
402: suffix: 1_continuous_trimmed
403: args: -dim 1 -continuous 1 -trimmed 1
405: test:
406: suffix: 2_simplex_discontinuous_trimmed
407: args: -dim 2 -tensor 0 -continuous 0 -trimmed 1
409: test:
410: suffix: 2_simplex_continuous_trimmed
411: args: -dim 2 -tensor 0 -continuous 1 -trimmed 1
413: test:
414: suffix: 2_tensor_discontinuous_trimmed
415: args: -dim 2 -tensor 1 -continuous 0 -trimmed 1
417: test:
418: suffix: 2_tensor_continuous_trimmed
419: args: -dim 2 -tensor 1 -continuous 1 -trimmed 1
421: test:
422: suffix: 3_simplex_discontinuous_trimmed
423: args: -dim 3 -tensor 0 -continuous 0 -trimmed 1
425: test:
426: suffix: 3_simplex_continuous_trimmed
427: args: -dim 3 -tensor 0 -continuous 1 -trimmed 1
429: test:
430: suffix: 3_tensor_discontinuous_trimmed
431: args: -dim 3 -tensor 1 -continuous 0 -trimmed 1
433: test:
434: suffix: 3_tensor_continuous_trimmed
435: args: -dim 3 -tensor 1 -continuous 1 -trimmed 1
437: test:
438: suffix: 3_wedge_discontinuous_trimmed
439: args: -dim 3 -tensor 2 -continuous 0 -trimmed 1
441: test:
442: suffix: 3_wedge_continuous_trimmed
443: args: -dim 3 -tensor 2 -continuous 1 -trimmed 1
445: TEST*/