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: PetscFunctionBegin;
19: PetscCall(PetscDualSpaceCreate(PETSC_COMM_SELF, &sp));
20: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, tensor ? PETSC_FALSE : PETSC_TRUE), &dm));
21: PetscCall(PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE));
22: PetscCall(PetscDualSpaceSetDM(sp, dm));
23: PetscCall(PetscDualSpaceSetOrder(sp, order));
24: PetscCall(PetscDualSpaceLagrangeSetContinuity(sp, PETSC_TRUE));
25: PetscCall(PetscDualSpaceLagrangeSetTensor(sp, tensor));
26: PetscCall(PetscDualSpaceSetFromOptions(sp));
27: PetscCall(PetscDualSpaceSetUp(sp));
28: PetscCall(PetscDualSpaceGetDimension(sp, &nFunc));
29: PetscCall(PetscDualSpaceGetSymmetries(sp, &perms, &flips));
30: if (!perms && !flips) {
31: PetscCall(PetscDualSpaceDestroy(&sp));
32: PetscCall(DMDestroy(&dm));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
35: PetscCall(PetscMalloc6(nFunc, &ids, nFunc, &idsCopy, nFunc, &idsCopy2, nFunc * dim, &vals, nFunc * dim, &valsCopy, nFunc * dim, &valsCopy2));
36: for (i = 0; i < nFunc; i++) ids[i] = idsCopy2[i] = i;
37: for (i = 0; i < nFunc; i++) {
38: PetscQuadrature q;
39: PetscInt numPoints, Nc, j;
40: const PetscReal *points;
41: const PetscReal *weights;
43: PetscCall(PetscDualSpaceGetFunctional(sp, i, &q));
44: PetscCall(PetscQuadratureGetData(q, NULL, &Nc, &numPoints, &points, &weights));
45: PetscCheck(Nc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support scalar quadrature, not %" PetscInt_FMT " components", Nc);
46: for (j = 0; j < dim; j++) vals[dim * i + j] = valsCopy2[dim * i + j] = (PetscScalar)points[j];
47: }
48: PetscCall(PetscDualSpaceGetNumDof(sp, &numDofs));
49: PetscCall(DMPlexGetDepth(dm, &depth));
50: PetscCall(DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure));
51: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
52: for (i = 0, offset = 0; i < closureSize; i++, offset += numDofs[depth]) {
53: PetscInt point = closure[2 * i], numFaces, j;
54: const PetscInt **pointPerms = perms ? perms[i] : NULL;
55: const PetscScalar **pointFlips = flips ? flips[i] : NULL;
56: PetscBool anyPrinted = PETSC_FALSE;
58: if (!pointPerms && !pointFlips) continue;
59: PetscCall(DMLabelGetValue(depthLabel, point, &depth));
60: {
61: DMPolytopeType ct;
62: /* The number of arrangements is no longer based on the number of faces */
63: PetscCall(DMPlexGetCellType(dm, point, &ct));
64: numFaces = DMPolytopeTypeGetNumArrangements(ct) / 2;
65: }
66: for (j = -numFaces; j < numFaces; j++) {
67: PetscInt k, l;
68: const PetscInt *perm = pointPerms ? pointPerms[j] : NULL;
69: const PetscScalar *flip = pointFlips ? pointFlips[j] : NULL;
71: for (k = 0; k < numDofs[depth]; k++) {
72: PetscInt kLocal = perm ? perm[k] : k;
74: idsCopy[kLocal] = ids[offset + k];
75: for (l = 0; l < dim; l++) valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
76: }
77: if (!printed && numDofs[depth] > 1) {
78: IS is;
79: Vec vec;
80: char name[256];
82: anyPrinted = PETSC_TRUE;
83: PetscCall(PetscSNPrintf(name, 256, "%" PetscInt_FMT "D, %s, Order %" PetscInt_FMT ", Point %" PetscInt_FMT " Symmetry %" PetscInt_FMT, dim, tensor ? "Tensor" : "Simplex", order, point, j));
84: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs[depth], idsCopy, PETSC_USE_POINTER, &is));
85: PetscCall(PetscObjectSetName((PetscObject)is, name));
86: PetscCall(ISView(is, PETSC_VIEWER_STDOUT_SELF));
87: PetscCall(ISDestroy(&is));
88: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dim, numDofs[depth] * dim, valsCopy, &vec));
89: PetscCall(PetscObjectSetName((PetscObject)vec, name));
90: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_SELF));
91: PetscCall(VecDestroy(&vec));
92: }
93: for (k = 0; k < numDofs[depth]; k++) {
94: PetscInt kLocal = perm ? perm[k] : k;
96: idsCopy2[offset + k] = idsCopy[kLocal];
97: for (l = 0; l < dim; l++) valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
98: }
99: for (k = 0; k < nFunc; k++) {
100: PetscCheck(idsCopy2[k] == ids[k], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Symmetry failure: %" PetscInt_FMT "D, %s, point %" PetscInt_FMT ", symmetry %" PetscInt_FMT ", order %" PetscInt_FMT ", functional %" PetscInt_FMT ": (%" PetscInt_FMT " != %" PetscInt_FMT ")", dim, tensor ? "Tensor" : "Simplex", point, j, order, k, ids[k], k);
101: for (l = 0; l < dim; l++) {
102: PetscCheck(valsCopy2[dim * k + l] == vals[dim * k + l], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Symmetry failure: %" PetscInt_FMT "D, %s, point %" PetscInt_FMT ", symmetry %" PetscInt_FMT ", order %" PetscInt_FMT ", functional %" PetscInt_FMT ", component %" PetscInt_FMT ": (%g != %g)", dim, tensor ? "Tensor" : "Simplex", point, j, order, k, l, (double)PetscAbsScalar(valsCopy2[dim * k + l]), (double)PetscAbsScalar(vals[dim * k + l]));
103: }
104: }
105: }
106: if (anyPrinted && !printed) printed = PETSC_TRUE;
107: }
108: PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure));
109: PetscCall(PetscFree6(ids, idsCopy, idsCopy2, vals, valsCopy, valsCopy2));
110: PetscCall(PetscDualSpaceDestroy(&sp));
111: PetscCall(DMDestroy(&dm));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: int main(int argc, char **argv)
116: {
117: PetscInt dim, order, tensor;
119: PetscFunctionBeginUser;
120: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
121: for (tensor = 0; tensor < 2; tensor++) {
122: for (dim = 1; dim <= 3; dim++) {
123: if (dim == 1 && tensor) continue;
124: for (order = 0; order <= (tensor ? 5 : 6); order++) PetscCall(CheckSymmetry(dim, order, tensor ? PETSC_TRUE : PETSC_FALSE));
125: }
126: }
127: PetscCall(PetscFinalize());
128: return 0;
129: }
131: /*TEST
132: test:
133: suffix: 0
134: TEST*/