Actual source code: ex5.c

  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);
 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++) PetscRandomGetValueReal(rand, &testSub[i]);
 66:       DMPlexComputeCellGeometryFEM(dm, point, NULL, testFull, J, NULL, &detJ);
 67:       /* project it into the full domain */
 68:       for (i = 0; i < dim; i++) {
 69:         for (j = 0; j < dim - 1; j++) testFull[i] += J[i * dim + j] * (testSub[j] - xi0[j]);
 70:       }
 71:       /* create a random vector in the full domain */
 72:       PetscFEGetDimension(fe, &nFull);
 73:       VecCreateSeq(PETSC_COMM_SELF, nFull, &vecFull);
 74:       VecGetArray(vecFull, &arrayFull);
 75:       for (i = 0; i < nFull; i++) PetscRandomGetValue(rand, &arrayFull[i]);
 76:       VecRestoreArray(vecFull, &arrayFull);
 77:       /* create a vector on the trace domain */
 78:       PetscFEGetDimension(traceFE, &nSub);
 79:       /* get the subset of the original finite element space that is supported on the trace space */
 80:       PetscDualSpaceGetSection(dualspace, &sectionFull);
 81:       PetscSectionSetUp(sectionFull);
 82:       /* get the trace degrees of freedom */
 83:       PetscMalloc1(nSub, &arraySub);
 84:       DMPlexVecGetClosure(dm, sectionFull, vecFull, point, &nSub, &arraySub);
 85:       /* get the tabulations */
 86:       PetscFECreateTabulation(traceFE, 1, 1, testSub, 0, &Tsub);
 87:       PetscFECreateTabulation(fe, 1, 1, testFull, 0, &Tfull);
 88:       for (i = 0; i < Nc; i++) {
 89:         outSub[i] = 0.0;
 90:         for (j = 0; j < nSub; j++) outSub[i] += Tsub->T[0][j * Nc + i] * arraySub[j];
 91:       }
 92:       VecGetArray(vecFull, &arrayFull);
 93:       err = 0.0;
 94:       for (i = 0; i < Nc; i++) {
 95:         PetscScalar diff;

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

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