Actual source code: ex5.c
petsc-3.10.5 2019-03-28
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);if (ierr) return ierr;
19: comm = PETSC_COMM_WORLD;
20: PetscOptionsBegin(comm,"","Options for subspace test","none");
21: PetscOptionsInt("-dim", "The spatial dimension","ex5.c",dim,&dim,NULL);
22: PetscOptionsBool("-simplex", "Test simplex element","ex5.c",simplex,&simplex,NULL);
23: PetscOptionsInt("-num_comp", "Number of components in space","ex5.c",Nc,&Nc,NULL);
24: PetscOptionsEnd();
25: DMShellCreate(comm,&dm);
26: PetscFECreateDefault(comm,dim,Nc,simplex,NULL,PETSC_DEFAULT,&fe);
27: DMDestroy(&dm);
28: PetscFEGetBasisSpace(fe,&space);
29: PetscSpaceGetNumComponents(space,&Nc);
30: PetscFEGetDualSpace(fe,&dualspace);
31: PetscDualSpaceGetHeightSubspace(dualspace,1,&dualsubspace);
32: PetscDualSpaceGetDM(dualspace,&dm);
33: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
34: if (cEnd > cStart) {
35: PetscInt coneSize;
37: DMPlexGetConeSize(dm,cStart,&coneSize);
38: if (coneSize) {
39: PetscFE traceFE;
40: const PetscInt *cone;
41: PetscInt point, nSub, nFull;
42: PetscReal xi0[3] = {-1., -1., -1.};
43: PetscScalar *outSub, *outFull;
44: PetscReal *testSub, *testFull;
45: PetscReal *Bsub, *Bfull;
46: PetscReal J[9], detJ;
47: PetscInt i, j;
48: PetscSection sectionFull;
49: Vec vecFull;
50: PetscScalar *arrayFull, *arraySub;
51: PetscReal err;
52: PetscRandom rand;
54: DMPlexGetCone(dm,cStart,&cone);
55: point = cone[0];
56: PetscFECreatePointTrace(fe,point,&traceFE);
57: PetscFESetUp(traceFE);
58: PetscFEViewFromOptions(traceFE,NULL,"-trace_fe_view");
59: PetscMalloc4(dim - 1,&testSub,dim,&testFull,Nc,&outSub,Nc,&outFull);
60: PetscRandomCreate(PETSC_COMM_SELF,&rand);
61: PetscRandomSetFromOptions(rand);
62: PetscRandomSetInterval(rand,-1.,1.);
63: /* create a random point in the trace domain */
64: for (i = 0; i < dim - 1; i++) {
65: PetscRandomGetValueReal(rand,&testSub[i]);
66: }
67: DMPlexComputeCellGeometryFEM(dm,point,NULL,testFull,J,NULL,&detJ);
68: /* project it into the full domain */
69: for (i = 0; i < dim; i++) {
70: for (j = 0; j < dim - 1; j++) testFull[i] += J[i * dim + j] * (testSub[j] - xi0[j]);
71: }
72: /* create a random vector in the full domain */
73: PetscFEGetDimension(fe,&nFull);
74: VecCreateSeq(PETSC_COMM_SELF,nFull,&vecFull);
75: VecGetArray(vecFull,&arrayFull);
76: for (i = 0; i < nFull; i++) {
77: PetscRandomGetValue(rand,&arrayFull[i]);
78: }
79: VecRestoreArray(vecFull,&arrayFull);
80: /* create a vector on the trace domain */
81: PetscFEGetDimension(traceFE,&nSub);
82: /* get the subset of the original finite element space that is supported on the trace space */
83: PetscDualSpaceCreateSection(dualspace,§ionFull);
84: PetscSectionSetUp(sectionFull);
85: /* get the trace degrees of freedom */
86: PetscMalloc1(nSub,&arraySub);
87: DMPlexVecGetClosure(dm,sectionFull,vecFull,point,&nSub,&arraySub);
88: /* get the tabulations */
89: PetscFEGetTabulation(traceFE,1,testSub,&Bsub,NULL,NULL);
90: PetscFEGetTabulation(fe,1,testFull,&Bfull,NULL,NULL);
91: for (i = 0; i < Nc; i++) {
92: outSub[i] = 0.0;
93: for (j = 0; j < nSub; j++) {
94: outSub[i] += Bsub[j * Nc + i] * arraySub[j];
95: }
96: }
97: VecGetArray(vecFull,&arrayFull);
98: err = 0.0;
99: for (i = 0; i < Nc; i++) {
100: PetscScalar diff;
102: outFull[i] = 0.0;
103: for (j = 0; j < nFull; j++) {
104: outFull[i] += Bfull[j * Nc + i] * arrayFull[j];
105: }
106: diff = outFull[i] - outSub[i];
107: err += PetscRealPart(PetscConj(diff) * diff);
108: }
109: err = PetscSqrtReal(err);
110: if (err > PETSC_SMALL) {
111: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Trace FE error %g\n",err);
112: }
113: VecRestoreArray(vecFull,&arrayFull);
114: PetscFERestoreTabulation(fe,1,testFull,&Bfull,NULL,NULL);
115: PetscFERestoreTabulation(traceFE,1,testSub,&Bsub,NULL,NULL);
116: /* clean up */
117: PetscFree(arraySub);
118: PetscSectionDestroy(§ionFull);
119: VecDestroy(&vecFull);
120: PetscRandomDestroy(&rand);
121: PetscFree4(testSub,testFull,outSub,outFull);
122: PetscFEDestroy(&traceFE);
123: }
124: }
125: PetscFEDestroy(&fe);
126: PetscFinalize();
127: return ierr;
128: }
130: /*TEST
131: test:
132: suffix: 0
133: args: -petscspace_degree 1 -trace_fe_view
134: TEST*/