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