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 DM dm; 9 PetscFE fe; 10 PetscSpace space; 11 PetscDualSpace dualspace, dualsubspace; 12 PetscInt dim = 2, Nc = 3, cStart, cEnd; 13 PetscBool simplex = PETSC_TRUE; 14 MPI_Comm comm; 15 16 PetscFunctionBeginUser; 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++) { PetscCall(PetscRandomGetValueReal(rand, &testSub[i])); } 65 PetscCall(DMPlexComputeCellGeometryFEM(dm, point, NULL, testFull, J, NULL, &detJ)); 66 /* project it into the full domain */ 67 for (i = 0; i < dim; i++) { 68 for (j = 0; j < dim - 1; j++) testFull[i] += J[i * dim + j] * (testSub[j] - xi0[j]); 69 } 70 /* create a random vector in the full domain */ 71 PetscCall(PetscFEGetDimension(fe, &nFull)); 72 PetscCall(VecCreateSeq(PETSC_COMM_SELF, nFull, &vecFull)); 73 PetscCall(VecGetArray(vecFull, &arrayFull)); 74 for (i = 0; i < nFull; i++) { PetscCall(PetscRandomGetValue(rand, &arrayFull[i])); } 75 PetscCall(VecRestoreArray(vecFull, &arrayFull)); 76 /* create a vector on the trace domain */ 77 PetscCall(PetscFEGetDimension(traceFE, &nSub)); 78 /* get the subset of the original finite element space that is supported on the trace space */ 79 PetscCall(PetscDualSpaceGetSection(dualspace, §ionFull)); 80 PetscCall(PetscSectionSetUp(sectionFull)); 81 /* get the trace degrees of freedom */ 82 PetscCall(PetscMalloc1(nSub, &arraySub)); 83 PetscCall(DMPlexVecGetClosure(dm, sectionFull, vecFull, point, &nSub, &arraySub)); 84 /* get the tabulations */ 85 PetscCall(PetscFECreateTabulation(traceFE, 1, 1, testSub, 0, &Tsub)); 86 PetscCall(PetscFECreateTabulation(fe, 1, 1, testFull, 0, &Tfull)); 87 for (i = 0; i < Nc; i++) { 88 outSub[i] = 0.0; 89 for (j = 0; j < nSub; j++) { outSub[i] += Tsub->T[0][j * Nc + i] * arraySub[j]; } 90 } 91 PetscCall(VecGetArray(vecFull, &arrayFull)); 92 err = 0.0; 93 for (i = 0; i < Nc; i++) { 94 PetscScalar diff; 95 96 outFull[i] = 0.0; 97 for (j = 0; j < nFull; j++) { outFull[i] += Tfull->T[0][j * Nc + i] * arrayFull[j]; } 98 diff = outFull[i] - outSub[i]; 99 err += PetscRealPart(PetscConj(diff) * diff); 100 } 101 err = PetscSqrtReal(err); 102 PetscCheck(err <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Trace FE error %g", (double)err); 103 PetscCall(VecRestoreArray(vecFull, &arrayFull)); 104 PetscCall(PetscTabulationDestroy(&Tfull)); 105 PetscCall(PetscTabulationDestroy(&Tsub)); 106 /* clean up */ 107 PetscCall(PetscFree(arraySub)); 108 PetscCall(VecDestroy(&vecFull)); 109 PetscCall(PetscRandomDestroy(&rand)); 110 PetscCall(PetscFree4(testSub, testFull, outSub, outFull)); 111 PetscCall(PetscFEDestroy(&traceFE)); 112 } 113 } 114 PetscCall(PetscFEDestroy(&fe)); 115 PetscCall(PetscFinalize()); 116 return 0; 117 } 118 119 /*TEST 120 test: 121 suffix: 0 122 args: -petscspace_degree 1 -trace_fe_view 123 TEST*/ 124