1 static char help[] = "Tests affine subspaces.\n\n";
2
3 #include <petscfe.h>
4 #include <petscdmplex.h>
5 #include <petscdmshell.h>
6
main(int argc,char ** argv)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;
16
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;
37
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;
54
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;
96
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 }
119
120 /*TEST
121 test:
122 suffix: 0
123 args: -petscspace_degree 1 -trace_fe_view
124 TEST*/
125