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;
 17:   PetscErrorCode    ierr;

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

 44:     PetscDualSpaceGetFunctional(sp,i,&q);
 45:     PetscQuadratureGetData(q,NULL,&Nc,&numPoints,&points,&weights);
 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:   PetscDualSpaceGetNumDof(sp,&numDofs);
 50:   DMPlexGetDepth(dm,&depth);
 51:   DMPlexGetTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);
 52:   DMPlexGetDepthLabel(dm,&depthLabel);
 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;

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

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

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

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

 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:   DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);
114:   PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);
115:   PetscDualSpaceDestroy(&sp);
116:   DMDestroy(&dm);
117:   return(0);
118: }

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

125:   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:         CheckSymmetry(dim,order,tensor ? PETSC_TRUE : PETSC_FALSE);
131:       }
132:     }
133:   }
134:   PetscFinalize();
135:   return ierr;
136: }

138: /*TEST
139:   test:
140:     suffix: 0
141: TEST*/