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;
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: PetscInt i, j;
49: PetscSection sectionFull;
50: Vec vecFull;
51: PetscScalar *arrayFull, *arraySub;
52: PetscReal err;
53: PetscRandom rand;
55: PetscCall(DMPlexGetCone(dm, cStart, &cone));
56: point = cone[0];
57: PetscCall(PetscFECreatePointTrace(fe, point, &traceFE));
58: PetscCall(PetscFESetUp(traceFE));
59: PetscCall(PetscFEViewFromOptions(traceFE, NULL, "-trace_fe_view"));
60: PetscCall(PetscMalloc4(dim - 1, &testSub, dim, &testFull, Nc, &outSub, Nc, &outFull));
61: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
62: PetscCall(PetscRandomSetFromOptions(rand));
63: PetscCall(PetscRandomSetInterval(rand, -1., 1.));
64: /* create a random point in the trace domain */
65: for (i = 0; i < dim - 1; i++) PetscCall(PetscRandomGetValueReal(rand, &testSub[i]));
66: PetscCall(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: PetscCall(PetscFEGetDimension(fe, &nFull));
73: PetscCall(VecCreateSeq(PETSC_COMM_SELF, nFull, &vecFull));
74: PetscCall(VecGetArray(vecFull, &arrayFull));
75: for (i = 0; i < nFull; i++) PetscCall(PetscRandomGetValue(rand, &arrayFull[i]));
76: PetscCall(VecRestoreArray(vecFull, &arrayFull));
77: /* create a vector on the trace domain */
78: PetscCall(PetscFEGetDimension(traceFE, &nSub));
79: /* get the subset of the original finite element space that is supported on the trace space */
80: PetscCall(PetscDualSpaceGetSection(dualspace, §ionFull));
81: PetscCall(PetscSectionSetUp(sectionFull));
82: /* get the trace degrees of freedom */
83: PetscCall(PetscMalloc1(nSub, &arraySub));
84: PetscCall(DMPlexVecGetClosure(dm, sectionFull, vecFull, point, &nSub, &arraySub));
85: /* get the tabulations */
86: PetscCall(PetscFECreateTabulation(traceFE, 1, 1, testSub, 0, &Tsub));
87: PetscCall(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: PetscCall(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);
103: PetscCheck(err <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Trace FE error %g", (double)err);
104: PetscCall(VecRestoreArray(vecFull, &arrayFull));
105: PetscCall(PetscTabulationDestroy(&Tfull));
106: PetscCall(PetscTabulationDestroy(&Tsub));
107: /* clean up */
108: PetscCall(PetscFree(arraySub));
109: PetscCall(VecDestroy(&vecFull));
110: PetscCall(PetscRandomDestroy(&rand));
111: PetscCall(PetscFree4(testSub, testFull, outSub, outFull));
112: PetscCall(PetscFEDestroy(&traceFE));
113: }
114: }
115: PetscCall(PetscFEDestroy(&fe));
116: PetscCall(PetscFinalize());
117: return 0;
118: }
120: /*TEST
121: test:
122: suffix: 0
123: args: -petscspace_degree 1 -trace_fe_view
124: TEST*/