1 static char help[] = "Tests dual space symmetry.\n\n";
2
3 #include <petscfe.h>
4 #include <petscdmplex.h>
5
CheckSymmetry(PetscInt dim,PetscInt order,PetscBool tensor)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;
18
19 PetscFunctionBegin;
20 ierr = PetscDualSpaceCreate(PETSC_COMM_SELF,&sp);CHKERRQ(ierr);
21 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF,dim,tensor ? PETSC_FALSE : PETSC_TRUE,&dm);CHKERRQ(ierr);
22 ierr = PetscDualSpaceSetType(sp,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr);
23 ierr = PetscDualSpaceSetDM(sp,dm);CHKERRQ(ierr);
24 ierr = PetscDualSpaceSetOrder(sp,order);CHKERRQ(ierr);
25 ierr = PetscDualSpaceLagrangeSetContinuity(sp,PETSC_TRUE);CHKERRQ(ierr);
26 ierr = PetscDualSpaceLagrangeSetTensor(sp,tensor);CHKERRQ(ierr);
27 ierr = PetscDualSpaceSetFromOptions(sp);CHKERRQ(ierr);
28 ierr = PetscDualSpaceSetUp(sp);CHKERRQ(ierr);
29 ierr = PetscDualSpaceGetDimension(sp,&nFunc);CHKERRQ(ierr);
30 ierr = PetscDualSpaceGetSymmetries(sp,&perms,&flips);CHKERRQ(ierr);
31 if (!perms && !flips) {
32 ierr = PetscDualSpaceDestroy(&sp);CHKERRQ(ierr);
33 ierr = DMDestroy(&dm);CHKERRQ(ierr);
34 PetscFunctionReturn(0);
35 }
36 ierr = PetscMalloc6(nFunc,&ids,nFunc,&idsCopy,nFunc,&idsCopy2,nFunc*dim,&vals,nFunc*dim,&valsCopy,nFunc*dim,&valsCopy2);CHKERRQ(ierr);
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;
43
44 ierr = PetscDualSpaceGetFunctional(sp,i,&q);CHKERRQ(ierr);
45 ierr = PetscQuadratureGetData(q,NULL,&Nc,&numPoints,&points,&weights);CHKERRQ(ierr);
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 ierr = PetscDualSpaceGetNumDof(sp,&numDofs);CHKERRQ(ierr);
50 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr);
51 ierr = DMPlexGetTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
52 ierr = DMPlexGetDepthLabel(dm,&depthLabel);CHKERRQ(ierr);
53 for (i = 0, offset = 0; i < closureSize; i++, offset += numDofs[depth]) {
54 PetscInt point = closure[2 * i], coneSize, j;
55 const PetscInt **pointPerms = perms ? perms[i] : NULL;
56 const PetscScalar **pointFlips = flips ? flips[i] : NULL;
57 PetscBool anyPrinted = PETSC_FALSE;
58
59 ierr = DMLabelGetValue(depthLabel,point,&depth);CHKERRQ(ierr);
60 ierr = DMPlexGetConeSize(dm,point,&coneSize);CHKERRQ(ierr);
61
62 if (!pointPerms && !pointFlips) continue;
63 for (j = -coneSize; j < coneSize; j++) {
64 PetscInt k, l;
65 const PetscInt *perm = pointPerms ? pointPerms[j] : NULL;
66 const PetscScalar *flip = pointFlips ? pointFlips[j] : NULL;
67
68 for (k = 0; k < numDofs[depth]; k++) {
69 PetscInt kLocal = perm ? perm[k] : k;
70
71 idsCopy[kLocal] = ids[offset + k];
72 for (l = 0; l < dim; l++) {
73 valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
74 }
75 }
76 if (!printed && numDofs[depth] > 1) {
77 IS is;
78 Vec vec;
79 char name[256];
80
81 anyPrinted = PETSC_TRUE;
82 ierr = PetscSNPrintf(name,256,"%DD, %s, Order %D, Point %D Symmetry %D",dim,tensor ? "Tensor" : "Simplex", order, point,j);CHKERRQ(ierr);
83 ierr = ISCreateGeneral(PETSC_COMM_SELF,numDofs[depth],idsCopy,PETSC_USE_POINTER,&is);CHKERRQ(ierr);
84 ierr = PetscObjectSetName((PetscObject)is,name);CHKERRQ(ierr);
85 ierr = ISView(is,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
86 ierr = ISDestroy(&is);CHKERRQ(ierr);
87 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dim,numDofs[depth]*dim,valsCopy,&vec);CHKERRQ(ierr);
88 ierr = PetscObjectSetName((PetscObject)vec,name);CHKERRQ(ierr);
89 ierr = VecView(vec,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
90 ierr = VecDestroy(&vec);CHKERRQ(ierr);
91 }
92 for (k = 0; k < numDofs[depth]; k++) {
93 PetscInt kLocal = perm ? perm[k] : k;
94
95 idsCopy2[offset + k] = idsCopy[kLocal];
96 for (l = 0; l < dim; l++) {
97 valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
98 }
99 }
100 for (k = 0; k < nFunc; k++) {
101 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);
102 for (l = 0; l < dim; l++) {
103 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);
104 }
105 }
106 }
107 if (anyPrinted && !printed) printed = PETSC_TRUE;
108 }
109 ierr = DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr);
110 ierr = PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);CHKERRQ(ierr);
111 ierr = PetscDualSpaceDestroy(&sp);CHKERRQ(ierr);
112 ierr = DMDestroy(&dm);CHKERRQ(ierr);
113 PetscFunctionReturn(0);
114 }
115
main(int argc,char ** argv)116 int main(int argc, char **argv)
117 {
118 PetscInt dim, order, tensor;
119 PetscErrorCode ierr;
120
121 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
122 for (tensor = 0; tensor < 2; tensor++) {
123 for (dim = 1; dim <= 3; dim++) {
124 if (dim == 1 && tensor) continue;
125 for (order = 0; order <= (tensor ? 5 : 6); order++) {
126 ierr = CheckSymmetry(dim,order,tensor ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr);
127 }
128 }
129 }
130 ierr = PetscFinalize();
131 return ierr;
132 }
133
134 /*TEST
135 test:
136 suffix: 0
137 TEST*/
138