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 PetscErrorCode ierr; 17 18 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 19 comm = PETSC_COMM_WORLD; 20 ierr = PetscOptionsBegin(comm,"","Options for subspace test","none");PetscCall(ierr); 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 ierr = PetscOptionsEnd();PetscCall(ierr); 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++) { 66 PetscCall(PetscRandomGetValueReal(rand,&testSub[i])); 67 } 68 PetscCall(DMPlexComputeCellGeometryFEM(dm,point,NULL,testFull,J,NULL,&detJ)); 69 /* project it into the full domain */ 70 for (i = 0; i < dim; i++) { 71 for (j = 0; j < dim - 1; j++) testFull[i] += J[i * dim + j] * (testSub[j] - xi0[j]); 72 } 73 /* create a random vector in the full domain */ 74 PetscCall(PetscFEGetDimension(fe,&nFull)); 75 PetscCall(VecCreateSeq(PETSC_COMM_SELF,nFull,&vecFull)); 76 PetscCall(VecGetArray(vecFull,&arrayFull)); 77 for (i = 0; i < nFull; i++) { 78 PetscCall(PetscRandomGetValue(rand,&arrayFull[i])); 79 } 80 PetscCall(VecRestoreArray(vecFull,&arrayFull)); 81 /* create a vector on the trace domain */ 82 PetscCall(PetscFEGetDimension(traceFE,&nSub)); 83 /* get the subset of the original finite element space that is supported on the trace space */ 84 PetscCall(PetscDualSpaceGetSection(dualspace,§ionFull)); 85 PetscCall(PetscSectionSetUp(sectionFull)); 86 /* get the trace degrees of freedom */ 87 PetscCall(PetscMalloc1(nSub,&arraySub)); 88 PetscCall(DMPlexVecGetClosure(dm,sectionFull,vecFull,point,&nSub,&arraySub)); 89 /* get the tabulations */ 90 PetscCall(PetscFECreateTabulation(traceFE,1,1,testSub,0,&Tsub)); 91 PetscCall(PetscFECreateTabulation(fe,1,1,testFull,0,&Tfull)); 92 for (i = 0; i < Nc; i++) { 93 outSub[i] = 0.0; 94 for (j = 0; j < nSub; j++) { 95 outSub[i] += Tsub->T[0][j * Nc + i] * arraySub[j]; 96 } 97 } 98 PetscCall(VecGetArray(vecFull,&arrayFull)); 99 err = 0.0; 100 for (i = 0; i < Nc; i++) { 101 PetscScalar diff; 102 103 outFull[i] = 0.0; 104 for (j = 0; j < nFull; j++) { 105 outFull[i] += Tfull->T[0][j * Nc + i] * arrayFull[j]; 106 } 107 diff = outFull[i] - outSub[i]; 108 err += PetscRealPart(PetscConj(diff) * diff); 109 } 110 err = PetscSqrtReal(err); 111 if (err > PETSC_SMALL) { 112 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Trace FE error %g",err); 113 } 114 PetscCall(VecRestoreArray(vecFull,&arrayFull)); 115 PetscCall(PetscTabulationDestroy(&Tfull)); 116 PetscCall(PetscTabulationDestroy(&Tsub)); 117 /* clean up */ 118 PetscCall(PetscFree(arraySub)); 119 PetscCall(VecDestroy(&vecFull)); 120 PetscCall(PetscRandomDestroy(&rand)); 121 PetscCall(PetscFree4(testSub,testFull,outSub,outFull)); 122 PetscCall(PetscFEDestroy(&traceFE)); 123 } 124 } 125 PetscCall(PetscFEDestroy(&fe)); 126 PetscCall(PetscFinalize()); 127 return 0; 128 } 129 130 /*TEST 131 test: 132 suffix: 0 133 args: -petscspace_degree 1 -trace_fe_view 134 TEST*/ 135