1 static char help[] = "Tests affine subspaces.\n\n"; 2 3 #include <petscfe.h> 4 #include <petscdmplex.h> 5 #include <petscdmshell.h> 6 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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 18 comm = PETSC_COMM_WORLD; 19 PetscOptionsBegin(comm,"","Options for subspace test","none"); 20 PetscCall(PetscOptionsRangeInt("-dim", "The spatial dimension","ex5.c",dim,&dim,NULL,1,3)); 21 PetscCall(PetscOptionsBool("-simplex", "Test simplex element","ex5.c",simplex,&simplex,NULL)); 22 PetscCall(PetscOptionsBoundedInt("-num_comp", "Number of components in space","ex5.c",Nc,&Nc,NULL,1)); 23 PetscOptionsEnd(); 24 PetscCall(DMShellCreate(comm,&dm)); 25 PetscCall(PetscFECreateDefault(comm,dim,Nc,simplex,NULL,PETSC_DEFAULT,&fe)); 26 PetscCall(DMDestroy(&dm)); 27 PetscCall(PetscFESetName(fe, "solution")); 28 PetscCall(PetscFEGetBasisSpace(fe,&space)); 29 PetscCall(PetscSpaceGetNumComponents(space,&Nc)); 30 PetscCall(PetscFEGetDualSpace(fe,&dualspace)); 31 PetscCall(PetscDualSpaceGetHeightSubspace(dualspace,1,&dualsubspace)); 32 PetscCall(PetscDualSpaceGetDM(dualspace,&dm)); 33 PetscCall(DMPlexGetHeightStratum(dm,0,&cStart,&cEnd)); 34 if (cEnd > cStart) { 35 PetscInt coneSize; 36 37 PetscCall(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 PetscTabulation Tsub, Tfull; 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; 53 54 PetscCall(DMPlexGetCone(dm,cStart,&cone)); 55 point = cone[0]; 56 PetscCall(PetscFECreatePointTrace(fe,point,&traceFE)); 57 PetscCall(PetscFESetUp(traceFE)); 58 PetscCall(PetscFEViewFromOptions(traceFE,NULL,"-trace_fe_view")); 59 PetscCall(PetscMalloc4(dim - 1,&testSub,dim,&testFull,Nc,&outSub,Nc,&outFull)); 60 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand)); 61 PetscCall(PetscRandomSetFromOptions(rand)); 62 PetscCall(PetscRandomSetInterval(rand,-1.,1.)); 63 /* create a random point in the trace domain */ 64 for (i = 0; i < dim - 1; i++) { 65 PetscCall(PetscRandomGetValueReal(rand,&testSub[i])); 66 } 67 PetscCall(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 PetscCall(PetscFEGetDimension(fe,&nFull)); 74 PetscCall(VecCreateSeq(PETSC_COMM_SELF,nFull,&vecFull)); 75 PetscCall(VecGetArray(vecFull,&arrayFull)); 76 for (i = 0; i < nFull; i++) { 77 PetscCall(PetscRandomGetValue(rand,&arrayFull[i])); 78 } 79 PetscCall(VecRestoreArray(vecFull,&arrayFull)); 80 /* create a vector on the trace domain */ 81 PetscCall(PetscFEGetDimension(traceFE,&nSub)); 82 /* get the subset of the original finite element space that is supported on the trace space */ 83 PetscCall(PetscDualSpaceGetSection(dualspace,§ionFull)); 84 PetscCall(PetscSectionSetUp(sectionFull)); 85 /* get the trace degrees of freedom */ 86 PetscCall(PetscMalloc1(nSub,&arraySub)); 87 PetscCall(DMPlexVecGetClosure(dm,sectionFull,vecFull,point,&nSub,&arraySub)); 88 /* get the tabulations */ 89 PetscCall(PetscFECreateTabulation(traceFE,1,1,testSub,0,&Tsub)); 90 PetscCall(PetscFECreateTabulation(fe,1,1,testFull,0,&Tfull)); 91 for (i = 0; i < Nc; i++) { 92 outSub[i] = 0.0; 93 for (j = 0; j < nSub; j++) { 94 outSub[i] += Tsub->T[0][j * Nc + i] * arraySub[j]; 95 } 96 } 97 PetscCall(VecGetArray(vecFull,&arrayFull)); 98 err = 0.0; 99 for (i = 0; i < Nc; i++) { 100 PetscScalar diff; 101 102 outFull[i] = 0.0; 103 for (j = 0; j < nFull; j++) { 104 outFull[i] += Tfull->T[0][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 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Trace FE error %g",err); 112 } 113 PetscCall(VecRestoreArray(vecFull,&arrayFull)); 114 PetscCall(PetscTabulationDestroy(&Tfull)); 115 PetscCall(PetscTabulationDestroy(&Tsub)); 116 /* clean up */ 117 PetscCall(PetscFree(arraySub)); 118 PetscCall(VecDestroy(&vecFull)); 119 PetscCall(PetscRandomDestroy(&rand)); 120 PetscCall(PetscFree4(testSub,testFull,outSub,outFull)); 121 PetscCall(PetscFEDestroy(&traceFE)); 122 } 123 } 124 PetscCall(PetscFEDestroy(&fe)); 125 PetscCall(PetscFinalize()); 126 return 0; 127 } 128 129 /*TEST 130 test: 131 suffix: 0 132 args: -petscspace_degree 1 -trace_fe_view 133 TEST*/ 134