Actual source code: ex4.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      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*/