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