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:   PetscDualSpaceCreate(PETSC_COMM_SELF, &sp);
 19:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, tensor ? PETSC_FALSE : PETSC_TRUE), &dm);
 20:   PetscDualSpaceSetType(sp, PETSCDUALSPACELAGRANGE);
 21:   PetscDualSpaceSetDM(sp, dm);
 22:   PetscDualSpaceSetOrder(sp, order);
 23:   PetscDualSpaceLagrangeSetContinuity(sp, PETSC_TRUE);
 24:   PetscDualSpaceLagrangeSetTensor(sp, tensor);
 25:   PetscDualSpaceSetFromOptions(sp);
 26:   PetscDualSpaceSetUp(sp);
 27:   PetscDualSpaceGetDimension(sp, &nFunc);
 28:   PetscDualSpaceGetSymmetries(sp, &perms, &flips);
 29:   if (!perms && !flips) {
 30:     PetscDualSpaceDestroy(&sp);
 31:     DMDestroy(&dm);
 32:     return 0;
 33:   }
 34:   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;

 42:     PetscDualSpaceGetFunctional(sp, i, &q);
 43:     PetscQuadratureGetData(q, NULL, &Nc, &numPoints, &points, &weights);
 45:     for (j = 0; j < dim; j++) vals[dim * i + j] = valsCopy2[dim * i + j] = (PetscScalar)points[j];
 46:   }
 47:   PetscDualSpaceGetNumDof(sp, &numDofs);
 48:   DMPlexGetDepth(dm, &depth);
 49:   DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
 50:   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;

 57:     if (!pointPerms && !pointFlips) continue;
 58:     DMLabelGetValue(depthLabel, point, &depth);
 59:     {
 60:       DMPolytopeType ct;
 61:       /* The number of arrangements is no longer based on the number of faces */
 62:       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;

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

 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];

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

 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++) {
100:         for (l = 0; l < dim; l++) {
102:         }
103:       }
104:     }
105:     if (anyPrinted && !printed) printed = PETSC_TRUE;
106:   }
107:   DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);
108:   PetscFree6(ids, idsCopy, idsCopy2, vals, valsCopy, valsCopy2);
109:   PetscDualSpaceDestroy(&sp);
110:   DMDestroy(&dm);
111:   return 0;
112: }

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

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

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