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 DM dm; 8 PetscDualSpace sp; 9 PetscInt nFunc, *ids, *idsCopy, *idsCopy2, i, closureSize, *closure = NULL, offset, depth; 10 DMLabel depthLabel; 11 PetscBool printed = PETSC_FALSE; 12 PetscScalar *vals, *valsCopy, *valsCopy2; 13 const PetscInt *numDofs; 14 const PetscInt ***perms = NULL; 15 const PetscScalar ***flips = NULL; 16 17 PetscFunctionBegin; 18 PetscCall(PetscDualSpaceCreate(PETSC_COMM_SELF, &sp)); 19 PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, tensor ? PETSC_FALSE : PETSC_TRUE), &dm)); 20 PetscCall(PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE)); 21 PetscCall(PetscDualSpaceSetDM(sp, dm)); 22 PetscCall(PetscDualSpaceSetOrder(sp, order)); 23 PetscCall(PetscDualSpaceLagrangeSetContinuity(sp, PETSC_TRUE)); 24 PetscCall(PetscDualSpaceLagrangeSetTensor(sp, tensor)); 25 PetscCall(PetscDualSpaceSetFromOptions(sp)); 26 PetscCall(PetscDualSpaceSetUp(sp)); 27 PetscCall(PetscDualSpaceGetDimension(sp, &nFunc)); 28 PetscCall(PetscDualSpaceGetSymmetries(sp, &perms, &flips)); 29 if (!perms && !flips) { 30 PetscCall(PetscDualSpaceDestroy(&sp)); 31 PetscCall(DMDestroy(&dm)); 32 PetscFunctionReturn(0); 33 } 34 PetscCall(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; 41 42 PetscCall(PetscDualSpaceGetFunctional(sp, i, &q)); 43 PetscCall(PetscQuadratureGetData(q, NULL, &Nc, &numPoints, &points, &weights)); 44 PetscCheck(Nc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support scalar quadrature, not %" PetscInt_FMT " components", Nc); 45 for (j = 0; j < dim; j++) vals[dim * i + j] = valsCopy2[dim * i + j] = (PetscScalar)points[j]; 46 } 47 PetscCall(PetscDualSpaceGetNumDof(sp, &numDofs)); 48 PetscCall(DMPlexGetDepth(dm, &depth)); 49 PetscCall(DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure)); 50 PetscCall(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; 56 57 if (!pointPerms && !pointFlips) continue; 58 PetscCall(DMLabelGetValue(depthLabel, point, &depth)); 59 { 60 DMPolytopeType ct; 61 /* The number of arrangements is no longer based on the number of faces */ 62 PetscCall(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; 69 70 for (k = 0; k < numDofs[depth]; k++) { 71 PetscInt kLocal = perm ? perm[k] : k; 72 73 idsCopy[kLocal] = ids[offset + k]; 74 for (l = 0; l < dim; l++) { valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.); } 75 } 76 if (!printed && numDofs[depth] > 1) { 77 IS is; 78 Vec vec; 79 char name[256]; 80 81 anyPrinted = PETSC_TRUE; 82 PetscCall(PetscSNPrintf(name, 256, "%" PetscInt_FMT "D, %s, Order %" PetscInt_FMT ", Point %" PetscInt_FMT " Symmetry %" PetscInt_FMT, dim, tensor ? "Tensor" : "Simplex", order, point, j)); 83 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs[depth], idsCopy, PETSC_USE_POINTER, &is)); 84 PetscCall(PetscObjectSetName((PetscObject)is, name)); 85 PetscCall(ISView(is, PETSC_VIEWER_STDOUT_SELF)); 86 PetscCall(ISDestroy(&is)); 87 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dim, numDofs[depth] * dim, valsCopy, &vec)); 88 PetscCall(PetscObjectSetName((PetscObject)vec, name)); 89 PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_SELF)); 90 PetscCall(VecDestroy(&vec)); 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++) { valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.); } 97 } 98 for (k = 0; k < nFunc; k++) { 99 PetscCheck(idsCopy2[k] == ids[k], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Symmetry failure: %" PetscInt_FMT "D, %s, point %" PetscInt_FMT ", symmetry %" PetscInt_FMT ", order %" PetscInt_FMT ", functional %" PetscInt_FMT ": (%" PetscInt_FMT " != %" PetscInt_FMT ")", dim, tensor ? "Tensor" : "Simplex", point, j, order, k, ids[k], k); 100 for (l = 0; l < dim; l++) { 101 PetscCheck(valsCopy2[dim * k + l] == vals[dim * k + l], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Symmetry failure: %" PetscInt_FMT "D, %s, point %" PetscInt_FMT ", symmetry %" PetscInt_FMT ", order %" PetscInt_FMT ", functional %" PetscInt_FMT ", component %" PetscInt_FMT ": (%g != %g)", dim, tensor ? "Tensor" : "Simplex", point, j, order, k, l, (double)PetscAbsScalar(valsCopy2[dim * k + l]), (double)PetscAbsScalar(vals[dim * k + l])); 102 } 103 } 104 } 105 if (anyPrinted && !printed) printed = PETSC_TRUE; 106 } 107 PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure)); 108 PetscCall(PetscFree6(ids, idsCopy, idsCopy2, vals, valsCopy, valsCopy2)); 109 PetscCall(PetscDualSpaceDestroy(&sp)); 110 PetscCall(DMDestroy(&dm)); 111 PetscFunctionReturn(0); 112 } 113 114 int main(int argc, char **argv) { 115 PetscInt dim, order, tensor; 116 117 PetscFunctionBeginUser; 118 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 119 for (tensor = 0; tensor < 2; tensor++) { 120 for (dim = 1; dim <= 3; dim++) { 121 if (dim == 1 && tensor) continue; 122 for (order = 0; order <= (tensor ? 5 : 6); order++) { PetscCall(CheckSymmetry(dim, order, tensor ? PETSC_TRUE : PETSC_FALSE)); } 123 } 124 } 125 PetscCall(PetscFinalize()); 126 return 0; 127 } 128 129 /*TEST 130 test: 131 suffix: 0 132 TEST*/ 133