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*/