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;

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

 38:     PetscCall(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:       PetscSection    sectionFull;
 49:       Vec             vecFull;
 50:       PetscScalar    *arrayFull, *arraySub;
 51:       PetscReal       err;
 52:       PetscRandom     rand;

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

 96:         outFull[i] = 0.0;
 97:         for (PetscInt j = 0; j < nFull; j++) outFull[i] += Tfull->T[0][j * Nc + i] * arrayFull[j];
 98:         diff = outFull[i] - outSub[i];
 99:         err += PetscRealPart(PetscConj(diff) * diff);
100:       }
101:       err = PetscSqrtReal(err);
102:       PetscCheck(err <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Trace FE error %g", (double)err);
103:       PetscCall(VecRestoreArray(vecFull, &arrayFull));
104:       PetscCall(PetscTabulationDestroy(&Tfull));
105:       PetscCall(PetscTabulationDestroy(&Tsub));
106:       /* clean up */
107:       PetscCall(PetscFree(arraySub));
108:       PetscCall(VecDestroy(&vecFull));
109:       PetscCall(PetscRandomDestroy(&rand));
110:       PetscCall(PetscFree4(testSub, testFull, outSub, outFull));
111:       PetscCall(PetscFEDestroy(&traceFE));
112:     }
113:   }
114:   PetscCall(PetscFEDestroy(&fe));
115:   PetscCall(PetscFinalize());
116:   return 0;
117: }

119: /*TEST
120:   test:
121:     suffix: 0
122:     args: -petscspace_degree 1 -trace_fe_view
123: TEST*/