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 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