1 static char help[] = "Tests dual space symmetry.\n\n"; 2 3 #include <petscfe.h> 4 #include <petscdmplex.h> 5 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 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(0); 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; 42 43 PetscCall(PetscDualSpaceGetFunctional(sp,i,&q)); 44 PetscCall(PetscQuadratureGetData(q,NULL,&Nc,&numPoints,&points,&weights)); 45 PetscCheckFalse(Nc != 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only support scalar quadrature, not %D 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; 57 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 = DMPolytopeTypeGetNumArrangments(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; 70 71 for (k = 0; k < numDofs[depth]; k++) { 72 PetscInt kLocal = perm ? perm[k] : k; 73 74 idsCopy[kLocal] = ids[offset + k]; 75 for (l = 0; l < dim; l++) { 76 valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.); 77 } 78 } 79 if (!printed && numDofs[depth] > 1) { 80 IS is; 81 Vec vec; 82 char name[256]; 83 84 anyPrinted = PETSC_TRUE; 85 PetscCall(PetscSNPrintf(name,256,"%DD, %s, Order %D, Point %D Symmetry %D",dim,tensor ? "Tensor" : "Simplex", order, point,j)); 86 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,numDofs[depth],idsCopy,PETSC_USE_POINTER,&is)); 87 PetscCall(PetscObjectSetName((PetscObject)is,name)); 88 PetscCall(ISView(is,PETSC_VIEWER_STDOUT_SELF)); 89 PetscCall(ISDestroy(&is)); 90 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,dim,numDofs[depth]*dim,valsCopy,&vec)); 91 PetscCall(PetscObjectSetName((PetscObject)vec,name)); 92 PetscCall(VecView(vec,PETSC_VIEWER_STDOUT_SELF)); 93 PetscCall(VecDestroy(&vec)); 94 } 95 for (k = 0; k < numDofs[depth]; k++) { 96 PetscInt kLocal = perm ? perm[k] : k; 97 98 idsCopy2[offset + k] = idsCopy[kLocal]; 99 for (l = 0; l < dim; l++) { 100 valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.); 101 } 102 } 103 for (k = 0; k < nFunc; k++) { 104 PetscCheckFalse(idsCopy2[k] != ids[k],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); 105 for (l = 0; l < dim; l++) { 106 PetscCheckFalse(valsCopy2[dim * k + l] != vals[dim * k + l],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); 107 } 108 } 109 } 110 if (anyPrinted && !printed) printed = PETSC_TRUE; 111 } 112 PetscCall(DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure)); 113 PetscCall(PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2)); 114 PetscCall(PetscDualSpaceDestroy(&sp)); 115 PetscCall(DMDestroy(&dm)); 116 PetscFunctionReturn(0); 117 } 118 119 int main(int argc, char **argv) 120 { 121 PetscInt dim, order, tensor; 122 123 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 124 for (tensor = 0; tensor < 2; tensor++) { 125 for (dim = 1; dim <= 3; dim++) { 126 if (dim == 1 && tensor) continue; 127 for (order = 0; order <= (tensor ? 5 : 6); order++) { 128 PetscCall(CheckSymmetry(dim,order,tensor ? PETSC_TRUE : PETSC_FALSE)); 129 } 130 } 131 } 132 PetscCall(PetscFinalize()); 133 return 0; 134 } 135 136 /*TEST 137 test: 138 suffix: 0 139 TEST*/ 140