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 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 ierr = PetscDualSpaceCreate(PETSC_COMM_SELF,&sp);CHKERRQ(ierr); 21 ierr = DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(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], numFaces, j; 55 const PetscInt **pointPerms = perms ? perms[i] : NULL; 56 const PetscScalar **pointFlips = flips ? flips[i] : NULL; 57 PetscBool anyPrinted = PETSC_FALSE; 58 59 if (!pointPerms && !pointFlips) continue; 60 ierr = DMLabelGetValue(depthLabel,point,&depth);CHKERRQ(ierr); 61 { 62 DMPolytopeType ct; 63 /* The number of arrangements is no longer based on the number of faces */ 64 ierr = DMPlexGetCellType(dm, point, &ct);CHKERRQ(ierr); 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; 71 72 for (k = 0; k < numDofs[depth]; k++) { 73 PetscInt kLocal = perm ? perm[k] : k; 74 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]; 84 85 anyPrinted = PETSC_TRUE; 86 ierr = PetscSNPrintf(name,256,"%DD, %s, Order %D, Point %D Symmetry %D",dim,tensor ? "Tensor" : "Simplex", order, point,j);CHKERRQ(ierr); 87 ierr = ISCreateGeneral(PETSC_COMM_SELF,numDofs[depth],idsCopy,PETSC_USE_POINTER,&is);CHKERRQ(ierr); 88 ierr = PetscObjectSetName((PetscObject)is,name);CHKERRQ(ierr); 89 ierr = ISView(is,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 90 ierr = ISDestroy(&is);CHKERRQ(ierr); 91 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dim,numDofs[depth]*dim,valsCopy,&vec);CHKERRQ(ierr); 92 ierr = PetscObjectSetName((PetscObject)vec,name);CHKERRQ(ierr); 93 ierr = VecView(vec,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 94 ierr = VecDestroy(&vec);CHKERRQ(ierr); 95 } 96 for (k = 0; k < numDofs[depth]; k++) { 97 PetscInt kLocal = perm ? perm[k] : k; 98 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 ierr = DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 114 ierr = PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);CHKERRQ(ierr); 115 ierr = PetscDualSpaceDestroy(&sp);CHKERRQ(ierr); 116 ierr = DMDestroy(&dm);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 int main(int argc, char **argv) 121 { 122 PetscInt dim, order, tensor; 123 PetscErrorCode ierr; 124 125 ierr = 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 ierr = CheckSymmetry(dim,order,tensor ? PETSC_TRUE : PETSC_FALSE);CHKERRQ(ierr); 131 } 132 } 133 } 134 ierr = PetscFinalize(); 135 return ierr; 136 } 137 138 /*TEST 139 test: 140 suffix: 0 141 TEST*/ 142