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,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], coneSize, j;
 55:     const PetscInt    **pointPerms = perms ? perms[i] : NULL;
 56:     const PetscScalar **pointFlips = flips ? flips[i] : NULL;
 57:     PetscBool         anyPrinted = PETSC_FALSE;

 59:     DMLabelGetValue(depthLabel,point,&depth);
 60:     DMPlexGetConeSize(dm,point,&coneSize);

 62:     if (!pointPerms && !pointFlips) continue;
 63:     for (j = -coneSize; j < coneSize; j++) {
 64:       PetscInt          k, l;
 65:       const PetscInt    *perm = pointPerms ? pointPerms[j] : NULL;
 66:       const PetscScalar *flip = pointFlips ? pointFlips[j] : NULL;

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

 71:         idsCopy[kLocal] = ids[offset + k];
 72:         for (l = 0; l < dim; l++) {
 73:           valsCopy[kLocal * dim + l] = vals[(offset + k) * dim + l] * (flip ? flip[kLocal] : 1.);
 74:         }
 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,"%DD, %s, Order %D, Point %D Symmetry %D",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++) {
 97:           valsCopy2[(offset + k) * dim + l] = valsCopy[kLocal * dim + l] * (flip ? PetscConj(flip[kLocal]) : 1.);
 98:         }
 99:       }
100:       for (k = 0; k < nFunc; k++) {
101:         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);
102:         for (l = 0; l < dim; l++) {
103:           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);
104:         }
105:       }
106:     }
107:     if (anyPrinted && !printed) printed = PETSC_TRUE;
108:   }
109:   DMPlexRestoreTransitiveClosure(dm,0,PETSC_TRUE,&closureSize,&closure);
110:   PetscFree6(ids,idsCopy,idsCopy2,vals,valsCopy,valsCopy2);
111:   PetscDualSpaceDestroy(&sp);
112:   DMDestroy(&dm);
113:   return(0);
114: }

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

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

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