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++) {
75: valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
76: }
77: }
78: if (!printed && numDofs[depth] > 1) {
79: IS is;
80: Vec vec;
81: char name[256];
83: anyPrinted = PETSC_TRUE;
84: PetscSNPrintf(name,256,"%DD, %s, Order %D, Point %D Symmetry %D",dim,tensor ? "Tensor" : "Simplex", order, point,j);
85: ISCreateGeneral(PETSC_COMM_SELF,numDofs[depth],idsCopy,PETSC_USE_POINTER,&is);
86: PetscObjectSetName((PetscObject)is,name);
87: ISView(is,PETSC_VIEWER_STDOUT_SELF);
88: ISDestroy(&is);
89: VecCreateSeqWithArray(PETSC_COMM_SELF,dim,numDofs[depth]*dim,valsCopy,&vec);
90: PetscObjectSetName((PetscObject)vec,name);
91: VecView(vec,PETSC_VIEWER_STDOUT_SELF);
92: VecDestroy(&vec);
93: }
94: for (k = 0; k < numDofs[depth]; k++) {
95: PetscInt kLocal = perm ? perm[k] : k;
97: idsCopy2[offset + k] = idsCopy[kLocal];
98: for (l = 0; l < dim; l++) {
99: valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
100: }
101: }
102: for (k = 0; k < nFunc; k++) {
104: for (l = 0; l < dim; l++) {
106: }
107: }
108: }
109: if (anyPrinted && !printed) printed = PETSC_TRUE;
110: }
111: DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);
112: PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);
113: PetscDualSpaceDestroy(&sp);
114: DMDestroy(&dm);
115: return 0;
116: }
118: int main(int argc, char **argv)
119: {
120: PetscInt dim, order, tensor;
122: PetscInitialize(&argc,&argv,NULL,help);
123: for (tensor = 0; tensor < 2; tensor++) {
124: for (dim = 1; dim <= 3; dim++) {
125: if (dim == 1 && tensor) continue;
126: for (order = 0; order <= (tensor ? 5 : 6); order++) {
127: CheckSymmetry(dim,order,tensor ? PETSC_TRUE : PETSC_FALSE);
128: }
129: }
130: }
131: PetscFinalize();
132: return 0;
133: }
135: /*TEST
136: test:
137: suffix: 0
138: TEST*/