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:       PetscTabulation Tsub, Tfull;
 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:       PetscDualSpaceGetSection(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:       PetscFECreateTabulation(traceFE,1,1,testSub,0,&Tsub);
 91:       PetscFECreateTabulation(fe,1,1,testFull,0,&Tfull);
 92:       for (i = 0; i < Nc; i++) {
 93:         outSub[i] = 0.0;
 94:         for (j = 0; j < nSub; j++) {
 95:           outSub[i] += Tsub->T[0][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] += Tfull->T[0][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:       PetscTabulationDestroy(&Tfull);
116:       PetscTabulationDestroy(&Tsub);
117:       /* clean up */
118:       PetscFree(arraySub);
119:       VecDestroy(&vecFull);
120:       PetscRandomDestroy(&rand);
121:       PetscFree4(testSub,testFull,outSub,outFull);
122:       PetscFEDestroy(&traceFE);
123:     }
124:   }
125:   PetscFEDestroy(&fe);
126:   PetscFinalize();
127:   return ierr;
128: }

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