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