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++) {
 75:           valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
 76:         }
 77:       }
 78:       if (!printed && numDofs[depth] > 1) {
 79:         IS   is;
 80:         Vec  vec;
 81:         char name[256];

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

 97:         idsCopy2[offset + k] = idsCopy[kLocal];
 98:         for (l = 0; l < dim; l++) {
 99:           valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
100:         }
101:       }
102:       for (k = 0; k < nFunc; k++) {
104:         for (l = 0; l < dim; l++) {
106:         }
107:       }
108:     }
109:     if (anyPrinted && !printed) printed = PETSC_TRUE;
110:   }
111:   DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);
112:   PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);
113:   PetscDualSpaceDestroy(&sp);
114:   DMDestroy(&dm);
115:   return 0;
116: }

118: int main(int argc, char **argv)
119: {
120:   PetscInt       dim, order, tensor;

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

135: /*TEST
136:   test:
137:     suffix: 0
138: TEST*/