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;

 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(PETSC_SUCCESS);
 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;

 43:     PetscCall(PetscDualSpaceGetFunctional(sp, i, &q));
 44:     PetscCall(PetscQuadratureGetData(q, NULL, &Nc, &numPoints, &points, &weights));
 45:     PetscCheck(Nc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support scalar quadrature, not %" PetscInt_FMT " 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;

 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 = DMPolytopeTypeGetNumArrangements(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;

 71:       for (k = 0; k < numDofs[depth]; k++) {
 72:         PetscInt kLocal = perm ? perm[k] : k;

 74:         idsCopy[kLocal] = ids[offset + k];
 75:         for (l = 0; l < dim; l++) valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
 76:       }
 77:       if (!printed && numDofs[depth] > 1) {
 78:         IS   is;
 79:         Vec  vec;
 80:         char name[256];

 82:         anyPrinted = PETSC_TRUE;
 83:         PetscCall(PetscSNPrintf(name, 256, "%" PetscInt_FMT "D, %s, Order %" PetscInt_FMT ", Point %" PetscInt_FMT " Symmetry %" PetscInt_FMT, dim, tensor ? "Tensor" : "Simplex", order, point, j));
 84:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numDofs[depth], idsCopy, PETSC_USE_POINTER, &is));
 85:         PetscCall(PetscObjectSetName((PetscObject)is, name));
 86:         PetscCall(ISView(is, PETSC_VIEWER_STDOUT_SELF));
 87:         PetscCall(ISDestroy(&is));
 88:         PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dim, numDofs[depth] * dim, valsCopy, &vec));
 89:         PetscCall(PetscObjectSetName((PetscObject)vec, name));
 90:         PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_SELF));
 91:         PetscCall(VecDestroy(&vec));
 92:       }
 93:       for (k = 0; k < numDofs[depth]; k++) {
 94:         PetscInt kLocal = perm ? perm[k] : k;

 96:         idsCopy2[offset + k] = idsCopy[kLocal];
 97:         for (l = 0; l < dim; l++) valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
 98:       }
 99:       for (k = 0; k < nFunc; k++) {
100:         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);
101:         for (l = 0; l < dim; l++) {
102:           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]));
103:         }
104:       }
105:     }
106:     if (anyPrinted && !printed) printed = PETSC_TRUE;
107:   }
108:   PetscCall(DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure));
109:   PetscCall(PetscFree6(ids, idsCopy, idsCopy2, vals, valsCopy, valsCopy2));
110:   PetscCall(PetscDualSpaceDestroy(&sp));
111:   PetscCall(DMDestroy(&dm));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: int main(int argc, char **argv)
116: {
117:   PetscInt dim, order, tensor;

119:   PetscFunctionBeginUser;
120:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
121:   for (tensor = 0; tensor < 2; tensor++) {
122:     for (dim = 1; dim <= 3; dim++) {
123:       if (dim == 1 && tensor) continue;
124:       for (order = 0; order <= (tensor ? 5 : 6); order++) PetscCall(CheckSymmetry(dim, order, tensor ? PETSC_TRUE : PETSC_FALSE));
125:     }
126:   }
127:   PetscCall(PetscFinalize());
128:   return 0;
129: }

131: /*TEST
132:   test:
133:     suffix: 0
134: TEST*/