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, §ionFull);
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*/