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

      1: static char help[] = "Tests affine subspaces.\n\n";

      3:  #include <petscfe.h>
      4:  #include <petscdmplex.h>
      5:  #include <petscdmshell.h>

      7: int main(int argc, char **argv)
      8: {
      9:   DM             dm;
     10:   PetscFE        fe;
     11:   PetscSpace     space;
     12:   PetscDualSpace dualspace, dualsubspace;
     13:   PetscInt       dim = 2, Nc = 3, cStart, cEnd;
     14:   PetscBool      simplex = PETSC_TRUE;
     15:   MPI_Comm       comm;

     18:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
     19:   comm = PETSC_COMM_WORLD;
     20:   PetscOptionsBegin(comm,"","Options for subspace test","none");
     21:   PetscOptionsRangeInt("-dim", "The spatial dimension","ex5.c",dim,&dim,NULL,1,3);
     22:   PetscOptionsBool("-simplex", "Test simplex element","ex5.c",simplex,&simplex,NULL);
     23:   PetscOptionsBoundedInt("-num_comp", "Number of components in space","ex5.c",Nc,&Nc,NULL,1);
     24:   PetscOptionsEnd();
     25:   DMShellCreate(comm,&dm);
     26:   PetscFECreateDefault(comm,dim,Nc,simplex,NULL,PETSC_DEFAULT,&fe);
     27:   DMDestroy(&dm);
     28:   PetscFESetName(fe, "solution");
     29:   PetscFEGetBasisSpace(fe,&space);
     30:   PetscSpaceGetNumComponents(space,&Nc);
     31:   PetscFEGetDualSpace(fe,&dualspace);
     32:   PetscDualSpaceGetHeightSubspace(dualspace,1,&dualsubspace);
     33:   PetscDualSpaceGetDM(dualspace,&dm);
     34:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
     35:   if (cEnd > cStart) {
     36:     PetscInt coneSize;

     38:     DMPlexGetConeSize(dm,cStart,&coneSize);
     39:     if (coneSize) {
     40:       PetscFE traceFE;
     41:       const PetscInt *cone;
     42:       PetscInt        point, nSub, nFull;
     43:       PetscReal       xi0[3] = {-1., -1., -1.};
     44:       PetscScalar     *outSub, *outFull;
     45:       PetscReal       *testSub, *testFull;
     46:       PetscReal       *Bsub, *Bfull;
     47:       PetscReal       J[9], detJ;
     48:       PetscInt        i, j;
     49:       PetscSection    sectionFull;
     50:       Vec             vecFull;
     51:       PetscScalar     *arrayFull, *arraySub;
     52:       PetscReal       err;
     53:       PetscRandom     rand;

     55:       DMPlexGetCone(dm,cStart,&cone);
     56:       point = cone[0];
     57:       PetscFECreatePointTrace(fe,point,&traceFE);
     58:       PetscFESetUp(traceFE);
     59:       PetscFEViewFromOptions(traceFE,NULL,"-trace_fe_view");
     60:       PetscMalloc4(dim - 1,&testSub,dim,&testFull,Nc,&outSub,Nc,&outFull);
     61:       PetscRandomCreate(PETSC_COMM_SELF,&rand);
     62:       PetscRandomSetFromOptions(rand);
     63:       PetscRandomSetInterval(rand,-1.,1.);
     64:       /* create a random point in the trace domain */
     65:       for (i = 0; i < dim - 1; i++) {
     66:         PetscRandomGetValueReal(rand,&testSub[i]);
     67:       }
     68:       DMPlexComputeCellGeometryFEM(dm,point,NULL,testFull,J,NULL,&detJ);
     69:       /* project it into the full domain */
     70:       for (i = 0; i < dim; i++) {
     71:         for (j = 0; j < dim - 1; j++) testFull[i] += J[i * dim + j] * (testSub[j] - xi0[j]);
     72:       }
     73:       /* create a random vector in the full domain */
     74:       PetscFEGetDimension(fe,&nFull);
     75:       VecCreateSeq(PETSC_COMM_SELF,nFull,&vecFull);
     76:       VecGetArray(vecFull,&arrayFull);
     77:       for (i = 0; i < nFull; i++) {
     78:         PetscRandomGetValue(rand,&arrayFull[i]);
     79:       }
     80:       VecRestoreArray(vecFull,&arrayFull);
     81:       /* create a vector on the trace domain */
     82:       PetscFEGetDimension(traceFE,&nSub);
     83:       /* get the subset of the original finite element space that is supported on the trace space */
     84:       PetscDualSpaceCreateSection(dualspace,&sectionFull);
     85:       PetscSectionSetUp(sectionFull);
     86:       /* get the trace degrees of freedom */
     87:       PetscMalloc1(nSub,&arraySub);
     88:       DMPlexVecGetClosure(dm,sectionFull,vecFull,point,&nSub,&arraySub);
     89:       /* get the tabulations */
     90:       PetscFEGetTabulation(traceFE,1,testSub,&Bsub,NULL,NULL);
     91:       PetscFEGetTabulation(fe,1,testFull,&Bfull,NULL,NULL);
     92:       for (i = 0; i < Nc; i++) {
     93:         outSub[i] = 0.0;
     94:         for (j = 0; j < nSub; j++) {
     95:           outSub[i] += Bsub[j * Nc + i] * arraySub[j];
     96:         }
     97:       }
     98:       VecGetArray(vecFull,&arrayFull);
     99:       err = 0.0;
    100:       for (i = 0; i < Nc; i++) {
    101:         PetscScalar diff;

    103:         outFull[i] = 0.0;
    104:         for (j = 0; j < nFull; j++) {
    105:           outFull[i] += Bfull[j * Nc + i] * arrayFull[j];
    106:         }
    107:         diff = outFull[i] - outSub[i];
    108:         err += PetscRealPart(PetscConj(diff) * diff);
    109:       }
    110:       err = PetscSqrtReal(err);
    111:       if (err > PETSC_SMALL) {
    112:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Trace FE error %g\n",err);
    113:       }
    114:       VecRestoreArray(vecFull,&arrayFull);
    115:       PetscFERestoreTabulation(fe,1,testFull,&Bfull,NULL,NULL);
    116:       PetscFERestoreTabulation(traceFE,1,testSub,&Bsub,NULL,NULL);
    117:       /* clean up */
    118:       PetscFree(arraySub);
    119:       PetscSectionDestroy(&sectionFull);
    120:       VecDestroy(&vecFull);
    121:       PetscRandomDestroy(&rand);
    122:       PetscFree4(testSub,testFull,outSub,outFull);
    123:       PetscFEDestroy(&traceFE);
    124:     }
    125:   }
    126:   PetscFEDestroy(&fe);
    127:   PetscFinalize();
    128:   return ierr;
    129: }

    131: /*TEST
    132:   test:
    133:     suffix: 0
    134:     args: -petscspace_degree 1 -trace_fe_view
    135: TEST*/