Actual source code: ex4.c
1: static char help[] = "Tests dual space symmetry.\n\n";
3: #include <petscfe.h>
4: #include <petscdmplex.h>
6: static PetscErrorCode CheckSymmetry(PetscInt dim, PetscInt order, PetscBool tensor)
7: {
8: DM dm;
9: PetscDualSpace sp;
10: PetscInt nFunc, *ids, *idsCopy, *idsCopy2, i, closureSize, *closure = NULL, offset, depth;
11: DMLabel depthLabel;
12: PetscBool printed = PETSC_FALSE;
13: PetscScalar *vals, *valsCopy, *valsCopy2;
14: const PetscInt *numDofs;
15: const PetscInt ***perms = NULL;
16: const PetscScalar ***flips = NULL;
18: PetscDualSpaceCreate(PETSC_COMM_SELF, &sp);
19: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, tensor ? PETSC_FALSE : PETSC_TRUE), &dm);
20: PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE);
21: PetscDualSpaceSetDM(sp, dm);
22: PetscDualSpaceSetOrder(sp, order);
23: PetscDualSpaceLagrangeSetContinuity(sp, PETSC_TRUE);
24: PetscDualSpaceLagrangeSetTensor(sp, tensor);
25: PetscDualSpaceSetFromOptions(sp);
26: PetscDualSpaceSetUp(sp);
27: PetscDualSpaceGetDimension(sp, &nFunc);
28: PetscDualSpaceGetSymmetries(sp, &perms, &flips);
29: if (!perms && !flips) {
30: PetscDualSpaceDestroy(&sp);
31: DMDestroy(&dm);
32: return 0;
33: }
34: PetscMalloc6(nFunc, &ids, nFunc, &idsCopy, nFunc, &idsCopy2, nFunc * dim, &vals, nFunc * dim, &valsCopy, nFunc * dim, &valsCopy2);
35: for (i = 0; i < nFunc; i++) ids[i] = idsCopy2[i] = i;
36: for (i = 0; i < nFunc; i++) {
37: PetscQuadrature q;
38: PetscInt numPoints, Nc, j;
39: const PetscReal *points;
40: const PetscReal *weights;
42: PetscDualSpaceGetFunctional(sp, i, &q);
43: PetscQuadratureGetData(q, NULL, &Nc, &numPoints, &points, &weights);
45: for (j = 0; j < dim; j++) vals[dim * i + j] = valsCopy2[dim * i + j] = (PetscScalar)points[j];
46: }
47: PetscDualSpaceGetNumDof(sp, &numDofs);
48: DMPlexGetDepth(dm, &depth);
49: DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
50: DMPlexGetDepthLabel(dm, &depthLabel);
51: for (i = 0, offset = 0; i < closureSize; i++, offset += numDofs[depth]) {
52: PetscInt point = closure[2 * i], numFaces, j;
53: const PetscInt **pointPerms = perms ? perms[i] : NULL;
54: const PetscScalar **pointFlips = flips ? flips[i] : NULL;
55: PetscBool anyPrinted = PETSC_FALSE;
57: if (!pointPerms && !pointFlips) continue;
58: DMLabelGetValue(depthLabel, point, &depth);
59: {
60: DMPolytopeType ct;
61: /* The number of arrangements is no longer based on the number of faces */
62: DMPlexGetCellType(dm, point, &ct);
63: numFaces = DMPolytopeTypeGetNumArrangments(ct) / 2;
64: }
65: for (j = -numFaces; j < numFaces; j++) {
66: PetscInt k, l;
67: const PetscInt *perm = pointPerms ? pointPerms[j] : NULL;
68: const PetscScalar *flip = pointFlips ? pointFlips[j] : NULL;
70: for (k = 0; k < numDofs[depth]; k++) {
71: PetscInt kLocal = perm ? perm[k] : k;
73: idsCopy[kLocal] = ids[offset + k];
74: for (l = 0; l < dim; l++) valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
75: }
76: if (!printed && numDofs[depth] > 1) {
77: IS is;
78: Vec vec;
79: char name[256];
81: anyPrinted = PETSC_TRUE;
82: PetscSNPrintf(name, 256, "%" PetscInt_FMT "D, %s, Order %" PetscInt_FMT ", Point %" PetscInt_FMT " Symmetry %" PetscInt_FMT, dim, tensor ? "Tensor" : "Simplex", order, point, j);
83: ISCreateGeneral(PETSC_COMM_SELF, numDofs[depth], idsCopy, PETSC_USE_POINTER, &is);
84: PetscObjectSetName((PetscObject)is, name);
85: ISView(is, PETSC_VIEWER_STDOUT_SELF);
86: ISDestroy(&is);
87: VecCreateSeqWithArray(PETSC_COMM_SELF, dim, numDofs[depth] * dim, valsCopy, &vec);
88: PetscObjectSetName((PetscObject)vec, name);
89: VecView(vec, PETSC_VIEWER_STDOUT_SELF);
90: VecDestroy(&vec);
91: }
92: for (k = 0; k < numDofs[depth]; k++) {
93: PetscInt kLocal = perm ? perm[k] : k;
95: idsCopy2[offset + k] = idsCopy[kLocal];
96: for (l = 0; l < dim; l++) valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
97: }
98: for (k = 0; k < nFunc; k++) {
100: for (l = 0; l < dim; l++) {
102: }
103: }
104: }
105: if (anyPrinted && !printed) printed = PETSC_TRUE;
106: }
107: DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
108: PetscFree6(ids, idsCopy, idsCopy2, vals, valsCopy, valsCopy2);
109: PetscDualSpaceDestroy(&sp);
110: DMDestroy(&dm);
111: return 0;
112: }
114: int main(int argc, char **argv)
115: {
116: PetscInt dim, order, tensor;
119: PetscInitialize(&argc, &argv, NULL, help);
120: for (tensor = 0; tensor < 2; tensor++) {
121: for (dim = 1; dim <= 3; dim++) {
122: if (dim == 1 && tensor) continue;
123: for (order = 0; order <= (tensor ? 5 : 6); order++) CheckSymmetry(dim, order, tensor ? PETSC_TRUE : PETSC_FALSE);
124: }
125: }
126: PetscFinalize();
127: return 0;
128: }
130: /*TEST
131: test:
132: suffix: 0
133: TEST*/